NUMAL Section 5.2.1.1.1.1.G

BEGIN SECTION : 5.2.1.1.1.1.G (February, 1979)

PROCEDURE : DIFFSYS.

AUTHORS   : R.BULIRSCH AND J.STOER.

CONTRIBUTOR: K.DEKKER.

INSTITUTE: MATHEMATICAL CENTRE.

RECEIVED: 731231.

BRIEF DESCRIPTION:

    DIFFSYS SOLVES AN INITIAL VALUE PROBLEM ,FOR A SYSTEM OF FIRST
    ORDER ORDINARY DIFFERENTIAL EQUATIONS DY / DX = F(X,Y).
    THE METHOD IS RECOMMENDED IF HIGH ACCURACY IS DESIRED.
    DIFFSYS IS NOT SUITED FOR STIFF EQUATIONS.

KEYWORDS:

    INITIAL VALUE PROBLEMS,
    SYSTEM OF FIRST ORDER ORDINARY DIFFERENTIAL EQUATIONS.

CALLING SEQUENCE:

    THE HEADING OF THE PROCEDURE DIFFSYS READS:
    "PROCEDURE" DIFFSYS(X,XE,N,Y,DERIVATIVE,AETA,RETA,S,H0,OUTPUT);
    "VALUE" N;
    "INTEGER" N;
    "REAL" X,XE,AETA,RETA,H0;
    "ARRAY" Y,S;
    "PROCEDURE" DERIVATIVE,OUTPUT;
    "CODE" 33180;

    THE MEANING OF THE FORMAL PARAMETERS IS:
    X:      <VARIABLE>;
            THE INDEPENDENT VARIABLE  ;
            ENTRY: THE INITIAL VALUE X0;
            EXIT : THE FINAL VALUE XE;
    XE:     <ARITHMETIC EXPRESSION>;
            THE FINAL VALUE OF X (XE>=X);
    N:      <ARITHMETIC EXPRESSION>;
            THE NUMBER OF EQUATIONS;
    Y:      <ARRAY IDENTIFIER>;
            "REAL" "ARRAY" Y[1:N];
            THE DEPENDENT VARIABLE;
            ENTRY: THE  INITIAL  VALUES  OF  THE SYSTEM OF DIFFERENTIAL
                   EQUATIONS AT X=X0;
            EXIT : THE FINAL VALUES OF THE SOLUTION AT X=XE;
    DERIVATIVE: <PROCEDURE IDENTIFIER>;
            THE HEADING OF THIS PROCEDURE READS:
            "PROCEDURE" DERIVATIVE(X,Y,DY); "REAL" X; "ARRAY" Y,DY;
            THIS  PROCEDURE  SHOULD  DELIVER  THE  RIGHT  HAND  SIDE OF
            THE I-TH DIFFERENTIAL EQUATION AT THE POINT (X,Y) AS DY[I],
            I=1,...,N;
    AETA:   <ARITHMETIC EXPRESSION>;
            REQUIRED  ABSOLUTE  PRECISION  IN THE  INTEGRATION PROCESS;
    RETA:   <ARITHMETIC EXPRESSION>;
            REQUIRED  RELATIVE  PRECISION  IN THE  INTEGRATION PROCESS;
    S:      <ARRAY IDENTIFIER>;
            "REAL" "ARRAY" S[1:N];
            THE ARRAY S IS USED TO CONTROL THE ACCURACY OF THE COMPUTED
            VALUES OF Y;
            ENTRY: IT IS ADVISABLE TO SET S[I]=0, I=1,...,N;
            EXIT : THE MAXIMUM VALUE  OF ABS(Y[I]), ENCOUNTERED  DURING
                   INTEGRATION, IF THIS VALUE EXCEEDS THE VALUE OF S[I]
                   ON ENTRY;
    H0:     <VARIABLE>;
            THE INITIAL STEP TO BE TAKEN;
    OUTPUT: <PROCEDURE IDENTIFIER>;
            THE HEADING OF THIS PROCEDURE READS;
            "PROCEDURE" OUTPUT;
            THIS  PROCEDURE  IS  CALLED  AT THE END OF EACH INTEGRATION
            STEP ; THE USER CAN ASK FOR OUTPUT OF SOME PARAMETERS , FOR
            EXAMPLE X, Y, S.

PROCEDURES USED: NONE.

REQUIRED CENTRAL MEMORY : CIRCA 28* N  MEMORY PLACES.

METHOD AND PERFORMANCE:

    THE PROCEDURE  DIFFSYS  IS  A SLIGHT  MODIFICATION OF THE ALGORITHM
    PUBLISHED BY BULIRSCH AND STOER (SEE REF[1]) . BY THIS MODIFICATION
    INTEGRATION  FROM X0  UNTIL XE  CAN BE  PERFORMED  BY  ONE  CALL OF
    DIFFSYS. A NUMBER OF INTEGRATION STEPS ARE TAKEN, STARTING WITH THE
    INITIAL STEP H0. IN EACH INTEGRATION STEP A NUMBER OF SOLUTIONS ARE
    COMPUTED BY MEANS OF THE MODIFIED MIDPOINT RULE . EXTRAPOLATION  IS
    USED TO IMPROVE THESE SOLUTIONS,UNTIL THE REQUIRED ACCURACY IS MET.
    AN INTEGRATION STEP  IS REJECTED, IF THE ACCURACY REQUIREMENTS  ARE
    NOT FULFILLED  AFTER NINE EXTRAPOLATION STEPS . IN THESE CASES  THE
    INTEGRATION STEP IS REJECTED , AND INTEGRATION IS TRIED  AGAIN WITH
    THE INTEGRATION STEP HALVED.
    THE ALGORITHM IS FOR EACH STEP A VARIABLE ORDER METHOD (THE HIGHEST
    ORDER IS 14 ), AND USES A VARIABLE NUMBER OF  FUNCTION EVALUATIONS,
    DEPENDING ON THE ORDER (MINIMUM IS 3, MAXIMUM IS 217).
    THE ALGORITHM IS LESS SENSITIVE TO TOO SMALL VALUES  OF THE INITIAL
    STEPSIZE THAN THE ORIGINAL ALGORITHM ( SEE REF [2] );  HOWEVER  BAD
    GUESSES REQUIRE STILL SOME MORE COMPUTATIONS.

REFERENCES:

    [1]. R.BULIRSCH AND J.STOER.
         NUMERICAL  TREATMENT  OF  ORDINARY  DIFFERENTIAL  EQUATIONS BY
         EXTRAPOLATION METHODS.
         NUMERISCHE MATHEMATIK, VOLUME 8, PAGE 1-13, 1965.

    [2]. PHYLLIS FOX.
         A  COMPARATIVE  STUDY  OF  COMPUTER  PROGRAMS  FOR INTEGRATING
         DIFFERENTIAL EQUATIONS.
         COMMUNICATIONS OF THE A.C.M., VOLUME 15, PAGE 941-948, 1972.

    [3]. T.E.HULL, W.H.ENRIGHT, B.M.FELLEN AND A.E.SEDGWICK.
         COMPARING   NUMERICAL   METHODS  FOR   ORDINARY   DIFFERENTIAL
         EQUATIONS.
         SIAM JOURNAL ON NUMERICAL ANALYSIS,VOLUME 9,PAGE 603-635,1972.

EXAMPLE OF USE:

    THE FOLLOWING  PROGRAM  ILLUSTRATES  THE COSTS  AND  THE ACCURACIES
    WHICH ARE OBTAINED WHEN SOLVING A SYSTEM OF DIFFERENTIAL  EQUATIONS
    ARISING FROM THE RESTRICTED PROBLEM OF THREE BODIES ( SEE REF[1] ).
    THE SOLUTION IS A CLOSED ORBIT WITH PERIOD T=6.192169331396.

    "BEGIN"
        "INTEGER" PASSES,K;
        "REAL" X,XE,TIME,TOL,H0;
        "REAL" "ARRAY" Y,S[1:4];

        "PROCEDURE" DER(X,Y,DY); "REAL" X; "ARRAY" Y,DY;
        "BEGIN" "REAL" MU,MU1,Y1,Y2,Y3,Y4,S1,S2;
            MU:=1/82.45; MU1:=1-MU;
            PASSES:=PASSES+1;
            Y1:=Y[1]; Y2:=DY[1]:=Y[2]; Y3:=Y[3]; Y4:=DY[3]:=Y[4];
            S1:=(Y1+MU)**2+Y3**2; S2:=(Y1-MU1)**2+Y3**2;
            S1:=S1*SQRT(S1); S2:=S2*SQRT(S2);
            DY[2]:=Y1+2*Y4-MU1*(Y1+MU)/S1-MU*(Y1-MU1)/S2;
            DY[4]:=Y3-2*Y2-MU1*Y3/S1-MU*Y3/S2
        "END";

        "PROCEDURE" OUT;
        "BEGIN" K:=K+1;
            "IF" X>=XE "THEN"
            OUTPUT(61,"("2(-5ZD),2(4B+Z.3DB3DB3DB3D),-5ZD.3D,/")",K,
            PASSES,Y[1],Y[3],CLOCK-TIME)
        "END";

        OUTPUT(61,"(""("    THIS LINE AND THE FOLLOWING TEXT IS ")"
        "("PRINTED BY THIS PROGRAM")",//,
        "("    THE RESULTS WITH DIFFSYS - H0=.2 - ARE: ")",/,
        "("     K   DER.EV.         Y[1]                 Y[3]      ")",
        "("       TIME")",/")");
        "FOR" TOL:="-4,"-6,"-8,"-10,"-12 "DO"
        "BEGIN" PASSES:=K:=0; X:=0; XE:=6.192169331396;
            Y[1]:=1.2; Y[2]:=Y[3]:=0; Y[4]:=-1.04935750983;
            S[1]:=S[2]:=S[3]:=S[4]:=0; H0:=.2; TIME:=CLOCK;
            DIFFSYS(X,XE,4,Y,DER,TOL,TOL,S,H0,OUT);
        "END"
    "END"

    THIS LINE AND THE FOLLOWING TEXT IS PRINTED BY THIS PROGRAM:

    THE RESULTS WITH DIFFSYS - H0=.2 - ARE:
     K   DER.EV.         Y[1]                  Y[3]             TIME
     30   2591    +1.320 357 347 741     -.032 645 454 836      5.686
     33   3414    +1.200 078 037 878     -.000 053 906 067      7.455
     37   4213    +1.200 003 282 801     -.000 002 363 741      9.267
     44   4618    +1.199 999 999 711     -.000 000 000 095     10.242
     56   6299    +1.200 000 000 003     -.000 000 000 090     13.827

SOURCE TEXT:
"CODE" 33180;