'**************************************************************** '* NUMERICAL SOLUTION OF A STIFF SYSTEM OF FIRST 0RDER ORDINARY * '* DIFFERENTIAL EQUATIONS Y':=F(X,Y) BY ROSENBROCK METHOD. * '* ------------------------------------------------------------ * '* SAMPLE RUN: * '* Example #1: * '* (Solve set of differential equations (N=2): * '* F(1) = Y(1) * Y(2) + COS(X) - HALF * SIN(TWO * X) * '* F(2) = Y(1) * Y(1) + Y(2) * Y(2) - (ONE + SIN(X)) * '* Find values of F(1), F(2) at X=1.5). * '* * '* SOLUTION AT X= 1.50000000000000E+0000 * '* Y(1) = 1.236006095804576 * '* Y(2) = -.1049268945803322 * '* * '* LAST STEP SIZE = 3.089293500117273D-04 * '* ERROR CODE = 1 * '* * '* Example #2: * '* (Solve set of differential equations (N=5): * '* F(1) = Y(2) * '* F(2) = Y(3) * '* F(3) = Y(4) * '* F(4) = Y(5) * '* F(5) = (45.0 * Y(3) * Y(4) * Y(5) - * '* 40.0 * Y(4) * Y(4) * Y(4)) / (9.0 * Y(3) * Y(3)) * '* Find values of F(1), F(2), ..., F(5) at X=1.5). * '* * '* SOLUTION AT X= 1.50000000000000E+0000 * '* Y(1) = 4.363967162542581 * '* Y(2) = 4.000019057753676 * '* Y(3) = 2.82847148934375 * '* Y(4) = 5.641335228805289D-05 * '* Y(5) = -3.77130489085171 * '* * '* LAST STEP SIZE = 7.626269049659622D-05 * '* ERROR CODE = 1 * '* ------------------------------------------------------------ * '* Ref.: From Numath Library By Tuan Dang Trong in Fortran 77 * '* [BIBLI 18]. * '* * '* Basic Release 1.0 By J-P Moreau, Paris * '* (www.jpmoreau.fr) * '**************************************************************** ' LIST OF USED SUBROUTINES (HERE INCLUDED): ' ======================================== ' 500 FCN DEFINE SYSTEM OF DIFFERENTIAL EQUATIONS (TWO EXAMPLES) ' 600 IMAX MAXIMUM OF TWO INTEGERS ' 610 IMIN MINIMUM OF TWO INTEGERS ' 620 XMAX MAXIMUM OF TWO REAL NUMBERS ' 630 XMIN MINIMUM OF TWO REAL NUMBERS ' 650 SIGN EMULATION OF FUNCTION SIGN OF FORTRAN ' 1000 ROS4 NUMERICAL SOLUTION OF A STIFF (OR DIFFERENTIAL ALGEBRAIC) ' SYSTEM OF FIRST 0RDER ORDINARY DIFFERENTIAL EQUATIONS ' MY'=F(X,Y) ' 2000 RO4COR CORE INTEGRATOR FOR ROS4 ' 3000 SHAMP DEFINE CONSTANTS A21... TO D4 (METHOD 1) ' 3001 GRK4A DEFINE CONSTANTS A21... TO D4 (METHOD 3) ' 3002 GRK4T DEFINE CONSTANTS A21... TO D4 (METHOD 2 USED HERE) ' 3003 VELDD DEFINE CONSTANTS A21... TO D4 (METHOD 4) ' 3004 VELDS DEFINE CONSTANTS A21... TO D4 (METHOD 5) ' 3005 LSTAB DEFINE CONSTANTS A21... TO D4 (METHOD 6) ' 4000 DECB MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION OF A BANDED ' MATRIX WITH LOWER BANDWIDTH MLDE AND UPPER BANDWIDTH MUE ' 4500 DECA GENERAL MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION ' 5000 SOLB SOLUTION OF A BANDED LINEAR SYSTEM, E*X = B ' (E IS A TRIANGULARIZED MATRIX OBTAINED FROM DECB). ' 5500 SOL SOLUTION OF A GENERAL LINEAR SYSTEM, E*X = B ' (E IS A TRIANGULARIZED MATRIX OBTAINED FROM DECA). ' -------------------------------------------------------------------------- 'NOTE: THE BANDED MATRIX BRANCH HAS NOT BEEN TESTED HERE, HOWEVER ' IS FULLY IMPLEMENTED. DefDbl A-H, O-Z DefInt I-N 'constants NMX = 30 'Maximum size of temporary vectors TMP1 to TMP11 and ITMP HALF = 0.5 ONE = 1# TEN = 10# TWO = 2# XNINE = 9# ZERO = 0# Dim B(NMX), F(NMX), Y(NMX), YY(NMX) Dim FJAC(5, 5), E(5, 5), FMAS(5, 5) 'variables for statistics (optional use) 'XNFCN,XNSTEP,XNJAC,XNACCPT,NREJCT,XNDEC,XNSOL '(Long integers are simulated by real numbers). 'begin main program 'Initialize parameters (see 1000 ROS4) N = 2 'DIMENSION OF THE SYSTEM (N=5 for example #2) IFCN = 1 'FCN(N,X,Y,F) MAY DEPEND ON X X = ZERO 'INITIAL X-VALUE XEND = 1.5 'FINAL X-VALUE (XEND-X MAY BE POSITIVE OR NEGATIVE) H = 0.001 'INITIAL STEP SIZE GUESS (0.01 FOR example #2) RTOL = 0.00000001 'RELATIVE XERROR TOLERANCE (HERE SCALAR) ATOL = 0.0000000001 'ABSOLUTE XERROR TOLERANCE (HERE SCALAR) ITOL = 0 'BOTH RTOL AND ATOL ARE SCALARS IJAC = 0 'JACOBIAN IS COMPUTED INTERNALLY BY FINITE 'DIFFERENCES, Subroutine "JAC" IS NEVER CALLED MLJAC = N 'JACOBIAN IS A FULL MATRIX. THE LINEAR ALGEBRA 'IS DONE BY FULL-MATRIX GAUSS-ELIMINATION IDFX = 0 'DF/DX IS COMPUTED INTERNALLY BY FINITE 'DIFFERENCES, Subroutine "DFX" IS NEVER CALLED IMAS = 0 'M IS SUPPOSED TO BE THE IDENTITY 'MATRIX, Subroutine "MAS" IS NEVER CALLED MLMAS = N 'MLMAS=N: THE FULL MATRIX CASE. THE LINEAR ALGEBRA 'IS DONE BY FULL-MATRIX GAUSS-ELIMINATION IOUT = 0 'Subroutine SOLOUT IS NEVER CALLED LE1 = N 'IF MLJAC=N (FULL JACOBIAN) LJAC = N 'IF MLJAC=N (FULL JACOBIAN) LMAS = 0 'IF IMAS=0 LIWORK = N + 2 'DECLARED LENGTH OF ARRAY "IWORK" LWORK = N * (LJAC + LMAS + LE1 + 8) + 5 'DECLARED LENGTH OF ARRAY "WORK" Dim WORK(LWORK) Dim IWORK(LIWORK) 'Temporary vectors Dim TMP1(NMX), TMP2(NMX), TMP3(NMX), TMP4(NMX), TMP5(NMX), TMP6(NMX) Dim TMP7(NMX), TMP8(NMX), TMP9(NMX), TMP10(NMX), TMP11(NMX) Dim ITMP(NMX) For I = 1 To LWORK WORK(I) = ZERO 'This triggers default values (see 1000 ROS4) Next I For I = 1 To LIWORK IWORK(I) = 0 Next I Y(1) = HALF 'INITIAL VALUES FOR Y Y(2) = HALF 'In example #2, Y(1) = Y(2) = ... = Y(5) = ONE Cls Print Print " Computing..." 'call Rosenbrock SUBROUTINE with appropriate parameters GoSub 1000 'call ROS4(N,IFCN,X,Y,XEND,H, 'RTOL,ATOL,ITOL, 'IJAC,MLJAC,MUJAC,IDFX, 'IMAS,MLMAS,MUMAS, 'IOUT,WORK,LWORK,IWORK,LIWORK,IDID) 'print results Cls Print Print " SOLUTION AT X="; X For I = 1 To N Print " Y("; I; ") = "; Y(I) Next I Print Print " LAST STEP SIZE ="; H Print " ERROR CODE = "; IDID Print INPUT "", RR$ End 'of main program 'define example #1 500 'FCN(N,XX,YY,F) F(1) = YY(1) * YY(2) + Cos(XX) - HALF * Sin(TWO * XX) F(2) = YY(1) * YY(1) + YY(2) * YY(2) - (ONE + Sin(XX)) Return 'define example #2 '500 'FCN(N,XX,YY,F) ' F(1) = YY(2); ' F(2) = YY(3); ' F(3) = YY(4); ' F(4) = YY(5); ' F(5) = (45.0 * YY(3) * YY(4) * YY(5) - ' 40.0 * YY(4) * YY(4) * YY(4)) / (NINE * YY(3) * YY(3)) 'return 600 'IMAX(ia,ib) If ia > ib Then IMAX = ia Else IMAX = ib End If Return 610 'IMIN(ia,ib) If ia < ib Then IMIN = ia Else IMIN = ib End If Return 620 'XMAX(xa,xb) If xa > xb Then XMAX = xa Else XMAX = xb End If Return 630 'XMIN(xa,xb) If xa < xb Then XMIN = xa Else XMIN = xb End If Return 650 'SIGN(xa,xb) If xb < 0 Then SIGN = -Abs(xa) Else SIGN = Abs(xa) End If Return '********************************************************************** 1000 'ROS4 (N,IFCN,X,Y,XEND,H,RTOL,ATOL,ITOL,IJAC,MLJAC,MUJAC,IDFX, ' IMAS,MLMAS,MUMAS,IOUT,WORK,LWORK,IWORK,LIWORK,IDID) ' --------------------------------------------------------------------- ' NUMERICAL SOLUTION OF A STIFF (OR DIFFERENTIAL ALGEBRAIC) ' SYSTEM OF FIRST 0RDER ORDINARY DIFFERENTIAL EQUATIONS MY'=F(X,Y). ' THIS IS AN EMBEDDED ROSENBROCK METHOD OF ORDER (3)4 ' (WITH STEP SIZE CONTROL). ' C.F. SECTION IV.7 ' ' AUTHORS: E. HAIRER AND G. WANNER ' UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES ' CH-1211 GENEVE 24, SWITZERLAND ' E-MAIL: HAIRER@CGEUGE51.BITNET, WANNER@CGEUGE51.BITNET ' ' THIS CODE IS PART OF THE BOOK: ' E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL ' EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS. ' SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS, ' SPRINGER-VERLAG (1990) ' ' VERSION OF OCTOBER 12, 1990 ' ' INPUT PARAMETERS ' ---------------- ' N DIMENSION OF THE SYSTEM ' ' FCN NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE ' VALUE OF F(X,Y): ' Subroutine FCN(N,X,Y,F) ' dim Y(N), F(N) ' F(1)=... ETC. ' ' IFCN GIVES INFORMATION ON FCN: ' IFCN=0: F(X,Y) INDEPENDENT OF X (AUTONOMOUS) ' IFCN=1: F(X,Y) MAY DEPEND ON X (NON-AUTONOMOUS) ' ' X INITIAL X-VALUE ' ' Y(N) INITIAL VALUES FOR Y ' ' XEND FINAL X-VALUE (XEND-X MAY BE POSITIVE OR NEGATIVE) ' ' H INITIAL STEP SIZE GUESS; ' FOR STIFF EQUATIONS WITH INITIAL TRANSIENT, ' H=1.D0/(NORM OF F'), USUALLY 1.D-2 OR 1.D-3, IS GOOD. ' THIS CHOICE IS NOT VERY IMPORTANT, THE CODE QUICKLY ' ADAPTS ITS STEP SIZE. STUDY THE CHOSEN VALUES FOR A FEW ' STEPS IN SUBROUTINE "SOLOUT", WHEN YOU ARE NOT SURE. ' (IF H=0.D0, THE CODE PUTS H=1.D-6). ' ' RTOL,ATOL RELATIVE AND ABSOLUTE XERROR TOLERANCES. THEY ' CAN BE BOTH SCALARS OR ELSE BOTH VECTORS OF LENGTH N. ' ' ITOL SWITCH FOR RTOL AND ATOL: ' ITOL=0: BOTH RTOL AND ATOL ARE SCALARS. ' THE CODE KEEPS, ROUGHLY, THE LOCAL XERROR OF ' Y(I) BELOW RTOL*ABS(Y(I))+ATOL ' ITOL=1: BOTH RTOL AND ATOL ARE VECTORS. ' THE CODE KEEPS THE LOCAL XERROR OF Y(I) BELOW ' RTOL(I)*ABS(Y(I))+ATOL(I). ' ' JAC NAME (EXTERNAL) OF THE SUBROUTINE WHICH COMPUTES ' THE PARTIAL DERIVATIVES OF F(X,Y) WITH RESPECT TO Y ' (THIS ROUTINE IS ONLY CALLED IF IJAC=1). ' FOR IJAC=1, THIS SUBROUTINE MUST HAVE THE FORM: ' SUBROUTINE JAC(N,X,Y,DFY,LDFY) ' DIM Y(N),DFY(LDFY,N) ' DFY(1,1)= ... ' LDFY, THE COLUMN-LENGTH OF THE ARRAY, IS ' FURNISHED BY THE CALLING PROGRAM. ' IF MLJAC = N, THE JACOBIAN IS SUPPOSED TO ' BE FULL AND THE PARTIAL DERIVATIVES ARE ' STORED IN DFY AS ' DFY(I,J) = PARTIAL F(I) / PARTIAL Y(J) ' ELSE, THE JACOBIAN IS TAKEN AS BANDED AND ' THE PARTIAL DERIVATIVES ARE STORED ' DIAGONAL-WISE AS ' DFY(I-J+MUJAC+1,J) = PARTIAL F(I) / PARTIAL Y(J). ' ' IJAC SWITCH FOR THE COMPUTATION OF THE JACOBIAN: ' IJAC=0: JACOBIAN IS COMPUTED INTERNALLY BY FINITE ' DIFFERENCES, SUBROUTINE "JAC" IS NEVER CALLED. ' IJAC=1: JACOBIAN IS SUPPLIED BY SUBROUTINE JAC. ' ' MLJAC SWITCH FOR THE BANDED STRUCTURE OF THE JACOBIAN: ' MLJAC=N: JACOBIAN IS A FULL MATRIX. THE LINEAR ' ALGEBRA IS DONE BY FULL-MATRIX GAUSS-ELIMINATION. ' 0<=MLJAC<N: MLJAC IS THE LOWER BANDWITH OF JACOBIAN ' MATRIX (>= NUMBER OF NON-ZERO DIAGONALS BELOW ' THE MAIN DIAGONAL). ' ' MUJAC UPPER BANDWITH OF JACOBIAN MATRIX (>= NUMBER OF NON- ' ZERO DIAGONALS ABOVE THE MAIN DIAGONAL). ' NEED NOT BE DEFINED IF MLJAC=N. ' ' DFX NAME (EXTERNAL) OF THE SUBROUTINE WHICH COMPUTES ' THE PARTIAL DERIVATIVES OF F(X,Y) WITH RESPECT TO X ' (THIS ROUTINE IS ONLY CALLED IF IDFX=1 AND IFCN=1). ' OTHERWISE, THIS SUBROUTINE MUST HAVE THE FORM ' SUBROUTINE DFX(N,X,Y,FX) ' DIM Y(N),FX(N) ' FX(1)= ... ' ' IDFX SWITCH FOR THE COMPUTATION OF THE DF/DX: ' IDFX=0: DF/DX IS COMPUTED INTERNALLY BY FINITE ' DIFFERENCES, SUBROUTINE "DFX" IS NEVER CALLED. ' IDFX=1: DF/DX IS SUPPLIED BY SUBROUTINE DFX. ' ' ---- MAS,IMAS,MLMAS, AND MUMAS HAVE ANALOG MEANINGS ----- ' ---- FOR THE "MASS MATRIX" (THE MATRIX "M" OF SECTION IV.8): - ' ' MAS NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE MASS- ' MATRIX M. ' IF IMAS=0, THIS MATRIX IS ASSUMED TO BE THE IDENTITY ' MATRIX AND NEEDS NOT TO BE DEFINED; ' IF IMAS=1, THE SUBROUTINE MAS IS OF THE FORM ' SUBROUTINE MAS(N,AM,LMAS) ' DIM AM(LMAS,N) ' AM(1,1)= .... ' IF (MLMAS.EQ.N) THE MASS-MATRIX IS STORED ' AS FULL MATRIX LIKE ' AM(I,J) = M(I,J) ' ELSE, THE MATRIX IS TAKEN AS BANDED AND STORED ' DIAGONAL-WISE AS ' AM(I-J+MUMAS+1,J) = M(I,J). ' ' IMAS GIVES INFORMATION ON THE MASS-MATRIX: ' IMAS=0: M IS SUPPOSED TO BE THE IDENTITY ' MATRIX, MAS IS NEVER CALLED. ' IMAS=1: MASS-MATRIX IS SUPPLIED. ' ' MLMAS SWITCH FOR THE BANDED STRUCTURE OF THE MASS-MATRIX: ' MLMAS=N: THE FULL MATRIX CASE. THE LINEAR ' ALGEBRA IS DONE BY FULL-MATRIX GAUSS-ELIMINATION. ' 0<=MLMAS<N: MLMAS IS THE LOWER BANDWITH OF THE ' MATRIX (>= NUMBER OF NON-ZERO DIAGONALS BELOW ' THE MAIN DIAGONAL). ' MLMAS IS SUPPOSED TO BE <= MLJAC. ' ' MUMAS UPPER BANDWITH OF MASS-MATRIX (>= NUMBER OF NON- ' ZERO DIAGONALS ABOVE THE MAIN DIAGONAL). ' NEED NOT BE DEFINED IF MLMAS=N. ' MUMAS IS SUPPOSED TO BE <= MUJAC. ' ' SOLOUT NAME (EXTERNAL) OF SUBROUTINE PROVIDING THE ' NUMERICAL SOLUTION DURING INTEGRATION. ' IF IOUT=0, NO INTERPOLATION SUBROUTINE IS NECESSARY. ' IF IOUT=1, IT IS CALLED AFTER EVERY SUCCESSFUL STEP. ' IT MUST HAVE THE FORM ' SUBROUTINE SOLOUT (NR,XOLD,X,Y,N,IRTRN) ' DIM Y(N) ' .... ' SOLOUT FURNISHES THE SOLUTION "Y" AT THE NR-TH ' GRID-POINT "X" (THEREBY THE INITIAL VALUE IS ' THE FIRST GRID-POINT). ' "IRTRN" SERVES TO INTXERRUPT THE INTEGRATION. IF IRTRN ' IS SET <0, ROS4 RETURNS TO THE CALLING PROGRAM. ' ' IOUT GIVES INFORMATION ON THE SUBROUTINE SOLOUT: ' IOUT=0: SUBROUTINE IS NEVER CALLED ' IOUT=1: SUBROUTINE IS USED FOR OUTPUT ' ' WORK ARRAY OF WORKING SPACE OF LENGTH "LWORK". ' SERVES AS WORKING SPACE FOR ALL VECTORS AND MATRICES. ' "LWORK" MUST BE AT LEAST ' N*(LJAC+LMAS+LE1+8)+5 ' WHERE ' LJAC=N IF MLJAC=N (FULL JACOBIAN) ' LJAC=MLJAC+MUJAC+1 IF MLJAC<N (BANDED JAC.) ' AND ' LMAS=0 IF IMAS=0 ' LMAS=N IF IMAS=1 AND MLMAS=N (FULL) ' LMAS=MLMAS+MUMAS+1 IF MLMAS<N (BANDED MASS-M.) ' AND ' LE1=N IF MLJAC=N (FULL JACOBIAN) ' LE1=2*MLJAC+MUJAC+1 IF MLJAC<N (BANDED JAC.). ' ' IN THE USUAL CASE WHERE THE JACOBIAN IS FULL AND THE ' MASS-MATRIX IS THE INDENTITY (IMAS=0), THE MINIMUM ' STORAGE REQUIREMENT IS ' LWORK = 2*N*N+8*N+5. ' ' LWORK DECLARED LENGHT OF ARRAY "WORK". ' ' IWORK INTEGER WORKING SPACE OF LENGTH "LIWORK". ' "LIWORK" MUST BE AT LEAST N+2. ' ' LIWORK DECLARED LENGHT OF ARRAY "IWORK". ' ' ---------------------------------------------------------------------- ' ' SOPHISTICATED SETTING OF PARAMETERS ' ----------------------------------- ' SEVERAL PARAMETERS OF THE CODE ARE TUNED TO MAKE IT WORK ' WELL. THEY MAY BE DEFINED BY SETTING WORK(1),..,WORK(5) ' AS WELL AS IWORK(1),IWORK(2) DIFFERENT FROM ZERO. ' FOR ZERO INPUT, THE CODE CHOOSES DEFAULT VALUES: ' ' IWORK(1) THIS IS THE MAXIMAL NUMBER OF ALLOWED STEPS. ' THE DEFAULT VALUE (FOR IWORK(1)=0) IS 100000. ' ' IWORK(2) SWITCH FOR THE CHOICE OF THE COEFFICIENTS ' IF IWORK(2).EQ.1 METHOD OF SHAMPINE ' IF IWORK(2).EQ.2 METHOD GRK4T OF KAPS-RENTROP ' IF IWORK(2).EQ.3 METHOD GRK4A OF KAPS-RENTROP ' IF IWORK(2).EQ.4 METHOD OF VAN VELDHUIZEN (GAMMA=1/2) ' IF IWORK(2).EQ.5 METHOD OF VAN VELDHUIZEN ("D-STABLE") ' IF IWORK(2).EQ.6 AN L-STABLE METHOD ' THE DEFAULT VALUE (FOR IWORK(2)=0) IS IWORK(2)=2. ' ' WORK(1) UROUND, THE ROUNDING UNIT, DEFAULT 1D-16. ' ' WORK(2) MAXIMAL STEP SIZE, DEFAULT XEND-X. ' ' WORK(3), WORK(4) PARAMETERS FOR STEP SIZE SELECTION ' THE NEW STEP SIZE IS CHOSEN SUBJECT TO THE RESTRICTION ' WORK(3) <= HNEW/HOLD <= WORK(4) ' DEFAULT VALUES: WORK(3)=0.2D0, WORK(4)=6.D0 ' ' WORK(5) AVOID THE HUMP: AFTER TWO CONSECUTIVE STEP REJECTIONS ' THE STEP SIZE IS MULTIPLIED BY WORK(5) ' DEFAULT VALUES: WORK(5)=0.1# ' '----------------------------------------------------------------------- ' ' OUTPUT PARAMETERS ' ----------------- ' X X-VALUE WHERE THE SOLUTION IS COMPUTED ' (AFTER SUCCESSFUL RETURN X=XEND) ' ' Y(N) SOLUTION AT X ' ' H PREDICTED STEP SIZE OF THE LAST ACCEPTED STEP ' ' IDID REPORTS ON SUCCESSFULNESS UPON RETURN: ' IDID=1 COMPUTATION SUCCESSFUL, ' IDID=-1 COMPUTATION UNSUCCESSFUL. ' ' --------------------------------------------------------- ' *** *** *** *** *** *** *** *** *** *** *** *** *** ' DECLARATIONS ' *** *** *** *** *** *** *** *** *** *** *** *** *** ' IAUTNMS,IMPLCT,JBAND,IARRET: Boolean (here integers). ' -------------------------------------------------------------------- ' --- THESE COMMON VARIABLES CAN BE USED FOR STATISTICS: ' --- XNFCN NUMBER OF FUNCTION EVALUATIONS (THOSE FOR NUMERICAL ' EVALUATION OF THE JACOBIAN ARE NOT COUNTED) ' --- XNJAC NUMBER OF JACOBIAN EVALUATIONS (EITHER ANALYTICALLY ' OR NUMERICALLY) ' --- XNSTEP NUMBER OF COMPUTED STEPS ' --- XNACCPT NUMBER OF ACCEPTED STEPS ' --- NREJCT NUMBER OF REJECTED STEPS (AFTER AT LEAST ONE STEP ' HAS BEEN ACCEPTED) ' --- XNDEC NUMBER OF LU-DECOMPOSITIONS (N-DIMENSIONAL MATRIX) ' --- XNSOL NUMBER OF FORWARD-BACKWARD SUBSTITUTIONS ' -------------------------------------------------------------------- ' LDE,LDJAC,LDMAS,LDMAS2,METH: Integer ' XNMAX: LongInt (here double) ' FAC1,FAC2,FACREJ,HMAX,UROUND: Double ' IEYNEW,IEDY1,IEDY,IEAK1,IEAK2,IEAK3,IEAK4, ' IEFX,IEJAC,IEMAS, IEE, ISTORE: Integer; ' IEIP: Integer ' *** *** *** *** *** *** *** ' SETTING THE PARAMETERS ' *** *** *** *** *** *** *** XNFCN = ZERO XNJAC = ZERO XNSTEP = ZERO XNACCPT = ZERO NREJCT = 0 XNDEC = ZERO XNSOL = ZERO IARRET = 0 ' ------- XNMAX , THE MAXIMAL NUMBER OF STEPS ----- If IWORK(1) = 0 Then XNMAX = 100000# Else XNMAX = 1# * IWORK(1) If XNMAX <= ZERO Then Print " WRONG INPUT, IWORK(1)= "; IWORK(1) IARRET = 1 End If End If ' -------- METH COEFFICIENTS OF THE METHOD ------ If IWORK(2) = 0 Then METH = 2 Else METH = IWORK(2) If METH <= 0 Or METH >= 7 Then Print " CURIOUS INPUT, IWORK(2)="; IWORK(2) IARRET = 1 End If End If ' -------- UROUND, SMALLEST NUMBER SATISFYING ONE + UROUND > ONE --- If WORK(1) = ZERO Then UROUND = 1E-16 Else UROUND = WORK(1) If UROUND <= 0.00000000000001 Or UROUND >= ONE Then Print " COEFFICIENTS HAVE 16 DIGITS, UROUND="; WORK(1) IARRET = 1 End If End If ' -------- MAXIMAL STEP SIZE ----------- If WORK(2) = ZERO Then HMAX = XEND - X Else HMAX = WORK(2) End If ' ------- FAC1,FAC2 PARAMETERS FOR STEP SIZE SELECTION ------ If WORK(3) = ZERO Then FAC1 = 5# Else FAC1 = ONE / WORK(3) End If If WORK(4) = ZERO Then FAC2 = ONE / 6# Else FAC2 = ONE / WORK(4) End If ' ------- FACREJ FOR THE HUMP ------- If WORK(5) = ZERO Then FACREJ = 0.1 Else FACREJ = WORK(5) End If ' --------- CHECK IF TOLERANCES ARE O.K. --------- If ITOL = 0 Then If ATOL <= ZERO Or RTOL <= TEN * UROUND Then Print " TOLERANCES ARE TOO SMALL." IARRET = 1 End If Else 'Multiple tolerances not implemented here 'For I=1 to N 'IF ATOL(I) <= ZERO OR RTOL(I) <= TEN*UROUND THEN ' print " TOLERANCE(";I;") IS TOO SMALL." ' IARRET=1 'end if End If ' *** *** *** *** *** *** *** *** *** *** *** *** *** ' COMPUTATION OF ARRAY ENTRIES ' *** *** *** *** *** *** *** *** *** *** *** *** *** ' ---- AUTONOMOUS, IMPLICIT, BANDED OR NOT ? If IFCN = 0 Then IAUTNMS = 1 Else IAUTNMS = 0 End If If IMAS <> 0 Then IMPLCT = 1 Else IMPLCT = 0 End If If MLJAC <> N Then JBAND = 1 Else JBAND = 0 End If IARRET = 0 ' -------- COMPUTATION OF THE ROW-DIMENSIONS OF THE 2-ARRAYS ------ ' -- JACOBIAN -- If JBAND <> 0 Then LDJAC = MLJAC + MUJAC + 1 Else LDJAC = N End If ' -- MATRIX E FOR LINEAR ALGEBRA -- If JBAND <> 0 Then LDE = 2 * MLJAC + MUJAC + 1 Else LDE = N End If ' -- MASS MATRIX -- If IMPLCT <> 0 Then If MLMAS <> N Then LDMAS = MLMAS + MUMAS + 1 Else LDMAS = N End If ' ------ BANDWITH OF "MAS" NOT LARGER THAN BANDWITH OF "JAC" ------- If MLMAS > MLJAC Or MUMAS > MUJAC Then Print " BANDWITH OF MAS MUST NOT BE LARGER THAN BANDWITH OF JAC." IARRET = 1 End If Else LDMAS = 0 End If ia = 1: ib = LDMAS: GoSub 600 LDMAS2 = IMAX ' ------- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN WORK ------ IEYNEW = 6 IEDY1 = IEYNEW + N IEDY = IEDY1 + N IEAK1 = IEDY + N IEAK2 = IEAK1 + N IEAK3 = IEAK2 + N IEAK4 = IEAK3 + N IEFX = IEAK4 + N IEJAC = IEFX + N IEMAS = IEJAC + N * LDJAC IEE = IEMAS + N * LDMAS ' ------ TOTAL STORAGE REQUIREMENT ----------- ISTORE = IEE + N * LDE - 1 If ISTORE > LWORK Then Print " INSUFFICIENT STORAGE FOR WORK, MINIMUM LWORK="; ISTORE IARRET = 1 End If ' ------- ENTRY POINTS FOR INTEGER WORKSPACE ------ IEIP = 3 ' --------- TOTAL REQUIREMENT --------------- ISTORE = IEIP + N - 1 If ISTORE > LIWORK Then Print " INSUFF. STORAGE FOR IWORK, MINIMUM LIWORK="; ISTORE IARRET = 1 End If ' ------ WHEN A FAIL HAS OCCURED, WE RETURN WITH IDID=-1 ------ If IARRET <> 0 Then IDID = -1 Return End If ' ------------- prepare arguments of RO4COR -------------------------- ' Here, appropriate parts of WORK (IWORK) are put in temporary vectors ' TMP1 to TMP11 (ITMP) to simulate Fortran arguments of RO4COR. For I = 1 To N If I + IEYNEW - 1 <= N Then TMP1(I) = WORK(I + IEYNEW - 1) Else TMP1(I) = ZERO End If Next I For I = 1 To N If I + IEDY1 - 1 <= N Then TMP2(I) = WORK(I + IEDY1 - 1) Else TMP2(I) = ZERO End If Next I For I = 1 To N If I + IEDY - 1 <= N Then TMP3(I) = WORK(I + IEDY - 1) Else TMP3(I) = ZERO End If Next I For I = 1 To N If I + IEAK1 - 1 <= N Then TMP4(I) = WORK(I + IEAK1 - 1) Else TMP4(I) = ZERO End If Next I For I = 1 To N If I + IEAK2 - 1 <= N Then TMP5(I) = WORK(I + IEAK2 - 1) Else TMP5(I) = ZERO End If Next I For I = 1 To N If I + IEAK3 - 1 <= N Then TMP6(I) = WORK(I + IEAK3 - 1) Else TMP6(I) = ZERO End If Next I For I = 1 To N If I + IEAK4 - 1 <= N Then TMP7(I) = WORK(I + IEAK4 - 1) Else TMP7(I) = ZERO End If Next I For I = 1 To N If I + IEFX - 1 <= N Then TMP8(I) = WORK(I + IEFX - 1) Else TMP8(I) = ZERO End If Next I For I = 1 To N If I + IEJAC - 1 <= N Then TMP9(I) = WORK(I + IEJAC - 1) Else TMP9(I) = ZERO End If Next I For I = 1 To N If I + IEE - 1 <= N Then TMP10(I) = WORK(I + IEE - 1) Else TMP10(I) = ZERO End If Next I For I = 1 To No If I + IEMAS - 1 <= N Then TMP11(I) = WORK(I + IEMAS - 1) Else TMP11(I) = ZERO End If Next I For I = 1 To N If I + IEIP - 1 <= N Then ITMP(I) = IWORK(I + IEIP - 1) Else ITMP(I) = 0 End If Next I ' -------- CALL TO CORE INTEGRATOR ------------ GoSub 2000 'call RO4COR(N,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,IJAC, 'MLJAC,MUJAC,IDFX,MLMAS,MUMAS,IOUT,IDID, 'XNMAX,UROUND,METH,FAC1,FAC2,FACREJ,IAUTNMS, 'IMPLCT,JBAND,LDJAC,LDE,LDMAS2,TMP1,TMP2, 'TMP3,TMP4,TMP5,TMP6,TMP7,TMP8,TMP9,TMP10, 'TMP11, ITMP) Return 'from 1000 ROS4 ' --------- ... AND HERE IS THE CORE INTEGRATOR ---------- 2000 'Subroutine RO4COR(N,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,IJAC, 'MLJAC,MUJAC,IDFX,MLMAS,MUMAS,IOUT,IDID, 'XNMAX,UROUND,METH,FAC1,FAC2,FACREJ,IAUTNMS, 'IMPLCT,JBAND,LDJAC,LDE,LDMAS2,TMP1,TMP2, 'TMP3,TMP4,TMP5,TMP6,TMP7,TMP8, TMP9,TMP10, 'TMP11,ITMP) ' ---------------------------------------------------------- ' CORE INTEGRATOR FOR ROS4 ' PARAMETERS SAME AS IN ROS4 WITH WORKSPACE ADDED ' ---------------------------------------------------------- ' DECLARATIONS ' ---------------------------------------------------------- 'Labels 2001,2002, 2012, 2014, 2079, 2080 ' I,J,K,L,MBB,MBDIAG,MBJAC,MDIAG,MDIFF,MLDE,MUE: Integer ' IREJECT,IRJECT2: Boolean (here integers) ' A21,A31,A32,C21,C31,C32,C41,C42,C43: Double ' B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4: Double ' DELT,HMAXN,HOPT,POSNEG,XDELT,XXOLD,YSAFE: Double ' IRTRN,LBEG,LDEND,MD,MUJACJ,MUJACP,NSING:Integer; ' FAC,HC21,HC31,HC32,HC41,HC42,HC43,SUM: Double ' I1,I2,IB,INFO,J1,J2,MADD: Integer; ' XERR, HD1,HD2,HD3,HD4, HNEW,S,SK: Double ' ------- restore Fortran parameters FJAC, E, FMAS ---------- For J = 1 To N For I = 1 To N FJAC(I, J) = TMP9(I + J - 1) E(I, J) = TMP10(I + J - 1) FMAS(I, J) = TMP11(I + J - 1) Next I Next J ' ------- COMPUTE MASS MATRIX FOR IMPLICIT CASE ---------- ' IF IMPLCT <> 0 CALL MAS(N,FMAS,LDMAS2) ' (Not provided here). ' ---- PREPARE BANDWIDTHS ----- If JBAND <> 0 Then MLDE = MLJAC MUE = MUJAC MBJAC = MLJAC + MUJAC + 1 MBB = MLMAS + MUMAS + 1 MDIAG = MLDE + MUE + 1 MBDIAG = MUMAS + 1 MDIFF = MLDE + MUE - MUMAS End If ' *** *** *** *** *** *** *** ' INITIALISATIONS ' *** *** *** *** *** *** *** } 'POSNEG=SIGN(ONE,XEND-X) xa = ONE: xb = XEND - X: GoSub 650: POSNEG = SIGN If METH = 1 Then GoSub 3000 'call SHAMP(A21,A31,A32,C21,C31,C32,C41,C42,C43, 'B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4) If METH = 2 Then GoSub 3002 'call GRK4T(A21,A31,A32,C21,C31,C32,C41,C42,C43, 'B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4) If METH = 3 Then GoSub 3001 'call GRK4A(A21,A31,A32,C21,C31,C32,C41,C42,C43, 'B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4) If METH = 4 Then GoSub 3003 'call VELDS(A21,A31,A32,C21,C31,C32,C41,C42,C43, 'B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4) If METH = 5 Then GoSub 3004 'call VELDD(A21,A31,A32,C21,C31,C32,C41,C42,C43, 'B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4) If METH = 6 Then GoSub 3005 'call LSTAB(A21,A31,A32,C21,C31,C32,C41,C42,C43, 'B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4) ' --- INITIAL PREPARATIONS --- 'HMAXN=XMIN(ABS(HMAX),ABS(XEND-X)) xa = Abs(HMAX): xb = Abs(XEND - X): GoSub 630 HMAXN = XMIN 'H=XMIN(XMAX(1D-10,ABS(H)),HMAXN) xa = 0.0000000001: xb = Abs(H): GoSub 620 xa = XMAX: xb = HMAXN: GoSub 630: H = XMIN 'H=SIGN(H,POSNEG) xa = H: xb = POSNEG: GoSub 650: H = SIGN IREJECT = 0 NSING = 0 IRTRN = 1 XXOLD = X ' IF IOUT <> 0 THEN ' call SOLOUT(XNACCPT+1,XXOLD,X,Y,N,IRTRN) ' (not provided here). ' END IF If IRTRN < 0 Then GoTo 2079 ' --- BASIC INTEGRATION STEP --- 2001 If XNSTEP > XNMAX Or X + 0.1 * H = X Or Abs(H) <= UROUND Then GoTo 2079 If (X - XEND) * POSNEG + UROUND > ZERO Then H = HOPT IDID = 1 Return 'Normal return End If HOPT = H If (X + H - XEND) * POSNEG > ZERO Then H = XEND - X 'call FCN(N,X,Y,TMP2) XX = X For I = 1 To N YY(I) = Y(I) Next I GoSub 500 For I = 1 To N TMP2(I) = F(I) Next I XNFCN = XNFCN + 1# ' *** *** *** *** *** *** *** ' COMPUTATION OF THE JACOBIAN ' *** *** *** *** *** *** *** XNJAC = XNJAC + ONE If IJAC = 0 Then ' --- COMPUTE JACOBIAN MATRIX NUMERICALLY --- If JBAND <> 0 Then ' --- JACOBIAN IS BANDED --- MUJACP = MUJAC + 1 'MD=IMIN(MBJAC,N) ia = MBJAC: ib = N: GoSub 610: MD = IMIN For K = 1 To MD J = K 2012 TMP5(J) = Y(J) 'TMP6(J)=SQR(UROUND*XMAX(1D-5,ABS(Y(J)))) xa = 0.00001: xb = Y(J): GoSub 620 TMP6(J) = Sqr(UROUND * XMAX) Y(J) = Y(J) + TMP6(J) J = J + MD If J <= N Then GoTo 2012 'call FCN(N,X,Y,TMP4) GoSub 500 For I = 1 To N TMP4(I) = F(I) Next I J = K 'LBEG=IMAX(1,J-MUJAC) ia = 1: ib = MUJAC: GoSub 600: LBEG = IMAX 2014 'LDEND=IMIN(N,J+MLJAC) ia = N: ib = J + MLJAC: GoSub 610: LDEND = IMIN Y(J) = TMP5(J) MUJACJ = MUJACP - J For L = LBEG To LDEND FJAC(L + MUJACJ, J) = (TMP4(L) - TMP2(L)) / TMP6(J) Next L J = J + MD LBEG = LDEND + 1 If J <= N Then GoTo 2014 Next K Else ' --- JACOBIAN IS FULL --- For I = 1 To N YSAFE = Y(I) 'DELT=SQR(UROUND*XMAX(1D-5,ABS(YSAFE))) xa = 0.00001: xb = Abs(YSAFE): GoSub 620 DELT = Sqr(UROUND * XMAX) Y(I) = YSAFE + DELT 'call FCN(N,X,Y,TMP4) XX = X For II = 1 To N YY(II) = Y(II) Next II GoSub 500 For II = 1 To N TMP4(II) = F(II) Next II For J = 1 To N FJAC(J, I) = (TMP4(J) - TMP2(J)) / DELT Next J Y(I) = YSAFE Next I MLJAC = N End If Else ' --- COMPUTE JACOBIAN MATRIX ANALYTICALLY --- ' JAC(N,X,Y,FJAC,LDJAC) ' (Not provided here). End If If IAUTNMS = 0 Then If IDTMP8 = 0 Then ' --- COMPUTE NUMERICALLY THE DERIVATIVE WITH RESPECT TO X --- 'DELT=SQR(UROUND*XMAX(1D-5,ABS(X))) xa = 0.00001: xb = Abs(X): GoSub 620 DELT = Sqr(UROUND * XMAX) XDELT = X + DELT 'call FCN(N,XDELT,Y,TMP4) XX = XDELT: GoSub 500 For I = 1 To N TMP4(I) = F(I) Next I For J = 1 To N TMP8(J) = (TMP4(J) - TMP2(J)) / DELT Next J ' ELSE ' --- COMPUTE ANALYTICALLY THE DERIVATIVE WITH RESPECT TO X --- ' CALL DTMP8(N,X,Y,TMP8) ' (Not provided here). End If End If ' *** *** *** *** *** *** *** ' COMPUTE THE STAGES ' *** *** *** *** *** *** *** 2002 XNDEC = XNDEC + ONE HC21 = C21 / H HC31 = C31 / H HC32 = C32 / H HC41 = C41 / H HC42 = C42 / H HC43 = C43 / H FAC = ONE / (H * GAMMA) If IMPLCT <> 0 Then If JBAND <> 0 Then ' --- THE MATRIX E (B IS A JBAND MATRIX, JACOBIAN A JBAND MATRIX) --- For J = 1 To N 'I1=IMAX(1,MUJAC+2-J) ia = 1: ib = MUJAC + 2 - J: GoSub 600: I1 = IMAX 'I2=IMIN(MBJAC,N+MUJAC+1-J) ia = MBJAC: ib = N + MUJAC + 1 - J: GoSub 610: I2 = IMIN For I = I1 To I2 E(I + MLDE, J) = -FJAC(I, J) Next I Next J For J = 1 To N 'I1=IMAX(1,MUMAS+2-J) ia = 1: ib = MUMAS + 2 - J: GoSub 600: I1 = IMAX 'I2=IMIN(MBB,N+MUMAS+1-J) ia = MBB: ib = N + MUMAS + 1 - J: GoSub 610: I2 = IMIN For I = I1 To I2 ib = I + MDIFF E(ib, J) = E(ib, J) + FAC * FMAS(I, J) Next I Next J GoSub 4000 'call DECB(N,E,MLDE,MUE,ITMP,INFO) If INFO <> 0 Then GoTo 2080 If IAUTNMS <> 0 Then ' --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE ' --- 1) THE DIFFERENTIAL EQUATION IS IN IMPLICIT FORM ' --- 2) THE MATRIX B AND THE JACOBIAN OF F ARE BANDED ' --- 3) THE DIFFERENTIAL EQUATION IS AUTONOMOUS. For I = 1 To N TMP4(I) = TMP2(I) Next I For I = 1 To N B(I) = TMP4(I) Next I GoSub 5000 'call SOLB(N,E,MLDE,MUE,TMP4,ITMP) For I = 1 To N TMP4(I) = B(I) Next I For I = 1 To N TMP1(I) = Y(I) + A21 * TMP4(I) Next I 'call FCN(N,X,TMP1,TMP3) XX = X For I = 1 To N YY(I) = TMP1(I) Next I GoSub 500 For I = 1 To N TMP3(I) = F(I) Next I For I = 1 To N TMP1(I) = HC21 * TMP4(I) Next I For I = 1 To N Sum = ZERO 'J1=IMAX(1,I-MLMAS) ia = 1: ib = I - MLMAS: GoSub 600: J1 = IMAX 'J2=IMIN(N,I+MUMAS) ia = N: ib = I + MUMAS: GoSub 610: J2 = IMIN For J = J1 To J2 Sum = Sum + FMAS(I - J + MBDIAG, J) * TMP1(J) Next J TMP5(I) = Sum + TMP3(I) Next I For I = 1 To N B(I) = TMP5(I) Next I GoSub 5000 'call SOLB(N,E,MLDE,MUE,TMP5,ITMP) For I = 1 To N TMP5(I) = B(I) Next I For I = 1 To N TMP1(I) = Y(I) + A31 * TMP4(I) + A32 * TMP5(I) Next I 'call FCN(N,X,TMP1,TMP3) For I = 1 To N YY(I) = TMP1(I) Next I GoSub 500 For I = 1 To N TMP3(I) = F(I) Next I For I = 1 To N TMP1(I) = HC31 * TMP4(I) + HC32 * TMP5(I) Next I For I = 1 To N Sum = ZERO 'J1=IMAX(1,I-MLMAS) ia = 1: ib = I - MLMAS: GoSub 600: J1 = IMAX 'J2=IMIN(N,I+MUMAS) ia = N: ib = I + MUMAS: GoSub 610: J2 = IMIN For J = J1 To J2 Sum = Sum + FMAS(I - J + MBDIAG, J) * TMP1(J) Next J TMP6(I) = Sum + TMP3(I) Next I For I = 1 To N B(I) = TMP6(I) Next I GoSub 5000 'call SOLB(N,E,MLDE,MUE,TMP6,ITMP) For I = 1 To N TMP6(I) = B(I) Next I For I = 1 To N TMP1(I) = HC41 * TMP4(I) + HC42 * TMP5(I) + HC43 * TMP6(I) Next I For I = 1 To N Sum = ZERO 'J1=IMAX(1,I-MLMAS) ia = 1: ib = I - MLMAS: GoSub 600: J1 = IMAX 'J2=IMIN(N,I+MUMAS) ia = N: ib = I + MUMAS: GoSub 610: J2 = IMIN For J = J1 To J2 Sum = Sum + FMAS(I - J + MBDIAG, J) * TMP1(J) Next J TMP7(I) = Sum + TMP3(I) Next I For I = 1 To N B(I) = TMP7(I) Next I GoSub 5000 'call SOLB(N,E,MLDE,MUE,TMP7,ITMP) For I = 1 To N TMP7(I) = B(I) Next I Else ' --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE ' --- 1) THE DIFFERENTIAL EQUATION IS IN IMPLICIT FORM ' --- 2) THE MATRIX B AND THE JACOBIAN OF F ARE BANDED ' --- 3) THE DIFFERENTIAL EQUATION IS NON-AUTONOMOUS. HD1 = H * D1 HD2 = H * D2 HD3 = H * D3 HD4 = H * D4 For I = 1 To N TMP4(I) = TMP2(I) + HD1 * TMP8(I) Next I For I = 1 To N B(I) = TMP4(I) Next I GoSub 5000 'call SOLB(N,E,MLDE,MUE,TMP4,ITMP) For I = 1 To N TMP4(I) = B(I) Next I For I = 1 To N TMP1(I) = Y(I) + A21 * TMP4(I) Next I 'call FCN(N,X+C2*H,TMP1,TMP3) XX = X + C2 * H For I = 1 To N YY(I) = TMP1(I) Next I GoSub 500 For I = 1 To N TMP3(I) = F(I) Next I For I = 1 To N TMP1(I) = HC21 * TMP4(I) Next I For I = 1 To N Sum = ZERO 'J1=IMAX(1,I-MLMAS) ia = 1: ib = I - MLMAS: GoSub 600: J1 = IMAX 'J2=IMIN(N,I+MUMAS) ia = N: ib = I + MUMAS: GoSub 610: J2 = IMIN For J = J1 To J2 Sum = Sum + FMAS(I - J + MBDIAG, J) * TMP1(J) Next J TMP5(I) = Sum + TMP3(I) + HD2 * TMP8(I) Next I For I = 1 To N B(I) = TMP5(I) Next I GoSub 5000 'call SOLB(N,E,MLDE,MUE,TMP5,ITMP) For I = 1 To N TMP5(I) = B(I) Next I For I = 1 To N TMP1(I) = Y(I) + A31 * TMP4(I) + A32 * TMP5(I) Next I 'call FCN(N,X+C3*H,TMP1,TMP3) XX = X + C3 * H For I = 1 To N YY(I) = TMP1(I) Next I GoSub 500 For I = 1 To N TMP3(I) = F(I) Next I For I = 1 To N TMP1(I) = HC31 * TMP4(I) + HC32 * TMP5(I) Next I For I = 1 To N Sum = ZERO 'J1=IMAX(1,I-MLMAS) ia = 1: ib = I - MLMAS: GoSub 600: J1 = IMAX 'J2=IMIN(N,I+MUMAS) ia = N: ib = I + MUMAS: GoSub 610: J2 = IMIN For J = J1 To J2 Sum = Sum + FMAS(I - J + MBDIAG, J) * TMP1(J) Next J TMP6(I) = Sum + TMP3(I) + HD3 * TMP8(I) Next I For I = 1 To N B(I) = TMP6(I) Next I GoSub 5000 'call SOLB(N,E,MLDE,MUE,TMP6,ITMP) For I = 1 To N TMP6(I) = B(I) Next I For I = 1 To N TMP1(I) = HC41 * TMP4(I) + HC42 * TMP5(I) + HC43 * TMP6(I) Next I For I = 1 To N Sum = ZERO 'J1=IMAX(1,I-MLMAS) ia = 1: ib = I - MLMAS: GoSub 600: J1 = IMAX 'J2=IMIN(N,I+MUMAS) ia = N: ib = I + MUMAS: GoSub 610: J2 = IMIN For J = J1 To J2 Sum = Sum + FMAS(I - J + MBDIAG, J) * TMP1(J) Next J TMP7(I) = Sum + TMP3(I) + HD4 * TMP8(I) Next I For I = 1 To N B(I) = TMP7(I) Next I GoSub 5000 'call SOLB(N,E,MLDE,MUE,TMP7,ITMP) For I = 1 To N TMP7(I) = B(I) Next I End If Else If MLMAS <> N Then ' --- THE MATRIX E (B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX) --- MADD = MUMAS + 1 For J = 1 To N For I = 1 To N E(I, J) = -FJAC(I, J) Next I Next J For J = 1 To N 'I1=IMAX(1,J-MUMAS) ia = 1: ib = J - MUMAS: GoSub 600: I1 = IMAX 'I2=IMIN(N,J+MLMAS) ia = N: ib = J + MLMAS: GoSub 610: I2 = IMIN For I = I1 To I2 E(I, J) = E(I, J) + FAC * FMAS(I - J + MADD, J) Next I Next J GoSub 4500 'call DECA(N,LDE,E,ITMP,INFO) If INFO <> 0 Then GoTo 2080 If IAUTNMS <> 0 Then ' --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE ' --- 1) THE DIFFERENTIAL EQUATION IS IN IMPLICIT FORM ' --- 2) THE MATRIX B IS BANDED BUT THE JACOBIAN OF F IS NOT ' --- 3) THE DIFFERENTIAL EQUATION IS AUTONOMOUS. For I = 1 To N TMP4(I) = TMP2(I) Next I For I = 1 To N B(I) = TMP4(I) Next I GoSub 5500 'call SOL(N,E,TMP4,ITMP) For I = 1 To N TMP4(I) = B(I) Next I For I = 1 To N TMP1(I) = Y(I) + A21 * TMP4(I) Next I 'call FCN(N,X,TMP1,TMP3) XX = X For I = 1 To N YY(I) = TMP1(I) Next I GoSub 500 For I = 1 To N TMP3(I) = F(I) Next I For I = 1 To N TMP1(I) = HC21 * TMP4(I) Next I For I = 1 To N Sum = TMP3(I) 'J1=IMAX(1,I-MLMAS) ia = 1: ib = I - MLMAS: GoSub 600: J1 = IMAX 'J2=IMIN(N,I+MUMAS) ia = N: ib = I + MUMAS: GoSub 610: J2 = IMIN For J = J1 To J2 Sum = Sum + FMAS(I - J + MADD, J) * TMP1(J) Next J TMP5(I) = Sum Next I For I = 1 To N B(I) = TMP5(I) Next I GoSub 5500 'call SOL(N,E,TMP5,ITMP) For I = 1 To N TMP5(I) = B(I) Next I For I = 1 To N TMP1(I) = Y(I) + A31 * TMP4(I) + A32 * TMP5(I) Next I 'call FCN(N,X,TMP1,TMP3) For I = 1 To N YY(I) = TMP1(I) Next I GoSub 500 For I = 1 To N TMP3(I) = F(I) Next I For I = 1 To N TMP1(I) = HC31 * TMP4(I) + HC32 * TMP5(I) Next I For I = 1 To N Sum = TMP3(I) 'J1=IMAX(1,I-MLMAS) ia = 1: ib = I - MLMAS: GoSub 600: J1 = IMAX 'J2=IMIN(N,I+MUMAS) ia = N: ib = I + MUMAS: GoSub 610: J2 = IMIN For J = J1 To J2 Sum = Sum + FMAS(I - J + MADD, J) * TMP1(J) Next J TMP6(I) = Sum Next I For I = 1 To N B(I) = TMP6(I) Next I GoSub 5500 'call SOL(N,E,TMP6,ITMP) For I = 1 To N TMP6(I) = B(I) Next I For I = 1 To N TMP1(I) = HC41 * TMP4(I) + HC42 * TMP5(I) + HC43 * TMP6(I) Next I For I = 1 To N Sum = TMP3(I) 'J1=IMAX(1,I-MLMAS) ia = 1: ib = I - MLMAS: GoSub 600: J1 = IMAX 'J2=IMIN(N,I+MUMAS) ia = N: ib = I + MUMAS: GoSub 610: J2 = IMIN For J = J1 To J2 Sum = Sum + FMAS(I - J + MADD, J) * TMP1(J) Next J TMP7(I) = Sum Next I For I = 1 To N B(I) = TMP7(I) Next I GoSub 5500 'call SOL(N,E,TMP7,ITMP) For I = 1 To N TMP7(I) = B(I) Next I Else ' --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE ' --- 1) THE DIFFERENTIAL EQUATION IS IN IMPLICIT FORM ' --- 2) THE MATRIX B IS BANDED BUT THE JACOBIAN OF F IS NOT ' --- 3) THE DIFFERENTIAL EQUATION IS NON-AUTONOMOUS. HD1 = H * D1 HD2 = H * D2 HD3 = H * D3 HD4 = H * D4 For I = 1 To N TMP4(I) = TMP2(I) + HD1 * TMP8(I) Next I For I = 1 To N B(I) = TMP4(I) Next I GoSub 5500 'call SOL(N,E,TMP4,ITMP) For I = 1 To N TMP4(I) = B(I) Next I For I = 1 To N TMP1(I) = Y(I) + A21 * TMP4(I) Next I 'call FCN(N,X+C2*H,TMP1,TMP3) XX = X + C2 * H For I = 1 To N YY(I) = TMP1(I) Next I GoSub 500 For I = 1 To N TMP3(I) = F(I) Next I For I = 1 To N TMP1(I) = HC21 * TMP4(I) Next I For I = 1 To N Sum = TMP3(I) + HD2 * TMP8(I) 'J1=IMAX(1,I-MLMAS) ia = 1: ib = I - MLMAS: GoSub 600: J1 = IMAX 'J2=IMIN(N,I+MUMAS) ia = N: ib = I + MUMAS: GoSub 610: J2 = IMIN For J = J1 To J2 Sum = Sum + FMAS(I - J + MADD, J) * TMP1(J) Next J TMP5(I) = Sum Next I For I = 1 To N B(I) = TMP5(I) Next I GoSub 5500 'call SOL(N,E,TMP5,ITMP) For I = 1 To N TMP5(I) = B(I) Next I For I = 1 To N TMP1(I) = Y(I) + A31 * TMP4(I) + A32 * TMP5(I) Next I 'call FCN(N,X+C3*H,TMP1,TMP3) XX = X + C3 * H For I = 1 To N YY(I) = TMP1(I) Next I GoSub 500 For I = 1 To N TMP3(I) = F(I) Next I For I = 1 To N TMP1(I) = HC31 * TMP4(I) + HC32 * TMP5(I) Next I For I = 1 To N Sum = TMP3(I) + HD3 * TMP8(I) 'J1=IMAX(1,I-MLMAS) ia = 1: ib = I - MLMAS: GoSub 600: J1 = IMAX 'J2=IMIN(N,I+MUMAS) ia = N: ib = I + MUMAS: GoSub 610: J2 = IMIN For J = J1 To J2 Sum = Sum + FMAS(I - J + MADD, J) * TMP1(J) Next J TMP6(I) = Sum Next I For I = 1 To N B(I) = TMP5(I) Next I GoSub 5500 'call SOL(N,E,TMP5,ITMP) For I = 1 To N TMP5(I) = B(I) Next I For I = 1 To N TMP1(I) = HC41 * TMP4(I) + HC42 * TMP5(I) + HC43 * TMP6(I) Next I For I = 1 To N Sum = TMP3(I) + HD4 * TMP8(I) 'J1=IMAX(1,I-MLMAS) ia = 1: ib = I - MLMAS: GoSub 600: J1 = IMAX 'J2=IMIN(N,I+MUMAS) ia = N: ib = I + MUMAS: GoSub 610: J2 = IMIN For J = J1 To J2 Sum = Sum + FMAS(I - J + MADD, J) * TMP1(J) Next J TMP7(I) = Sum Next I For I = 1 To N B(I) = TMP7(I) Next I GoSub 5500 'call SOL(N,E,TMP7,ITMP) For I = 1 To N TMP7(I) = B(I) Next I End If Else ' --- THE MATRIX E (B IS A FULL MATRIX, JACOBIAN A FULL OR BANDED MATRIX) --- If MLJAC = N Then For J = 1 To N For I = 1 To N E(I, J) = FMAS(I, J) * FAC - FJAC(I, J) Next I Next J Else MADD = MUJAC + 1 For J = 1 To N For I = 1 To N E(I, J) = FMAS(I, J) * FAC Next I Next J For J = 1 To N 'I1=IMAX(1,J-MUJAC) ia = 1: ib = J - MUJAC: GoSub 600: I1 = IMAX 'I2=IMIN(N,J+MLJAC) ia = N: ib = J + MLJAC: GoSub 610: I2 = IMIN For I = I1 To I2 E(I, J) = E(I, J) - FJAC(I - J + MADD, J) Next I Next J End If GoSub 4500 'call DECA(N,E,ITMP,INFO) If INFO <> 0 Then GoTo 2080 If IAUTNMS <> 0 Then ' --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE ' --- 1) THE DIFFERENTIAL EQUATION IS IN IMPLICIT FORM ' --- 2) THE MATRIX B IS NOT BANDED ' --- 3) THE DIFFERENTIAL EQUATION IS AUTONOMOUS. For I = 1 To N TMP4(I) = TMP2(I) Next I For I = 1 To N B(I) = TMP4(I) Next I GoSub 5500 'call SOL(N,E,TMP4,ITMP) For I = 1 To N TMP4(I) = B(I) Next I For I = 1 To N TMP1(I) = Y(I) + A21 * TMP4(I) Next I 'call FCN(N,X,TMP1,TMP3) XX = X For I = 1 To N YY(I) = TMP1(I) Next I GoSub 500 For I = 1 To N TMP3(I) = F(I) Next I For I = 1 To N TMP1(I) = HC21 * TMP4(I) Next I For I = 1 To N Sum = TMP3(I) For J = 1 To N Sum = Sum + FMAS(I, J) * TMP1(J) Next J TMP5(I) = Sum Next I For I = 1 To N B(I) = TMP5(I) Next I GoSub 5500 'call SOL(N,E,TMP5,ITMP) For I = 1 To N TMP5(I) = B(I) Next I For I = 1 To N TMP1(I) = Y(I) + A31 * TMP4(I) + A32 * TMP5(I) Next I 'call FCN(N,X,TMP1,TMP3) For I = 1 To N YY(I) = TMP1(I) Next I GoSub 500 For I = 1 To N TMP3(I) = F(I) Next I For I = 1 To N TMP1(I) = HC31 * TMP4(I) + HC32 * TMP5(I) Next I For I = 1 To N Sum = TMP3(I) For J = 1 To N Sum = Sum + FMAS(I, J) * TMP1(J) Next J TMP6(I) = Sum Next I For I = 1 To N B(I) = TMP6(I) Next I GoSub 5500 'call SOL(N,E,TMP6,ITMP) For I = 1 To N TMP6(I) = B(I) Next I For I = 1 To N TMP1(I) = HC41 * TMP4(I) + HC42 * TMP5(I) + HC43 * TMP6(I) Next I For I = 1 To N Sum = TMP3(I) For J = 1 To N Sum = Sum + FMAS(I, J) * TMP1(J) Next J TMP7(I) = Sum Next I For I = 1 To N B(I) = TMP7(I) Next I GoSub 5500 'call SOL(N,E,TMP7,ITMP) For I = 1 To N TMP7(I) = B(I) Next I Else ' --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE ' --- 1) THE DIFFERENTIAL EQUATION IS IN IMPLICIT FORM ' --- 2) THE MATRIX B IS NOT BANDED ' --- 3) THE DIFFERENTIAL EQUATION IS NON-AUTONOMOUS. HD1 = H * D1 HD2 = H * D2 HD3 = H * D3 HD4 = H * D4 For I = 1 To N TMP4(I) = TMP2(I) + HD1 * TMP8(I) Next I For I = 1 To N B(I) = TMP4(I) Next I GoSub 5500 'call SOL(N,E,TMP4,ITMP) For I = 1 To N TMP4(I) = B(I) Next I For I = 1 To N TMP1(I) = Y(I) + A21 * TMP4(I) Next I 'call FCN(N,X+C2*H,TMP1,TMP3) XX = X + C2 * H For I = 1 To N YY(I) = TMP1(I) Next I GoSub 500 For I = 1 To N TMP3(I) = F(I) Next I For I = 1 To N TMP1(I) = HC21 * TMP4(I) Next I For I = 1 To N Sum = TMP3(I) + HD2 * TMP8(I) For J = 1 To N Sum = Sum + FMAS(I, J) * TMP1(J) Next J TMP5(I) = Sum Next I For I = 1 To N B(I) = TMP5(I) Next I GoSub 5500 'call SOL(N,E,TMP5,ITMP) For I = 1 To N TMP5(I) = B(I) Next I For I = 1 To N TMP1(I) = Y(I) + A31 * TMP4(I) + A32 * TMP5(I) Next I 'call FCN(N,X+C3*H,TMP1,TMP3) XX = X + C3 * H For I = 1 To N YY(I) = TMP1(I) Next I GoSub 500 For I = 1 To N TMP3(I) = F(I) Next I For I = 1 To N TMP1(I) = HC31 * TMP4(I) + HC32 * TMP5(I) Next I For I = 1 To N Sum = TMP3(I) + HD3 * TMP8(I) For J = 1 To N Sum = Sum + FMAS(I, J) * TMP1(J) Next J TMP6(I) = Sum Next I For I = 1 To N B(I) = TMP6(I) Next I GoSub 5500 'call SOL(N,E,TMP6,ITMP) For I = 1 To N TMP6(I) = B(I) Next I For I = 1 To N TMP1(I) = HC41 * TMP4(I) + HC42 * TMP5(I) + HC43 * TMP6(I) Next I For I = 1 To N Sum = TMP3(I) + HD4 * TMP8(I) For J = 1 To N Sum = Sum + FMAS(I, J) * TMP1(J) Next J TMP7(I) = Sum Next I For I = 1 To N B(I) = TMP7(I) Next I GoSub 5500 'call SOL(N,E,TMP7,ITMP) For I = 1 To N TMP7(I) = B(I) Next I End If End If End If Else If JBAND <> 0 Then ' --- THE MATRIX E (B=IDENTITY, JACOBIAN A BANDED MATRIX) --- For J = 1 To N 'I1=IMAX(1,MUJAC+2-J) ia = 1: ib = MUJAC + 2 - J: GoSub 600: I1 = IMAX 'I2=IMIN(MBJAC,N+MUJAC+1-J) ia = MBJAC: ib = N + MUJAC + 1 - J: GoSub 610: I2 = IMIN For I = I1 To I2 E(I + MLDE, J) = -FJAC(I, J) Next I E(MDIAG, J) = E(MDIAG, J) + FAC Next J GoSub 4000 'call DECB(N,E,MLDE,MUE,ITMP,INFO) If INFO <> 0 Then GoTo 2080 If IAUTNMS <> 0 Then ' --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE ' --- 1) THE DIFFERENTIAL EQUATION IS IN EXPLICIT FORM ' --- 2) THE JACOBIAN OF THE PROBLEM IS A BANDED MATRIX ' --- 3) THE DIFFERENTIAL EQUATION IS AUTONOMOUS. For I = 1 To N TMP4(I) = TMP2(I) Next I For I = 1 To N B(I) = TMP4(I) Next I GoSub 5000 'call SOLB(N,E,MLDE,MUE,TMP4,ITMP) For I = 1 To N TMP4(I) = B(I) Next I For I = 1 To N TMP1(I) = Y(I) + A21 * TMP4(I) Next I 'call FCN(N,X,TMP1,TMP3) XX = X For I = 1 To N YY(I) = TMP1(I) Next I GoSub 500 For I = 1 To N TMP3(I) = F(I) Next I For I = 1 To N TMP5(I) = TMP3(I) + HC21 * TMP4(I) Next I For I = 1 To N B(I) = TMP5(I) Next I GoSub 5000 'call SOLB(N,E,MLDE,MUE,TMP5,ITMP) For I = 1 To N TMP5(I) = B(I) Next I For I = 1 To N TMP1(I) = Y(I) + A31 * TMP4(I) + A32 * TMP5(I) Next I 'call FCN(N,X,TMP1,TMP3) For I = 1 To N YY(I) = TMP1(I) Next I GoSub 500 For I = 1 To N TMP3(I) = F(I) Next I For I = 1 To N TMP6(I) = TMP3(I) + HC31 * TMP4(I) + HC32 * TMP5(I) Next I For I = 1 To N B(I) = TMP6(I) Next I GoSub 5000 'call SOLB(N,E,MLDE,MUE,TMP6,ITMP) For I = 1 To N TMP6(I) = B(I) Next I For I = 1 To N TMP7(I) = TMP3(I) + HC41 * TMP4(I) + HC42 * TMP5(I) + HC43 * TMP6(I) Next I For I = 1 To N B(I) = TMP7(I) Next I GoSub 5000 'call SOLB(N,E,MLDE,MUE,TMP7,ITMP) For I = 1 To N TMP7(I) = B(I) Next I Else ' --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE ' --- 1) THE DIFFERENTIAL EQUATION IS IN EXPLICIT FORM ' --- 2) THE JACOBIAN OF THE PROBLEM IS A BANDED MATRIX ' --- 3) THE DIFFERENTIAL EQUATION IS NON-AUTONOMOUS. HD1 = H * D1 HD2 = H * D2 HD3 = H * D3 HD4 = H * D4 For I = 1 To N TMP4(I) = TMP2(I) + HD1 * TMP8(I) Next I For I = 1 To N B(I) = TMP4(I) Next I GoSub 5000 'call SOLB(N,E,MLDE,MUE,TMP4,ITMP) For I = 1 To N TMP4(I) = B(I) Next I For I = 1 To N TMP1(I) = Y(I) + A21 * TMP4(I) Next I 'call FCN(N,X+C2*H,TMP1,TMP3) XX = X + C2 * H For I = 1 To N YY(I) = TMP1(I) Next I GoSub 500 For I = 1 To N TMP3(I) = F(I) Next I For I = 1 To N TMP5(I) = TMP3(I) + HD2 * TMP8(I) + HC21 * TMP4(I) Next I For I = 1 To N B(I) = TMP5(I) Next I GoSub 5000 'call SOLB(N,E,MLDE,MUE,TMP5,ITMP) For I = 1 To N TMP5(I) = B(I) Next I For I = 1 To N TMP1(I) = Y(I) + A31 * TMP4(I) + A32 * TMP5(I) Next I 'call FCN(N,X+C3*H,TMP1,TMP3) XX = X + C3 * H For I = 1 To N YY(I) = TMP1(I) Next I GoSub 500 For I = 1 To N TMP3(I) = F(I) Next I For I = 1 To N TMP6(I) = TMP3(I) + HD3 * TMP8(I) + HC31 * TMP4(I) + HC32 * TMP5(I) Next I For I = 1 To N B(I) = TMP6(I) Next I GoSub 5000 'call SOLB(N,E,MLDE,MUE,TMP6,ITMP) For I = 1 To N TMP6(I) = B(I) Next I For I = 1 To N TMP7(I) = TMP3(I) + HD4 * TMP8(I) + HC41 * TMP4(I) + HC42 * TMP5(I) + HC43 * TMP6(I) Next I For I = 1 To N B(I) = TMP7(I) Next I GoSub 5000 'call SOLB(N,E,MLDE,MUE,TMP7,ITMP) For I = 1 To N TMP7(I) = B(I) Next I End If Else ' --- THE MATRIX E (B=IDENTITY, JACOBIAN A FULL MATRIX) --- For J = 1 To N For I = 1 To N E(I, J) = -FJAC(I, J) Next I E(J, J) = E(J, J) + FAC Next J GoSub 4500 'call DECA(N,LDE,E,ITMP,INFO) If INFO <> 0 Then GoTo 2080 If IAUTNMS <> 0 Then ' --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE ' --- 1) THE DIFFERENTIAL EQUATION IS IN EXPLICIT FORM ' --- 2) THE JACOBIAN OF THE PROBLEM IS A FULL MATRIX ' --- 3) THE DIFFERENTIAL EQUATION IS AUTONOMOUS. For I = 1 To N TMP4(I) = TMP2(I) Next I For I = 1 To N B(I) = TMP4(I) Next I GoSub 5500 'call SOL(N,E,TMP4,ITMP) For I = 1 To N TMP4(I) = B(I) Next I For I = 1 To N TMP1(I) = Y(I) + A21 * TMP4(I) Next I 'call FCN(N,X,TMP1,TMP3) XX = X For I = 1 To N YY(I) = TMP1(I) Next I GoSub 500 For I = 1 To N TMP3(I) = F(I) Next I For I = 1 To N TMP5(I) = TMP3(I) + HC21 * TMP4(I) Next I For I = 1 To N B(I) = TMP5(I) Next I GoSub 5500 'call SOL(N,E,TMP5,ITMP) For I = 1 To N TMP5(I) = B(I) Next I For I = 1 To N TMP1(I) = Y(I) + A31 * TMP4(I) + A32 * TMP5(I) Next I 'call FCN(N,X,TMP1,TMP3) For I = 1 To N YY(I) = TMP1(I) Next I GoSub 500 For I = 1 To N TMP3(I) = F(I) Next I For I = 1 To N TMP6(I) = TMP3(I) + HC31 * TMP4(I) + HC32 * TMP5(I) Next I For I = 1 To N B(I) = TMP6(I) Next I GoSub 5500 'call SOL(N,E,TMP6,ITMP) For I = 1 To N TMP6(I) = B(I) Next I For I = 1 To N TMP7(I) = TMP3(I) + HC41 * TMP4(I) + HC42 * TMP5(I) + HC43 * TMP6(I) Next I For I = 1 To N B(I) = TMP7(I) Next I GoSub 5500 'call SOL(N,E,TMP7,ITMP) For I = 1 To N TMP7(I) = B(I) Next I Else ' --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE ' --- 1) THE DIFFERENTIAL EQUATION IS IN EXPLICIT FORM ' --- 2) THE JACOBIAN OF THE PROBLEM IS A FULL MATRIX ' --- 3) THE DIFFERENTIAL EQUATION IS NON-AUTONOMOUS. HD1 = H * D1 HD2 = H * D2 HD3 = H * D3 HD4 = H * D4 For I = 1 To N TMP4(I) = TMP2(I) + HD1 * TMP8(I) Next I For I = 1 To N B(I) = TMP4(I) Next I GoSub 5500 'call SOL(N,E,TMP4,ITMP) For I = 1 To N TMP4(I) = B(I) Next I For I = 1 To N TMP1(I) = Y(I) + A21 * TMP4(I) Next I 'call FCN(N,X+C2*H,TMP1,TMP3) XX = X + C2 * H For I = 1 To N YY(I) = TMP1(I) Next I GoSub 500 For I = 1 To N TMP3(I) = F(I) Next I For I = 1 To N TMP5(I) = TMP3(I) + HD2 * TMP8(I) + HC21 * TMP4(I) Next I For I = 1 To N B(I) = TMP5(I) Next I GoSub 5500 'call SOL(N,E,TMP5,ITMP) For I = 1 To N TMP5(I) = B(I) Next I For I = 1 To N TMP1(I) = Y(I) + A31 * TMP4(I) + A32 * TMP5(I) Next I 'call FCN(N,X+C3*H,TMP1,TMP3) XX = X + C3 * H For I = 1 To N YY(I) = TMP1(I) Next I GoSub 500 For I = 1 To N TMP3(I) = F(I) Next I For I = 1 To N TMP6(I) = TMP3(I) + HD3 * TMP8(I) + HC31 * TMP4(I) + HC32 * TMP5(I) Next I For I = 1 To N B(I) = TMP6(I) Next I GoSub 5500 'call SOL(N,E,TMP6,ITMP) For I = 1 To N TMP6(I) = B(I) Next I For I = 1 To N TMP7(I) = TMP3(I) + HD4 * TMP8(I) + HC41 * TMP4(I) + HC42 * TMP5(I) + HC43 * TMP6(I) Next I For I = 1 To N B(I) = TMP7(I) Next I GoSub 5500 'call SOL(N,E,TMP7,ITMP) For I = 1 To N TMP7(I) = B(I) Next I End If End If End If XNSOL = XNSOL + 4# XNFCN = XNFCN + 2# ' *** *** *** *** *** *** *** ' ERROR ESTIMATION ' *** *** *** *** *** *** *** XNSTEP = XNSTEP + ONE ' ------------ NEW SOLUTION --------------- For I = 1 To N TMP1(I) = Y(I) + B1 * TMP4(I) + B2 * TMP5(I) + B3 * TMP6(I) + B4 * TMP7(I) Next I ' ------------ COMPUTE ERROR ESTIMATION ---------------- XERR = ZERO For I = 1 To N S = E1 * TMP4(I) + E2 * TMP5(I) + E3 * TMP6(I) + E4 * TMP7(I) If ITOL = 0 Then 'SK = ATOL + RTOL * XMAX(ABS(Y(I)), ABS(TMP1(I))) xa = Abs(Y(I)): xb = Abs(TMP1(I)): GoSub 620 SK = ATOL + RTOL * XMAX Else ' Multiple tolerances not implemented here. ' SK = ATOL(I) + RTOL(I) * XMAX(ABS(Y(I)), ABS(TMP1(I))) End If XERR = XERR + Sqr(Abs(S / SK)) Next I XERR = Sqr(XERR / N) ' --- COMPUTATION OF HNEW ' --- WE REQUIRE 0.2<=HNEW/H<=6.0 'FAC=XMAX(FAC2,XMIN(FAC1,XERR^(0.25#/0.9#))) xa = FAC1: xb = XERR ^ (0.25 / 0.9): GoSub 630 xa = FAC2: xb = XMIN: GoSub 620: FAC = XMAX HNEW = H / FAC ' *** *** *** *** *** *** *** ' IS THE XERROR SMALL ENOUGH ? ' *** *** *** *** *** *** *** If XERR <= ONE Then ' --- STEP IS ACCEPTED --- XNACCPT = XNACCPT + ONE For I = 1 To N Y(I) = TMP1(I) Next I XXOLD = X X = X + H If IOUT <> 0 Then 'SOLOUT(XNACCPT+1,XXOLD,X,Y,N,IRTRN) End If If IRTRN < 0 Then GoTo 2079 If Abs(HNEW) > HMAXN Then HNEW = POSNEG * HMAXN If IREJECT <> 0 Then 'HNEW=POSNEG*XMIN(ABS(HNEW),ABS(H)) xa = Abs(HNEW): xb = Abs(H): GoSub 630 HNEW = POSNEG * XMIN End If IREJECT = 0 IRJECT2 = 0 H = HNEW GoTo 2001 Else ' --- STEP IS IREJECTED --- If IRJECT2 <> 0 Then HNEW = H * FACREJ If IREJECT <> 0 Then IRJECT2 = 1 IREJECT = 1 H = HNEW If XNACCPT >= ONE Then NREJCT = NREJCT + 1 GoTo 2002 End If ' --- EXIT SECTION --- 2080 Print " MATRIX E IS SINGULAR, INFO = "; INFO NSING = NSING + 1 If NSING >= 5 Then GoTo 2079 H = H * HALF GoTo 2002 2079 Print Print " EXIT OF ROS4 AT X="; X; " H="; H IDID = -1 Return 'from RO4COR 3000 'SHAMP(A21,A31,A32,C21,C31,C32,C41,C42,C43,B1,B2,B3,B4, 'E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4) A21 = 2# A31 = 48# / 25# A32 = 6# / 25# C21 = -8# C31 = 372# / 25# C32 = 12# / 5# C41 = -112# / 125# C42 = -54# / 125# C43 = -2# / 5# B1 = 19# / 9# B2 = 1# / 2# B3 = 25# / 108# B4 = 125# / 108# E1 = 17# / 54# E2 = 7# / 36# E3 = 0# E4 = 125# / 108# GAMMA = HALF C2 = 1# C3 = 0.6 D1 = 0.5 D2 = -1.5 D3 = 2.42 D4 = 0.116 Return 3001 'GRK4A(A21,A31,A32,C21,C31,C32,C41,C42,C43,B1,B2,B3,B4, 'E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4) A21 = 1.10886075949367 A31 = 2.37708526198336 A32 = 0.185011498889969 C21 = -4.92018840239764 C31 = 1.05558868604858 C32 = 3.35181726766894 C41 = 3.84686900704931 C42 = 3.42710924126818 C43 = -2.16240884875326 B1 = 1.84568324040584 B2 = 0.13697968943605 B3 = 0.712909778329156 B4 = 0.632911392405063 E1 = 4.83187017720177E-02 E2 = -0.647110865104951 E3 = 0.218687666050024 E4 = -0.632911392405063 GAMMA = 0.395 C2 = 0.438 C3 = 0.87 D1 = 0.395 D2 = -0.372672395484092 D3 = 6.62919654457149E-02 D4 = 0.434094696256863 Return 3002 'GRK4T(A21,A31,A32,C21,C31,C32,C41,C42,C43,B1,B2,B3,B4, 'E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4) A21 = 2# A31 = 4.52470820737312 A32 = 4.16352878859765 C21 = -5.07167533877632 C31 = 6.02015272865079 C32 = 0.159750684672712 C41 = -1.85634361868611 C42 = -8.50538085817983 C43 = -2.08407513602319 B1 = 3.95750374664078 B2 = 4.62489238836331 B3 = 0.617477263875011 B4 = 1.28261294526904 E1 = 2.302155402933 E2 = 3.07363448539262 E3 = -0.873280801804503 E4 = -1.28261294526904 GAMMA = 0.231 C2 = 0.462 C3 = 0.880208333333333 D1 = 0.231 D2 = -0.039629667752443 D3 = 0.550778939578913 D4 = -5.53509845705276E-02 Return 3003 'VELDS(A21,A31,A32,C21,C31,C32,C41,C42,C43,B1,B2,B3,B4, 'E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4) ' --- METHOD GIVEN BY VAN VELDHUIZEN --- A21 = 2# A31 = 1.75 A32 = 0.25 C21 = -8# C31 = -8# C32 = -1# C41 = 0.5 C42 = -0.5 C43 = 2# B1 = 1.33333333333333 B2 = 0.666666666666667 B3 = -1.33333333333333 B4 = 1.33333333333333 E1 = -0.333333333333333 E2 = -0.333333333333333 E3 = -0# E4 = -1.33333333333333 GAMMA = 0.5 C2 = 1# C3 = 0.5 D1 = 0.5 D2 = -1.5 D3 = -0.75 D4 = 0.25 Return 3004 'VELDD(A21,A31,A32,C21,C31,C32,C41,C42,C43,B1,B2,B3,B4, 'E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4) ' --- METHOD GIVEN BY VAN VELDHUIZEN --- A21 = 2# A31 = 4.81223436269544 A32 = 4.57814695674784 C21 = -5.33333333333333 C31 = 6.10052967884825 C32 = 1.80473679737843 C41 = -2.54051545663475 C42 = -9.44374632891521 C43 = -1.98847175321599 B1 = 4.28933925465454 B2 = 5.03609848285141 B3 = 0.608573642067392 B4 = 1.35595894120115 E1 = 2.17567278753176 E2 = 2.95091122257574 E3 = -0.785974454488743 E4 = -1.35595894120115 GAMMA = 0.225708114822568 C2 = 0.451416229645136 C3 = 0.875592894601846 D1 = 0.225708114822568 D2 = -4.59940350268058E-02 D3 = 0.517759050494408 D4 = -3.80562393805443E-02 Return 3005 'LSTAB(A21,A31,A32,C21,C31,C32,C41,C42,C43,B1,B2,B3,B4, 'E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4) ' --- AN L-STABLDE METHOD --- A21 = 2# A31 = 1.86794363780392 A32 = 0.234444971139916 C21 = -7.13761503641231 C31 = 2.58070808795146 C32 = 0.651595007644798 C41 = -2.13714899438253 C42 = -0.321466969123763 C43 = -0.694974250178178 B1 = 2.25557007341874 B2 = 0.287049326218679 B3 = 0.435317943184018 B4 = 1.09350225240916 E1 = -0.281543193214115 E2 = -7.27619912493892E-02 E3 = -0.108219620149531 E4 = -1.09350225240916 GAMMA = 0.57282 C2 = 1.14564 C3 = 0.65521686381559 D1 = 0.57282 D2 = -1.76919389131923 D3 = 0.759263343792048 D4 = -0.104902108710045 Return 4000 'Subroutine DECB(N, E, MLDE, MUE, ITMP, INFO) ' Labels 4007, 4020, 4030, 4045, 4060, 4070, 4080 '----------------------------------------------------------------------- ' MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION OF A BANDED ' MATRIX WITH LOWER BANDWIDTH MLDE AND UPPER BANDWIDTH MUE ' INPUTS: ' N ORDER OF THE ORIGINAL MATRIX A. ' E CONTAINS THE MATRIX IN BAND STORAGE. THE COLUMNS ' OF THE MATRIX ARE STORED IN THE COLUMNS OF E AND ' THE DIAGONALS OF THE MATRIX ARE STORED IN ROWS ' MLDE+1 THROUGH 2*MLDE+MUE+1 OF E. ' MLDE LOWER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED). ' MUE UPPER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED). ' OUTPUTS: ' E AN UPPER TRIANGULAR MATRIX IN BAND STORAGE AND ' THE MULTIPLIERS WHICH WERE USED TO OBTAIN IT. ' ITMP INDEX VECTOR OF PIVOT INDICES. ' ITMP(N) (-1)^(NUMBER OF INTERCHANGES) OR O. ' INFO = 0 IF MATRIX A IS NONSINGULAR, OR = K IF FOUND TO BE ' SINGULAR AT STAGE K. ' NOTE: USE SOLB TO OBTAIN SOLUTION OF LINEAR SYSTEM. ' DETERM(E) = ITMP(N)*E(MD,1)*E(MD,2)*...*E(MD,N) WITH MD=MLDE+MUE+1. ' IF ITMP(N)=O, E IS SINGULAR, SOLB WILL DIVIDE BY ZERO. ' ' REFERENCE.. ' THIS IS A MODIFICATION OF ' C. B. MOLDER, ALGORITHM 423, LINEAR EQUATION SOLVER, ' C.A.C.M. 15 (1972), P. 274. '----------------------------------------------------------------------- INFO = 0 ITMP(N) = 1 MD = MLDE + MUE + 1 MD1 = MD + 1 JU = 0 If MLDE = 0 Then GoTo 4070 If N = 1 Then GoTo 4070 If N < MUE + 2 Then GoTo 4007 For J = MUE + 2 To N For I = 1 To MLDE E(I, J) = ZERO Next I Next J 4007 NM1 = N - 1 For K = 1 To NM1 KP1 = K + 1 M = MD 'MDL = IMIN(MLDE,N-K) + MD ia = MLDE: ib = N - K: GoSub 610 MDL = IMIN + MD For I = MD1 To MDL If Abs(E(I, K)) > Abs(E(M, K)) Then M = I Next I ITMP(K) = M + K - MD t = E(M, K) If M = MD Then GoTo 4020 ITMP(N) = -ITMP(N) E(M, K) = E(MD, K) E(MD, K) = t 4020 If t = ZERO Then GoTo 4080 t = ONE / t For I = MD1 To MDL E(I, K) = -E(I, K) * t Next I 'JU = IMIN(IMAX(JU,MUE+ITMP(K)),N) ia = JU: ib = MUE + ITMP(K): GoSub 600 ia = IMAX: ib = N: GoSub 610: JU = IMIN MM = MD If JU < KP1 Then GoTo 4060 For J = KP1 To JU M = M - 1 MM = MM - 1 t = E(M, J) If M = MM Then GoTo 4030 E(M, J) = E(MM, J) E(MM, J) = t 4030 If t = ZERO Then GoTo 4045 JK = J - K For I = MD1 To MDL IJK = I - JK E(IJK, J) = E(IJK, J) + E(I, K) * t Next I 4045 Next J 4060 Next K 4070 K = N If E(MD, N) = ZERO Then GoTo 4080 Return 4080 INFO = K ITMP(N) = 0 Return 'DECB 4500 'Subroutine DECA(N, E, ITMP, INFO) ' REAL DOUBLE PRECISION VERSION OF DEC ' Labels: 4520, 4550, 4570, 4580 '----------------------------------------------------------------------- ' GENERAL MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION. ' INPUT.. ' N = ORDER OF MATRIX. ' E = MATRIX TO BE TRIANGULARIZED. ' OUTPUT.. ' E(I,J), I <= J = UPPER TRIANGULAR FACTOR, U. ' E(I,J), I > J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I - L. ' ITMP(K), K < N = INDEX OF K-TH PIVOT ROW. ' ITMP(N) = (-1)^(NUMBER OF INTERCHANGES) OR O. ' INFO = 0 IF MATRIX E IS NONSINGULAR, OR K IF FOUND TO BE ' SINGULAR AT STAGE K. ' NOTE: USE SOL TO OBTAIN SOLUTION OF LINEAR SYSTEM. ' DETERM(E) = ITMP(N)*E(1,1)*E(2,2)*...*E(N,N). ' IF ITMP(N)=O, A IS SINGULAR, SOL WILL DIVIDE BY ZERO. ' ' REFERENCE.. ' C. B. MOLDER, ALGORITHM 423, LINEAR EQUATION SOLVER, ' C.A.C.M. 15 (1972), P. 274. '------------------------------------------------------------------------ INFO = 0 ITMP(N) = 1 If N = 1 Then GoTo 4570 NM1 = N - 1 For K = 1 To NM1 KP1 = K + 1 M = K For I = KP1 To N If Abs(E(I, K)) > Abs(E(M, K)) Then M = I Next I ITMP(K) = M t = E(M, K) If M = K Then GoTo 4520 ITMP(N) = -ITMP(N) E(M, K) = E(K, K) E(K, K) = t 4520 If t = ZERO Then GoTo 4580 t = ONE / t For I = KP1 To N E(I, K) = -E(I, K) * t Next I For J = KP1 To N t = E(M, J) E(M, J) = E(K, J) E(K, J) = t If t = ZERO Then GoTo 4550 For I = KP1 To N E(I, J) = E(I, J) + E(I, K) * t Next I 4550 Next J Next K 4570 K = N If E(N, N) = ZERO Then GoTo 4580 Return 4580 INFO = K ITMP(N) = 0 Return 5000 'Subroutine SOLB(N, E, MLDE, MUE, B, ITMP) ' Labels 5025, 5050 '----------------------------------------------------------------------- ' SOLUTION OF A BANDED LINEAR SYSTEM, E*X = B. ' INPUTS: ' N ORDER OF MATRIX A. ' E TRIANGULARIZED MATRIX OBTAINED FROM DECB. ' MLDE LOWER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED). ' MUE UPPER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED). ' B RIGHT HAND SIDE VECTOR. ' ITMP PIVOT VECTOR OBTAINED FROM DECB. ' DO NOT USE IF DECB HAS SET INFO <> 0. ' OUTPUT: ' B SOLUTION VECTOR. '----------------------------------------------------------------------- MD = MLDE + MUE + 1 MD1 = MD + 1 MDM = MD - 1 NM1 = N - 1 If MLDE = 0 Then GoTo 5025 If N = 1 Then GoTo 5050 For K = 1 To NM1 M = ITMP(K) t = B(M) B(M) = B(K) B(K) = t 'MDL = IMIN(MLDE,N-K) + MD ia = MLDE: ib = N - K: GoSub 610 MDL = IMIN + MD For I = MD1 To MDL IMD = I + K - MD B(IMD) = B(IMD) + E(I, K) * t Next I Next K 5025 For KB = 1 To NM1 K = N + 1 - KB B(K) = B(K) / E(MD, K) t = -B(K) KMD = MD - K 'LM = IMAX(1,KMD+1) ia = 1: ib = KMD + 1: GoSub 600 LM = IMAX For I = LM To MDM IMD = I - KMD B(IMD) = B(IMD) + E(I, K) * t Next I Next KB 5050 B(1) = B(1) / E(MD, 1) Return 5500 'Subroutine SOL (N, E, B, ITMP) ' Label: 5550 '----------------------------------------------------------------------- ' SOLUTION OF A GENERAL LINEAR SYSTEM, E*X = B. ' INPUTS: ' N = ORDER OF MATRIX. ' E = TRIANGULARIZED MATRIX OBTAINED FROM DECA. ' B = RIGHT HAND SIDE VECTOR. ' ITMP = PIVOT VECTOR OBTAINED FROM DEC. ' DO NOT USE IF DEC HAS SET INFO <> 0. ' OUTPUT: ' B = SOLUTION VECTOR. '----------------------------------------------------------------------- If N = 1 Then GoTo 5550 NM1 = N - 1 For K = 1 To NM1 KP1 = K + 1 M = ITMP(K) t = B(M) B(M) = B(K) B(K) = t For I = KP1 To N B(I) = B(I) + E(I, K) * t Next I Next K For KB = 1 To NM1 KM1 = N - KB K = KM1 + 1 B(K) = B(K) / E(K, K) t = -B(K) For I = 1 To KM1 B(I) = B(I) + E(I, K) * t Next I Next KB 5550 B(1) = B(1) / E(1, 1) Return 'end of file tros4.bas
Saturday, November 18, 2017
NUMERICAL SOLUTION OF A STIFF SYSTEM OF FIRST ORDER ORDINARY DIFFERENTIAL EQUATIONS Y':=F(X,Y) BY ROSENBROCK METHOD.
Subscribe to:
Posts (Atom)