Contents Previous Next Subchapters

Solving Stiff Ordinary Differential Equations
Syntax odestiff(function fvalordertitout, ...
         
maxstepminstepyiepsskip)
See Also ode4rk , odepade

Description
Uses a backward difference method to solve the differential equation
     d
     -- y(t) = f(ty(t))
     dt
with initial condition y(t) = yi at t = ti. The integer scalar order specifies the order of the backward difference method. The real or double-precision scalar ti specifies the initial value for t. The row vector tout has the same type as ti and specifies the values of t at which odestiff should calculate the value of y. The scalar maxstep has the same type as ti and specifies the maximum step size in t. The scalar minstep has the same type as ti and specifies the minimum step size in t. The column vector yi has the same type as ti and specifies the initial value of y. The column vector eps has the same type and dimension as yi and eps(i) is the desired accuracy for the i-th element of y(t). The integer scalar skip specifies the level of tracing inside of odestiff. If skip = 0, no tracing is done; otherwise, the value of y(t) is printed at time tout(skip), tout(2 skip) ... . The return value for odestiff has the same type and row dimension as yi, the same number of columns as tout, and its j-th column is the value of y at time tout(j).

fval(tinyinfoutdfout)
This function call sets fout the value of f(tinyin) and dfout to the value
     d
     -- f(tinyin)
     dy
The scalar tin has the same type as ti. The vector yin has the same type and dimension as yi. The input value of fout does not matter. Its output value has the same type and dimension and row dimension as yi. The input value of dfout does not matter. Its output value is a square matrix with the same type and row dimension as yi and its (i,j)-th element is the partial of the i-th element of f(ty) with respect to the j-th element of y.

Example
The following 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 fval(tin, yin, fout, dfout) begin
     fout   = {yin(2), -yin(1)}
     dfout  = {[0., 1.], [-1., 0.]}
     return
end
order   = 4
ti      = 0.
tf      = 2 * pi
maxstep = 1.
minstep = .01
yi      = {0., 1.}
eps     = {1e-4, 1e-4}
skip    = 0
tout    = seq(30)' * (tf - ti) / real(30)
yout    = odestiff(function fval, order, ti, tout, ...
                   maxstep, minstep, yi, eps, skip)
gplot(tout', yout')