comphy (version 1.0.5)

EulerODE: Euler method for systems of ODEs

Description

Solves a system of \(m\) first-order ODEs using the explicit Euler method.

Usage

EulerODE(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 EulerODE. 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 is accurate and stable when the stepsize h is relatively small. The local error is \(O(h^2)\), while the global error is \(O(h)\). Other numerical methods are generally used to calculate solutions with a higher accuracy.

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 <- EulerODE(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 <- EulerODE(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])