NUMAL Section 5.2.1.1.1.2.F

BEGIN SECTION : 5.2.1.1.1.2.F (October, 1975)

AUTHOR: B.LINDBERG.

CONTRIBUTOR: K.DEKKER.

INSTITUTE: MATHEMATICAL CENTRE.

RECEIVED: 741101.

BRIEF DESCRIPTION:

    IMPEX SOLVES AN INITIAL VALUE PROBLEM,GIVEN AS AN AUTONOMOUS SYSTEM
    OF FIRST ORDER  DIFFERENTIAL  EQUATIONS , BY MEANS OF  THE IMPLICIT
    MIDPOINT RULE WITH SMOOTHING AND EXTRAPOLATION.
    AUTOMATIC STEPSIZE CONTROL IS PROVIDED.
    IN PARTICULAR THIS METHOD IS SUITABLE FOR THE INTEGRATION OF  STIFF
    DIFFERENTIAL EQUATIONS.

KEYWORDS:

    DIFFERENTIAL EQUATIONS,
    INITIAL VALUE PROBLEMS,
    STIFF EQUATIONS,
    IMPLICIT MIDPOINT RULE,
    SMOOTHING,
    EXTRAPOLATION.

CALLING SEQUENCE:

    THE HEADING OF THE PROCEDURE IMPEX READS:
    "PROCEDURE" IMPEX (N, T0, TEND, Y0, DERIV, AVAILABLE, H0, HMAX,
                      PRESCH, EPS, WEIGHTS, UPDATE, FAIL, CONTROL);
    "VALUE" N;
    "INTEGER" N;
    "REAL" T0,TEND,H0,HMAX,EPS;
    "BOOLEAN" PRESCH,FAIL;
    "ARRAY" Y0,WEIGHTS;
    "BOOLEAN" "PROCEDURE" AVAILABLE;
    "PROCEDURE" DERIV,UPDATE,CONTROL;
    "CODE" 33135;

    IMPEX:   INTEGRATES THE SYSTEM OF DIFFERENTIAL EQUATIONS , GIVEN AS
             DY/DT=F(T,Y), FROM T0 UNTIL TEND.

    THE MEANING OF THE FORMAL PARAMETERS IS:
    N:      <ARITHMETIC EXPRESSION>;
            THE NUMBER OF EQUATIONS;
    T0:     <ARITHMETIC EXPRESSION>;
            THE INITIAL VALUE OF THE INDEPENDENT VARIABLE;

    TEND:   <ARITHMETIC EXPRESSION>;
            THE FINAL VALUE OF THE INDEPENDENT VARIABLE;
    Y0:     <ARRAY IDENTIFIER>;
            "REAL" "ARRAY" Y0[1:N];
            ENTRY: THE  INITIAL  VALUES  OF  THE SYSTEM OF DIFFERENTIAL
                   EQUATIONS: Y0[I] AT T=T0;
    DERIV:  <PROCEDURE" IDENTIFIER>;
            THE HEADING OF THIS PROCEDURE READS:
            "PROCEDURE" DERIV(T,Y,F,N);
            "INTEGER" N; "REAL" T; "ARRAY" Y,F;
            THIS PROCEDURE  SHOULD DELIVER  THE VALUE OF F(T,Y)  IN THE
            ARRAY F[1:N];
    AVAILABLE: <PROCEDURE IDENTIFIER>;
            THE HEADING OF THIS PROCEDURE READS:
            "BOOLEAN" "PROCEDURE" AVAILABLE(T,Y,A,N);
            "INTEGER" N; "REAL" T; "ARRAY" Y,A;
            IF AN  ANALYTIC  EXPRESSION OF  THE JACOBIAN MATRIX  AT THE
            POINT (T,Y) IS NOT AVAILABLE, THIS PROCEDURE SHOULD DELIVER
            THE VALUE "FALSE";
            OTHERWISE THE  PROCEDURE  SHOULD  DELIVER THE VALUE "TRUE",
            AND THE  JACOBIAN  MATRIX  SHOULD BE ASSIGNED TO THE  ARRAY
            A[1:N,1:N];
    H0:     <ARITHMETIC EXPRESSION>;
            THE INITIAL STEPSIZE;
    HMAX:   <ARITHMETIC EXPRESSION>;
            MAXIMAL  STEPSIZE  BY WHICH  THE INTEGRATION  IS PERFORMED;
    PRESCH: <BOOLEAN EXPRESSION>;
            INDICATOR FOR CHOICE OF STEPSIZE;
            THE STEPSIZE IS AUTOMATICALLY CONTROLLED IF PRESCH="FALSE";
            OTHERWISE  THE STEPSIZE  HAS TO BE PRESCRIBED , EITHER ONLY
            INITIALLY OR ALSO BY THE PROCEDURE CONTROL;
    EPS:    <ARITHMETIC EXPRESSION>;
            BOUND FOR THE ESTIMATE OF THE LOCAL ERROR;
    WEIGHTS: <ARRAY IDENTIFIER>;
            "REAL" "ARRAY" WEIGHTS[1:N];
            WEIGHTS FOR THE COMPUTATION  OF THE WEIGHTED EUCLIDEAN NORM
            OF THE ERRORS;
            ENTRY: INITIAL WEIGHTS;
            NOTE THAT THE CHOICE  WEIGHTS[I] = 1 IMPLIES  AN ESTIMATION
            OF THE ABSOLUTE ERROR , WHEREAS WEIGTHS[I] = Y[I] DEFINES A
            RELATIVE ERROR;
    UPDATE: <PROCEDURE IDENTIFIER>;
            THE HEADING OF THIS PROCEDURE READS:
            "PROCEDURE" UPDATE(WEIGHTS,Y2,N);
            "INTEGER" N; "ARRAY" WEIGHTS,Y2;
            THIS PROCEDURE  MAY CHANGE THE ARRAY WEIGHTS , ACCORDING TO
            THE VALUE OF AN APPROXIMATION FOR Y(T) , GIVEN IN THE ARRAY
            Y2[1:N];

    FAIL:   <BOOLEAN EXPRESSION>;
            EXIT :
            IF  THE PROCEDURE  FAILS TO SOLVE  THE SYSTEM OF EQUATIONS,
            FAIL WILL HAVE THE VALUE "TRUE" UPON EXIT;
            THIS MAY OCCUR  UPON DIVERGENCE  OF THE ITERATION  PROCESS,
            USED IN THE MIDPOINT RULE , WHILE INTEGRATION  IS PERFORMED
            WITH A USER DEFINED PRESCRIBED STEPSIZE;
    CONTROL: <PROCEDURE IDENTIFIER>;
            THE HEADING OF THIS PROCEDURE READS:
            "PROCEDURE" CONTROL(TPRINT,T,H,HNEW,Y,ERROR,N);
            "INTEGER" N; "REAL" TPRINT,T,H,HNEW; "ARRAY" Y,ERROR;
            CONTROL IS CALLED ON ENTRY OF IMPEX, AND FURTHER AS SOON AS
            THE INEQUALITY  TPRINT <= T  HOLDS;
            DURING A CALL OF CONTROL  PRINTING OF RESULTS AND
            CHANGE OF STEPSIZE (IF PRESCH = "TRUE")  IS THEN  POSSIBLE;
            THE MEANING OF THE FORMAL PARAMETERS IS:
            TPRINT: <VARIABLE>;
                    ENTRY: THE  VALUE  OF THE  INDEPENDENT  VARIABLE AT
                           WHICH A CALL OF CONTROL WAS DESIRED;
                    EXIT:  A NEW VALUE  (TPRINT>T)  AT WHICH  A CALL OF
                           CONTROL IS DESIRED;
            T:      <VARIABLE>;
                    THE ACTUAL VALUE OF THE INDEPENDENT VARIABLE, UP TO
                    WHICH INTEGRATION HAS BEEN PERFORMED;
            H:      <VARIABLE>;
                    HALVE THE ACTUAL STEPSIZE;
            HNEW:   <VARIABLE>;
                    THE NEW STEPSIZE;
                    IF PRESCH="TRUE", THEN THE USER MAY PRESCRIBE A NEW
                    STEPSIZE BY CHANGING HNEW;
            Y:      <ARRAY IDENTIFIER>;
                    "REAL" "ARRAY" Y[1:5,1:N];
                    THE VALUE  OF THE DEPENDENT VARIABLE  AND ITS FIRST
                    FOUR DIVIDED DIFFERENCES  AT THE POINT T  ARE GIVEN
                    IN THIS ARRAY;
           ERROR:   <ARRAY IDENTIFIER>;
                    "REAL" "ARRAY" ERROR[1:3];
                    THE ELEMENTS  OF THIS ARRAY  CONTAIN  THE FOLLOWING
                    ERRORS:
                    ERROR[1]: THE LOCAL ERROR;
                    ERROR[2]: THE GLOBAL ERROR OF SECOND ORDER IN H;
                    ERROR[3]: THE GLOBAL ERROR OF FOURTH ORDER IN H;
           N:       <VARIABLE>;
                    THE NUMBER OF EQUATIONS;
           EXAMPLE OF USE: SEE EXAMPLE OF USE OF THE PROCEDURE IMPEX;

DATA AND RESULTS:
    FOR DATA, SEE REF[1].
    THE RESULTS OF THE INTEGRATION ARE ATTAINABLE THROUGH THE PROCEDURE
    CONTROL , WHICH IS CALLED AT SPECIFIED, USER DEFINED, VALUES OF THE
    INDEPENDENT  VARIABLE . IN PARTICULAR , THE VALUES OF THE DEPENDENT
    VARIABLE  AT THE ENDPOINT OF INTEGRATION  ARE OBTAINED BY A CALL OF
    CONTROL WITH TPRINT=TEND.

PROCEDURES USED:

    INIVEC   = CP31010,
    INIMAT   = CP31011,
    MULVEC   = CP31020,
    MULROW   = CP31021,
    DUPVEC   = CP31030,
    DUPROWVEC= CP31032,
    DUPMAT   = CP31035,
    VECVEC   = CP34010,
    MATVEC   = CP34011,
    MATMAT   = CP34013,
    ELMVEC   = CP34020,
    ELMROW   = CP34024,
    DEC      = CP34300,
    SOL      = CP34051.

REQUIRED CENTRAL MEMORY:

    EXECUTION FIELD LENGTH: CIRCA  N * ( 23 + 2 * N ) (DECIMAL).

RUNNING TIME: DEPENDS  STRONGLY  ON THE DIFFERENTIAL EQUATION TO SOLVE.

LANGUAGE: ALGOL 60.

METHOD AND PERFORMANCE:

    THE INTEGRATION METHOD (REF[1]) IS BASED ON  THE COMPUTATION OF TWO
    INDEPENDENT SOLUTIONS Y(T,H) AND Y(T,H/2)  BY THE IMPLICIT MIDPOINT
    RULE. PASSIVE SMOOTHING  AND PASSIVE EXTRAPOLATION  IS PERFORMED TO
    OBTAIN STABILITY  AND HIGH ACCURACY . THE ALGORITHM  USES  FOR EACH
    STEP AT LEAST THREE FUNCTION EVALUATIONS, AND ON CHANGE OF STEPSIZE
    OR AT SLOW CONVERGENCE IN THE ITERATION PROCESS AN APPROXIMATION OF
    THE JACOBIAN MATRIX (COMPUTED BY DIVIDED DIFFERENCES  OR EXPLICITLY
    SPECIFIED  BY THE USER) . IF THE COMPUTED LOCAL ERROR  EXCEEDS  THE
    TOLERANCE , THE LAST STEP IS REJECTED. MOREOVER , TWO GLOBAL ERRORS
    ARE COMPUTED.

REFERENCES:

    [1]. B.LINDBERG.
         IMPEX 2 , A PROCEDURE  FOR THE SOLUTION  OF  SYSTEMS OF  STIFF
         DIFFERENTIAL EQUATIONS.
         ROYAL INSTITUTE OF TECHNOLOGY,STOCKHOLM. TRITA-NA-7303 (1973).

    [2]. (TO APPEAR)
         COLLOQUIUM STIFF DIFFERENTIAL EQUATIONS 3 (DUTCH).
         M.C. SYLLABUS 15.3 (1974), MATHEMATICAL CENTRE.

EXAMPLE OF USE:
    CONSIDER THE AUTONOMOUS SYSTEM OF DIFFERENTIAL EQUATIONS:
    DY[1]/DX = .2 * ( Y[2] - Y[1] ),
    DY[2]/DX = 10 * Y[1] - ( 60 - Y[3]/8 ) * Y[2] + Y[3]/8,
    DY[3]/DX = 1,
    WITH  INITIAL  CONDITIONS  AT  X=0 : Y[1]=Y[2]=Y[3]=0 (SEE REF[2]).
    THE SOLUTION  AT SEVERAL POINTS  IN THE INTERVAL  [0, 400]  MAY  BE
    OBTAINED BY THE FOLLOWING PROGRAM:
    (THE  SOLUTION  AT  X=400  IS: Y[1]=22.24222011,  Y[2]=27.11071335)

"BEGIN" "INTEGER" N,NFE,NJE,POINT;
    "REAL" T,TEND,EPS,HMAX,L,H2,TIME;
    "ARRAY" Y,SW[1:3],PRINT[1:5];
    "BOOLEAN" FAIL;

    "PROCEDURE" LIPEST(L,Y,EPS,T,F,N);
    "REAL" T,L,EPS; "ARRAY" Y; "INTEGER" N; "PROCEDURE" F;
    "BEGIN" "REAL" N1,N2; "INTEGER" I,IT; "ARRAY" F1,F2,Z,X[1:N];
        "REAL" "PROCEDURE" NORM(Y); "ARRAY" Y;
        NORM:=SQRT(VECVEC(1,N,0,Y,Y));
        DUPVEC(1,N,0,Z,Y);
        "FOR" I:=1 "STEP" 1 "UNTIL" N "DO"
        X[I]:="IF" Y[I]=0 "THEN" EPS "ELSE" (1+EPS)*Y[I];
        N1:=NORM(X)*EPS; F(T,X,F1,N);
        "FOR" IT:=1 "STEP" 1 "UNTIL" 5 "DO"
        "BEGIN" F(T,Z,F2,N);
            ELMVEC(1,N,0,F2,F1,-1);
            N2:=N1/NORM(F2);
            DUPVEC(1,N,0,Z,X); ELMVEC(1,N,0,Z,F2,N2)
        "END";
        F(T,Z,F2,N);
        ELMVEC(1,N,0,F2,F1,-1);
        L:=NORM(F2)/N1
    "END" LIPEST;

    "PROCEDURE" F(T,Y,F1,N); "VALUE" T; "REAL" T; "ARRAY" Y,F1;
    "INTEGER" N;
    "BEGIN" NFE:=NFE+1;
        F1[1]:=0.2*(Y[2]-Y[1]);
        F1[2]:=10*Y[1]-(60-.125*Y[3])*Y[2]+.125*Y[3];
        F1[3]:=1
    "END";

    "BOOLEAN" "PROCEDURE" AVAILABLE(T,Y,A,N);
    "INTEGER" N; "REAL" T; "ARRAY" Y,A;
    "BEGIN" NJE:=NJE+1; AVAILABLE:="TRUE";
        A[1,1]:=-.2; A[1,2]:=.2; A[1,3]:=A[3,1]:=A[3,2]:=A[3,3]:=0;
        A[2,1]:=10; A[2,2]:=.125*Y[3]-60; A[2,3]:=.125*(1+Y[2])
    "END"

    "PROCEDURE" UPDATE(SW,R1,N); "INTEGER" N; "ARRAY" SW,R1;
    "BEGIN" "REAL" S1,S2; "INTEGER" I;
        "FOR" I:=1 "STEP" 1 "UNTIL" N "DO"
        "BEGIN" S1:=1/SW[I]; S2:=ABS(R1[I]);
            "IF" S1<S2 "THEN" SW[I]:=1/S2
        "END"
    "END";

    "PROCEDURE" CONTROL(TP,T,H,HNEW,Y,ERR,N);
    "REAL" TP,T,H,HNEW; "ARRAY" Y,ERR; "INTEGER" N;
    "BEGIN" "INTEGER" I;
        "ARRAY" C[3:5],X[1:N];
        "REAL" S,S2,S3,S4,C1;
  NEXT: S:=(T-TP)/H;
        S2:=S*S; S3:=S2*S; S4:=S3*S;
        C[3]:=(S2-S)/2;
        C[4]:=-S3/6+S2/2-S/3;
        C[5]:=S4/24-S3/4+11*S2/24-S/4;
        "FOR" I:=1 "STEP" 1 "UNTIL" N "DO"
        X[I]:=Y[1,I]-S*Y[2,I]+C[3]*Y[3,I]+C[4]*Y[4,I]+C[5]*Y[5,I];
        OUTPUT(61,"("3ZD.2D2B,D.D"-D2B,2(+D.8D"-D2B),2(4ZD),3ZD.2D,
               /")",TP,ERR[3],X[1],X[2],NFE,NJE,CLOCK-TIME);
        "IF" TP<TEND "THEN"
        "BEGIN" POINT:=POINT+1; TP:=PRINT[POINT];
                "IF" TP<=T "THEN" "GOTO" NEXT
        "END"
    "END" CONTROL;

    N:=3; NJE:=NFE:=0; T:=0; TEND:=400; EPS:="-5; HMAX:=400;
    Y[1]:=Y[2]:=Y[3]:=0; SW[1]:=SW[2]:=SW[3]:=1;
    PRINT[1]:=.1; PRINT[2]:=1; PRINT[3]:=10; PRINT[4]:=100;
    PRINT[5]:=400; POINT:= 0;
    LIPEST(L,Y,"-5,T,F,N);
    H2:=(EPS*320)**(1/5)/(4*L);
    OUTPUT(61,"(""("EPS=")",D.2D"-D,/,"("INTERVAL OF INTEGRATION=(")",
    3ZD,"(",")",3ZD,"(")")",/,"("MAXIMALLY ALLOWED STEPSIZE=")",
    D.2D"-D,//")",EPS,T,TEND,HMAX);
    OUTPUT(61,"(""("LIPSCHCONST=")",BD.3D"+D,/,"("STARTING STEPSIZE")"
    "("=")",BD.2D"+D,/,"("FUNCTIONAL EVAL=")",4ZD,//")",L,H2,NFE);
    TIME:=CLOCK;
    OUTPUT(61,"(""("    X     ERROR       Y[1]            Y[2]")"
    "("         NFE  NJE   TIME")",/")");
    IMPEX(N,T,TEND,Y,F,AVAILABLE,H2,HMAX,"FALSE",EPS,SW,UPDATE,FAIL,
    CONTROL);
    OUTPUT(61,"("/"("NO OF FUNCTIONAL EVALUATIONS= ")",3ZD,/,
    "("NO OF JACOBEAN EVALUATIONS= ")",3ZD,/")",NFE,NJE)
"END"

THIS PROGRAM DELIVERS:

EPS=1.00"-5
INTERVAL OF INTEGRATION=(   0, 400)
MAXIMALLY ALLOWED STEPSIZE=4.00" 2

LIPSCHCONST= 6.003"+1
STARTING STEPSIZE= 1.32"-3
FUNCTIONAL EVAL=    7

    X     ERROR       Y[1]            Y[2]         NFE  NJE   TIME
   0.00  0.0" 0  +0.00000000" 0  +0.00000000" 0      7    0   0.01
   0.10  6.3"-7  +1.49614151"-6  +1.74013792"-4     46    4   0.72
   1.00  1.5"-6  +1.91041887"-4  +2.08361269"-3     85    8   1.48
  10.00  8.7"-7  +1.30147663"-2  +2.34487800"-2    119    9   1.99
 100.00  1.3"-5  +3.06302487"-1  +3.27552180"-1    225   13   3.47
 400.00  1.4"-5  +2.22406546" 1  +2.71090507" 1    556   30   7.51

NO OF FUNCTIONAL EVALUATIONS=  556
NO OF JACOBEAN EVALUATIONS=   30

SOURCE TEXT(S):
"CODE" 33135;