ODE Call

CALL ODE (r, "dername", c, t, h <*>, J="jacobian" <*>, EPS=eps <*>, "SAS-data-set" ) ;

The ODE subroutine performs numerical integration of first-order vector differential equations of the form

\[  \frac{d\mb {y}}{dt} = f(t,\mb {y}(t)) ~ ~ ~  \mbox{ with } ~  \mb {y}(0) = \mb {c}  \]

The ODE subroutine returns the following values:

r

is a numeric matrix that contains the results of the integration over connected subintervals. The number of columns in r is equal to the number of subintervals of integration as defined by the argument t. In case of any error in the integration on any subinterval, partial results are not reported in r.

The input arguments to the ODE subroutine are as follows:

dername

specifies the name of a module used to evaluate the integrand.

c

specifies an initial value vector for the variable $y$.

t

specifies a sorted vector that describes the limits of integration over connected subintervals. The simplest form of the vector t contains only the limits of the integration on one interval. The first component of t should contain the initial value, and the second component should be the final value of the independent variable. For more advanced usage of the ODE subroutine, the vector t can contain more than two components. The components of the vector must be sorted in ascending order. Two consecutive components of the vector t are interpreted as a subinterval. The ODE subroutine reports the final result of integration at the right endpoint of each subinterval. This information is vital if $f(\cdot )$ has internal points of discontinuity. To produce accurate solutions, it is essential that you provide the location of these points in the variable t. The continuity of the forcing function is vital to the internal control of error.

h

specifies a numeric vector that contains three components: the minimum allowable step size, hmin; the maximum allowable step size, hmax; and the initial step size to start the integration process, hinit.

jacobian

optionally specifies the name of a module that is used to evaluate the Jacobian analytically. The Jacobian is the matrix $J$, with

\[  J_{ij} = \frac{\partial f_ i}{\partial y_ j}  \]

If the jacobian module is not specified, the ODE subroutine uses a finite-difference method to approximate the Jacobian. The keyword for this option is J.

eps

specifies a scalar that indicates the required accuracy. It has a default value of 1E$-$4. The keyword for this option is EPS.

SAS-data-set

is an optional argument that specifies the name of a valid predefined SAS data set name. The data set is used to save the successful independent and dependent variables of the integration at each step. The keyword for this option is DATA.

The ODE subroutine is an adaptive, variable-order, variable-step-size, stiff integrator based on implicit backward-difference methods. See Aiken (1985), Bickart and Picel (1973), Donelson and Hansen (1971), Gaffney (1984), and Shampine (1978). The integrator is an implicit predictor-corrector method that locally attempts to maintain the prescribed precision eps relative to

\[  d = \max _{0 \leq t \leq T}(\|  y(t) \| _{\infty }, 1)  \]

As you can see from the expression, this quantity is dynamically updated during the integration process and can help you to understand the validity of the results reported by the subroutine.

A Linear Differential Equation

Consider the differential equation

\[  \frac{dy}{dt} = -ty \mbox{ with } y=0.5 \mbox{ at } t = 0  \]

The following statements attempt to find the solution at $t=1$:

/* Define the integrand */
start fun(t,y);
   v = -t*y;
   return(v);
finish;

/* Call ODE */
c   = 0.5;
t   = {0 1};
h   = {1E-12 1 1E-5};
call ode(r1, "FUN", c, t, h);
print r1[format=E21.14];

Figure 24.246: Solution to a Differential Equation at $t=1$

r1
3.03432290135600E-01


In this case, the integration is carried out over $(0, 1)$ to give the value of $y$ at $t=1$. The optional parameter eps has not been specified, so it is internally set to 1E$-$4. Also, the optional parameter jacobian has not been specified, so finite-difference methods are used to estimate the Jacobian. The accuracy of the answer can be increased by specifying eps. For example, set EPS=1E$-$7, as follows:

call ode(r2, "FUN", c, t, h) eps=1E-7;
print r2[format =E21.14];

Figure 24.247: A Solution with Increased Accuracy

r2
3.03265687354960E-01


Compare this value to $0.5e^{-0.5} = 3.03265329856310$E$-01$ and observe that the result is correct through the sixth decimal digit and has an error relative to 1 that is $O(1$E$-7)$.

If the solution was desired at 1 and 2 with an accuracy of 1E$-$7, you would use the following statements:

t   = {0 1 2};
h   = {1E-12 1 1E-5};
call ode(r3, "FUN", c, t, h) eps=1E-7;
print r3[format=E21.14];

Figure 24.248: A Solution at Two Times

r3
3.03265687354960E-01 6.76677185425360E-02


Note that r3 contains the solution at $t=1$ in the first column and at $t=2$ in the second column.

A Discontinuous Forcing Function

Now consider the smoothness of the forcing function $f(\cdot )$. For the purpose of estimating errors, adaptive methods require some degree of smoothness in the function $f(\cdot )$. If this smoothness is not present in $f(\cdot )$ over the interior and including the left endpoint of the subinterval, the reported result does not have the desired accuracy. The function $f(\cdot )$ must be at least continuous. If the function does not meet this requirement, you should specify the discontinuity as an intermediate point. For example, consider the differential equation

\[  \frac{dy}{dt} = \left\{  \begin{array}{ll} t &  \mbox{if } t<1 \\ 0.5t^2 &  \mbox{if } t \geq 1 \end{array} \right.  \]

To find the solution at $t=2$, use the following statements:

/* Define the integrand */
start fun(t,y);
    if t < 1 then v = t;
    else v = .5*t*t;
    return(v);
finish;

c   = 0;
t   = {0 2};
h   = {1E-12 1. 1E-5};
call ode(r1, "FUN", c, t, h) eps=1E-12;
print r1[format=E21.14];

Figure 24.249: Numerical Solution Across a Discontinuity

r1
1.66666626639430E+00


In the preceding case, the integration is carried out over a single interval, $(0,2)$. The optional parameter eps is specified to be 1E$-$12. The optional parameter jacobian is not specified, so finite-difference methods are used to estimate the Jacobian.

Note that the value of r1 does not have the required accuracy (it should contain a 12 decimal-place representation of 5/3), although no error message is produced. The reason is that the function is not continuous at the point $t=1$. Even the lowest-order method cannot produce a local reliable error estimate near the point of discontinuity. To avoid this problem, you can create subintervals so that the integration is carried out first over $(0,1)$ and then over $(1,2)$. The following statements implement this method:

c   = 0;
t   = {0 1 2};
h   = {1E-12 1 1E-5};
call ode(r2, "FUN", c, t, h) eps=1E-12;
print r2[format=E21.14];

Figure 24.250: Numerical Solution on Subintervals

r2
5.00000000003280E-01 1.66666666667280E+00


The variable r2 contains the solutions at both $t=1$ and $t=2$, and the errors are of the specified order. Although there is no interest in the solution at the point $t=1$, the advantage of specifying subintervals with no discontinuities is that the function $f(\cdot )$ is infinitely differentiable in each subinterval.

A Piecewise Continuous Forcing Function

When $f(\cdot )$ is continuous, the ODE subroutine can compute the integration to the specified precision, even if the function is defined piecewise. Consider the differential equation

\[  \frac{dy}{dt} = \left\{  \begin{array}{ll} t &  \mbox{if } t<1 \\ t^2 &  \mbox{if } t \geq 1 \end{array} \right.  \]

The following statements finds the solution at $t=2$. Since the function $f(\cdot )$ is continuous, the requirements for error control are satisfied.

/* Define the integrand */
start fun(t,y);
    if t < 1 then v = t;
    else v = t*t;
    return(v);
finish;

c   = 0.5;
t   = {0 2};
h   = {1E-12 1 1E-5};
call ode(r, "FUN", c, t, h) eps=1E-12;
print r[format=E21.14];

Figure 24.251: Numerical Solution Across a Discontinuity

r
3.33333333334290E+00


Comparing Numerical Integration with an Eigenvalue Decomposition

This example compares the ODE subroutine to an eigenvalue decomposition for stiff-linear systems. In the problem

\[  \frac{dy}{dt} = A y \;  \mbox{ with } \;  y(0)=c  \]

where $A$ is a symmetric constant matrix, the solution can be written in terms of the eigenvalue decomposition as

\[  y(t) = U e^{D t} U^{\prime } c  \]

where $U$ is the matrix of eigenvectors and $D$ is a diagonal matrix with the eigenvalues on its diagonal.

The following statements produce two solutions, one by using the ODE subroutine and the other by using the eigenvalue decomposition:

/* Define the integrand */
start fun(t,x) global(a,count);
   count = count+1;
   v = a*x;
   return(v);
finish;

/* Define the Jacobian */
start jac(t,x) global(a);
   return(a);
finish;

a = {-1000 -1 -2 -3,
        -1 -2  3 -1,
        -2  3 -4 -3,
        -3 -1 -3 -5 };

count = 0;
t = {0 1 2};
h = {1E-12  1 1E-5};
eps = 1E-9;
c = {1, 0, 0, 0 };
call ode(z, "FUN", c, t, h)  eps=eps j="JAC";
print z[format=E21.14];
print count;

Figure 24.252: Numerical Integration of a Linear System

z
1.85787365492010E-06 6.58581431443360E-06
-1.76251618648210E-03 -6.35540480231360E-03
-1.56685608329260E-03 -5.72429355490000E-03
1.01207978491490E-03 3.73655984699120E-03

count
437


/* Do the eigenvalue decomposition */
start eval(t) global(d,u,c);
   v = u*diag(exp(d*t))*u`*c;
   return(v);
finish;

call eigen(d,u,a);
free z1;
do i = 1 to nrow(t)*ncol(t)-1;
   z1 = z1 || (eval(t[i+1]));
end;
print z1[format=E21.14];

Figure 24.253: Analytic Solution of a Linear System

z
1.85787365492010E-06 6.58581431443360E-06
-1.76251618648210E-03 -6.35540480231360E-03
-1.56685608329260E-03 -5.72429355490000E-03
1.01207978491490E-03 3.73655984699120E-03

count
437

z1
1.85787839378720E-06 6.58580950202810E-06
-1.76251639451200E-03 -6.35540294085790E-03
-1.56685625917120E-03 -5.72429205508220E-03
1.01207878768800E-03 3.73655890904620E-03


Is this an $O(1$E$-9)$ result? Note that for the problem

\[  d = \max _{0 \leq t \leq T}(\| y(t)\| _\infty , 1) = 1  \]

with the 1E$-$6 result, the ODE subroutine should attempt to maintain an accuracy of 1E$-$9 relative to 1. Therefore, the 1E$-$6 result should have almost three correct decimal places. At $t=2$, the first component of z is 6.58597048842310E$-$06, while its more accurate value is 6.58580950203220E$-$06, showing an $O(1$E$-10)$ error.

Troubleshooting

The ODE subroutine can fail for problems with unusual qualitative properties, such as finite escape time in the interval of integration (that is, the solution goes towards infinity at some finite time). In such cases, try testing with different subintervals and different levels of accuracy to gain some qualitative information about the behavior of the solution of the differential equation.