Contents Previous Next Subchapters Current Chapters-> trapz gaussq quad2d gaussq2d gaussleg quadint int2inf ode4rk expmat ode45rk odepade odestiff Parent Chapters-> Omatrix6 integration odepade Search Tools-> contents reference index search

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') ```