Contents Previous Next Subchapters

Fifth Order Runge-Kutta Method With Fourth Order Error Control
, ,
Syntax [toutyoutsout] = ode45rk(function f, ...
t0tfy0tollevelsminsmaxsini)
Optional souttollevelsminsmaxsini
See Also ode23 , ode45 , ode4rk , odepade , odestiff

Description
Solves a system the following system of ordinary differential equations:
     dy
     -- = f(ty(t))  and  y(t0) = y0
     dt
using a fifth order Runge-Kutta method and using a fourth order Runge-Kutta method to estimate the error in the solution reference

f(ty)
This function call returns the value of f(ty(t)) where t is an integer, real or double-precision scalar, and y is an integer, real or double-precision column vector with the same length as y0. The return value is an integer, real, double-precision or complex vector with the same dimension as y0.

t0
is an integer, real or double-precision scalar specifying the initial value of t.

tf
is an integer, real or double-precision scalar specifying the final value of t.

y0
is an integer, real, double-precision or complex column vector specifying the value of y(t0).

tol
is an integer, real or double-precision vector that specifies the requested accuracy in the return values in yout. If tol is a column vector with the same length as y0, tol(j) specifies the accuracy for the elements in the j-th column of yout. All the elements of tol must be greater than zero. If tol is a scalar, it is equivalent to the case where it is a column vector with the same length as y0 and all the elements have the value specified by the scalar. If tol is not present in the syntax, the default value
     tol  = 1d-3
is used.

level
is a logical integer, real or double-precision scalar that specifies the level of tracing during solution of the differential equation. If level is equal to zero, no tracing is done. If level is equal to one, the current value of t, the current step size, and the transpose of the current value of y(t) are printed upon completion of each successful step. If level is not present in the syntax, the default value
     level  = 0
is used.

smin
is an integer, real or double-precision scalar that specifies the minimum step size; i.e., difference between successive elements of tout. If this argument is not present in the syntax, the default value
     smin  = (tf - t0) / 1d4
is used.

smax
is an integer, real or double-precision scalar that specifies the maximum step size; i.e., difference between successive elements of tout. If this argument is not present in the syntax, the default value
     smax  = (tf - t0) / 16
is used.

sini
is an integer, real or double-precision scalar that specifies the initial step size. It must satisfy the relations
     smin < sini < smax
If this argument is not present in the syntax, the default value
     sini  = smax
is used.

tout
is a monotone increasing column vector with its first element equal to t0, and its last element equal to tf. In addition,
     smin < tout   -  tout  > smax
                  i+1       i
If all the input values are real, tout is real, otherwise it is double-precision.

yout
is a matrix containing an approximate solution of the differential equation. Its row dimension is equal to the row dimension of tout and its column dimension is equal to the row dimension of y0. The ode45rk routine will try to accomplish the requested accuracy condition which is that for i between one and the row dimension of tout and j between one and the row dimension of y0,
     | yout  - y (tout(i)) | <= tol
           j    j                  j
If y0 is complex, yout is also complex. Otherwise yout is double-precision.

sout
is a scalar that suggests a value for sini when continuing solution of the differential equation using ode45rk. (A continuation is where the new value of t0 is equal to the old value of tf and the new value of y0 is approximately equal to y(tf).)

Example
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)
      1

     y (t) = cos(t)
      2


clear
#
function f(t, y) begin
     return {y(2), -y(1)}
end
t0           = 0.
tf           = 2. * pi
y0           = {0., 1.}
[tout, yout] =  ode45rk(function f, t0, tf, y0)
gplot(tout, yout)

Mlmode
In Mlmode , this function is accessed using the following syntax
     [
toutyout] = ode(ypt0tfy0tollevel)
where ode is ode23 or ode45 and the arguments tol, level are optional. In addition, each call of the form f(ty) in ode45rk is translated to a call of the form feval(ypty).

Suppose that the file yp.m in the current directory contains the text
     function dy = yp(t, y) 
     dy = [y(2); -y(1)];
If in Mlmode you execute the program below, the same plot will be generated as for the ode45rk example above.

clear
t0           = 0.;
tf           = 2. * pi;
y0           = [0.; 1.];
[tout, yout] =  ode45('yp', t0, tf, y0)
plot(tout, yout)

Reference
The method uses the Cash-Karp parameters in Section 16.2 of the second edition of Numerical Recipes in Fortran by William H. Press, Saul A. Teukolsky, William T. Vetterling, and Brian P. Flannery.