Contents Previous Next Subchapters

A Fourth Order Runge-Kutta ODE Solver
Syntax ode4rk(function ftitfyin)
See Also ode45rk , odepade , odestiff

Uses a fourth-order Runge-Kutta method to solve the differential equation
     -- y(t) = f(ty(t))
with initial condition y(ti) = yi. The return value is the numerical approximation for
     [y(ti + s), y(ti + 2 s), . . . , y(ti + n s)], 
where s = (tf - ti)/n is the step size of the Runge-Kutta method. The return value has the same type and number of rows as yi, and has n columns.

The scalar ti specifies the initial time point. If yi is complex, ti is double-precision; otherwise, ti has the same type as yi. The scalar tf has the same type as ti and specifies the final time point. The column vector yi is either real, double-precision, or complex and specifies the value of y(t) at t = ti. The integer scalar n specifies the number of integration steps between ti and tf.

The function call f(ty) returns the derivative of y at time t, provided t is a scalar with the same type as ti, and y is a vector with the same type and dimension as yi. The return value has the same type and dimension as yi.

The follows program plots the solution of the differential equation

     y'(t) =  y (t)     y (0) = 0
      1        2         1
     y'(t) = -y (t)     y (0) = 1
      2        1         2
The solution is

     y (t) = sin(t)

     y (t) = cos(t)

function f(t, y) begin
     return {y(2), -y(1)}
ti = 0.
tf = 2 * pi
yi = {0., 1.}
n  = 30
y  = ode4rk(function f, ti, tf, yi, n)
s  = (tf - ti) / real(n)
t  = tf + seq(n) * s
gplot(t, y')