|
|
|
/****************************************************************/
/* 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;