comphy (version 1.0.5)

HeunODE: Heun method for systems of ODEs

Description

Solves a system of \(m\) first-order ODEs using the Heun method (also known as the improved Euler method).

Usage

HeunODE(f, t0, tf, y0, h, ...)

Value

A list with elements t (time points) and y

(solution matrix). The first row of the matrix contains the initial values of y at time t0. Each column of the matrix contains the numerical solution for each one of the \(m\)

functions of the system of ODEs.

Arguments

f

A function of the form f(t,y) returning a numeric vector. It must be defined before using HeunODE. This function is the right hand side of the ODE, i.e. the gradient of the ODE system.

t0

Initial time.

tf

Final time.

y0

A numeric vector with initial values (length = \(m\)).

h

Step size.

...

Other parameters potentially needed by the gradient function.

Details

The method improves upon the Euler method by using an average of the slopes at the beginning and end of each time step. It is more accurate, with local error \(O(h^3)\) and global error \(O(h^2)\).

Examples

# IVP: \eqn{dy/dt=6-2y,\ y(0)=0}.
# Define gradient
f <- function(t,y) {dy <- 6-2*y; return(dy)}

# Solution interval
t0 <- 0
tf <- 2

# Initial condition
y0 <- 0

# Step
h <- 0.1

# Numerical solution
ltmp <- HeunODE(f,t0,tf,y0,h)

# Print grid
print(ltmp$t)

# Print numerical solution
print(ltmp$y)

# Example with two ODEs. 
# \eqn{dy_1/dt=y_1+2y_2}
# \eqn{dy_2/dt=(3/2)y_1-y_2}
# \eqn{y_1(0)=1, y_2(0)=-2}

# Define gradient
dy <- function(t,y) {
  dy1 <- y[1]+2*y[2] 
  dy2 <- 1.5*y[1]-y[2] 
  return(c(dy1,dy2))
}

# Solution interval
t0 <- 0
tf <- 2

# Initial conditions
y0 <- c(1,-2)

# Step
h <- 0.1

# Numerical solution
ltmp <- HeunODE(dy,t0,tf,y0,h)

# Print grid
print(ltmp$t)

# Print numerical solution y1
print(ltmp$y[,1])

# Print numerical solution y2
print(ltmp$y[,2])