comphy (version 1.0.5)

RK4ODE: Runge-Kutta 4th order method for systems of ODEs

Description

Solves a system of \(m\) first-order ODEs using the classical fourth-order Runge-Kutta method.

Usage

RK4ODE(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 RK4ODE. 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

This method achieves high accuracy by evaluating the gradient function four times per step. It has local error \(O(h^5)\) and global error \(O(h^4)\). It is one of the most widely used methods for solving initial value problems numerically.

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