Resources

ODE: A Differential Equation Solver

 /****************************************************************/
 /*          S A S   S A M P L E   L I B R A R Y                 */
 /*                                                              */
 /*    NAME: ODE                                                 */
 /*   TITLE: ODE: A Differential Equation Solver                 */
 /* PRODUCT: IML                                                 */
 /*  SYSTEM: ALL                                                 */
 /*    KEYS:                                                     */
 /*   PROCS: IML                                                 */
 /* SUPPORT: GHG                         UPDATE:                 */
 /*     REF:                                                     */
 /*    MISC:                                                     */
 /*                                                              */
 /****************************************************************/


 PROC IML  ;

 /*-------------------------------------------------------------*/
 /* EXAMPLE1: A SIMPLE PROGRAM TO SOLVE THE DIFFERENTIAL        */
 /*           EQUATION                                          */
 /*                                                             */
 /*           dx                                                */
 /*          --- = - t x     with x(0) = .5                     */
 /*           dt                                                */
 /*                                                             */
 /* THE USER MUST SUPPLY                                        */
 /*                                                             */
 /*          1) THE (ORDERED) LIST OF POINTS AT WHICH THE       */
 /*             SOLUTION IS TO BE  REPORTED.                    */
 /*             THE INTEGRATION IS TO BE CARRIED UP TO .3       */
 /*             AND ALSO REPORT THE RESULTS AT .1 and .2        */
 /*                                                             */
 /*          2) THE RIGHT HAND SIDE   (here it is : -t*x)       */
 /*                                                             */
 /*          3) THE STEP-SIZE VECTOR                            */
 /*             THE INTEGRATION IS TO BE CARRIED  WITH          */
 /*                 A MAXIMUM STEP SIZE OF 1                    */
 /*                 A MINIMUM STEP SIZE OF 1E-12                */
 /*                 AN INITIAL STEP SIZE OF 1E-5                */
 /*          THESE NUMBERS ARE SPECIFIED IN H (NOT ORDERED)     */
 /*                                                             */
 /* THE USER CAN SUPPLY                                         */
 /*                                                             */
 /*          1) THE JACOBIAN OF THE RIGHT HAND SIDE  (HERE : -T)*/
 /*          IF NOT SUPPLIED FINITE DIFFERENCE METHODS WILL     */
 /*          BE  USED TO ESTIMATE THE JACOBIAN                  */
 /*                                                             */
 /* RESULTS ARE                                                 */
 /*                                                             */
 /*          Z[1,1] =  4.975062395981900000E-01                 */
 /*          Z[1,2] =  4.900993366551400000E-01                 */
 /*          Z[1,3] =  4.779987409188300000E-01                 */
 /*-------------------------------------------------------------*/

 START FUN(T,X) GLOBAL(COUNT);       /* DEFINE THE RHS          */
     COUNT = COUNT+1;                /* WE COUNT THE NUMBER OF  */
                                     /* FUNCTION EVALUATIONS    */
     V = -T*X;
     RETURN(V);
 FINISH;

 START JAC(T,X) ;                    /* DEFINE THE JACOBIAN     */
     V = -T;
     RETURN(V);
 FINISH;


 T = { 0 .1 .2 .3};                  /* SPECIFY THE ORDERED LIST*/
                                     /* OF POINTS THAT DESCRIBE */
                                     /* THE INTERVALS           */
 C = .5;                             /* SPECIFY THE INITIAL     */
                                     /* VALUE                   */

 H = {1.E-12 1 1.E-5};               /* SPECIFY THE H VECTOR    */


 COUNT = 0;                          /* SET THE COUNTER         */

 CALL ODE(Z,"FUN",C,T,H) EPS = 1.E-12 J="JAC";
                                     /* CALL ODE                */

 PRINT COUNT;
 PRINT Z[FORMAT =E25.15];

 /*-------------------------------------------------------------*/
 /* EXAMPLE2: A SIMPLE PROGRAM TO SOLVE THE DIFFERENTIAL        */
 /*           EQUATION                                          */
 /*                                                             */
 /*           dx                        | 1 |                   */
 /*          --- =  A x     with x(0) = |   |                   */
 /*           dt                        | 2 |                   */
 /*                                                             */
 /*            |           |                                    */
 /*            | -1      1 |                                    */
 /*        A = |           |                                    */
 /*            |  1      2 |                                    */
 /*            |           |                                    */
 /*                                                             */
 /* THE USER MUST SUPPLY                                        */
 /*                                                             */
 /*          1) THE (ORDERED) LIST OF POINTS AT WHICH THE       */
 /*             SOLUTION IS TO BE  REPORTED.                    */
 /*             THE INTEGRATION IS TO BE CARRIED UP TO .3       */
 /*             AND ALSO REPORT THE RESULTS AT .1 AND .2        */
 /*                                                             */
 /*          2) THE RIGHT HAND SIDE   (HERE IT IS : -T*X)       */
 /*                                                             */
 /*          3) THE STEP-SIZE VECTOR                            */
 /*             THE INTEGRATION IS TO BE CARRIED  WITH          */
 /*                 A MAXIMUM STEP SIZE OF 1                    */
 /*                 A MINIMUM STEP SIZE OF 1E-12                */
 /*                 AN INITIAL STEP SIZE OF 1E-5                */
 /*          THESE NUMBERS ARE SPECIFIED IN H (NOT ORDERED)     */
 /*                                                             */
 /*          4) THE INITIAL VECTOR                              */
 /*                                                             */
 /* THE USER CAN SUPPLY                                         */
 /*                                                             */
 /*          1) THE JACOBIAN OF THE RIGHT HAND SIDE  (HERE : -T)*/
 /*          IF NOT SUPPLIED FINITE DIFFERENCE METHODS WILL     */
 /*          BE  USED TO ESTIMATE THE JACOBIAN                  */
 /* RESULTS                                                     */
 /*                                                             */
 /*          FOR T= 1 THE RESULT IS STORED IN Z[,1]             */
 /*                                                             */
 /*            Z[1,1] = 1.059276279257100000E+00                */
 /*            Z[2,1] = 7.554305542309500000E-01                */
 /*                                                             */
 /*          FOR T= 2 THE RESULT IS STORED IN Z[,2]             */
 /*                                                             */
 /*            Z[1,2] = 7.504499669077400000E-01                */
 /*            Z[2,2] = 4.711431015044000000E-01                */
 /*                                                             */
 /*          FOR T= 3 THE RESULT IS STORED IN Z[,3]             */
 /*                                                             */
 /*            Z[1,3] = 5.141930081253500000E-01                */
 /*            Z[2,3] = 3.183434053151600000E-01                */
 /*                                                             */
 /*-------------------------------------------------------------*/


 START FUN2(T,X) GLOBAL(A) ;           /* DEFINE THE RHS        */
                                       /* PASS A AS A PARAMETER */
     V = A*X;
     RETURN(V);
 FINISH;


 T = { 0  1  2  3 };                   /* DECLARE THE INTERVALS */

 C = { 1 , 2 };                        /* DECLARE THE INTIAL    */
                                       /* VECTOR C              */
 A = { -1  1 , 1  -2 } ;               /* DEFINE A              */

 H = {1.E-9 1 1.E-5};                  /* SET THE H VECTOR      */

 CALL ODE(Z,"FUN2",C,T,H) EPS = 1.E-5; /* CALL THE INTEGRATOR   */

 DO I = 1 TO 3;
   D = Z[,I];
   PRINT D[FORMAT =E25.15];
 END;