Contents Previous Next Subchapters

Solving An ODE Using The Matrix Exponential
Syntax odepade(A, tf, yi, n)
See Also ode4rk , odestiff

Description
Solves the differential equation dy/dt = A y with initial condition y(0) = yi, where A is a real, double-precision, or complex square matrix specifying the linear differential equation, tf is a scalar specifying the final time point (the initial time point is automatically 0), yi is a column vector specifying the value of y(t) at t = 0, and n is an integer scalar specifying the number of time points at which to return the value of y(t). If A is complex, tf must be double-precision; otherwise, tf must have the same type as A. The matrix A and the vector yi have the same number of rows.

The return value is the matrix
     [y(s), y(2 s), . . . , y(n s)]
that is, the j-th column of the return value is y(t) at t = j s, where s = tf / n. The return value has the same type and number of rows as yi, and has n columns.

This function uses the following Pade approximation for the matrix exponential of B, where B = s A
                    2       3   -1                2       3
     [       B     B       B   ]   [       B     B       B   ]
     [ I  -  -  +  --  -  ---  ]   [ I  +  -  +  --  +  ---  ]
     [       2     10     120  ]   [       2     10     120  ]
where I is the identity matrix.

Example
The following program plots the solution of the differential equation
     d  [ y ]   [ 0  1 ]   [ y ]
     -- [  1] = [      ] * [  1]
     dt [ y ]   [-1  0 ]   [ y ]
           2                  2
with the initial condition

     y (0) = 0
      1
     y (0) = 1
      2
The solution is

     y (t) = sin(t)
      1
     y (t) = cos(t)
      2

A  = {[0., 1.], [-1., 0.]}  
tf = 2 * pi
yi = {0., 1.}
n  = 30
y  = odepade(A, tf, yi, n)
s  = tf / real(n)
t  = seq(n) * s
gplot(t, y')