DECLARE SUB extrapol (row AS INTEGER, fhilf#(), y#(), n%, epsabs#, epsrel#, x#, h#, h0#, ahead AS INTEGER, index%, extmax AS INTEGER, bufol() AS INTEGER) DECLARE FUNCTION ibulstoe% (x#, xend#, n%, y#(), epsabs#, epsrel#, h#, hmax#, neu%, fmax AS INTEGER, aufrufe AS INTEGER) DECLARE SUB abmschritt (x0#, n%, h#, methode%, aufrufe AS INTEGER, f#(), tmp#(), kisum#(), y1#(), y2#(), diff#()) DECLARE SUB printvec (Name1 AS STRING, n%, v#()) DECLARE FUNCTION Max% (ia%, ib%) DECLARE FUNCTION Min% (ia%, ib%) DECLARE SUB copyvector (A#(), b#(), n%) DECLARE SUB init0vector (vektor#(), n%) DECLARE SUB incvector (ziel#(), quelle#(), factor#, n%) DECLARE SUB addvector (summe#(), summand1#(), summand2#(), factor#, n%) DECLARE SUB rkschritt (x0#, y0#(), n%, f#(), tmp#(), kisum#(), y1#(), y2#(), diff#(), h#) DECLARE FUNCTION irkstart% (x#, x0#, y#(), n%, xend#, h#, hmax#, newstep%, methode%, aufrufe AS INTEGER, f#(), tmp#(), kisum#(), y1#(), y2#(), diff#()) DECLARE FUNCTION ipraekorr% (x#, y#(), n%, xend#, h#, epsabs#, epsrel#, fmax AS INTEGER, aufrufe AS INTEGER, hmax#, neu%) DECLARE SUB gauss (mode%, n%, mat() AS DOUBLE, lumat() AS DOUBLE, perm() AS INTEGER, b#(), x#(), signd AS INTEGER, code AS INTEGER) DECLARE SUB gaudec (n%, mat() AS DOUBLE, lumat() AS DOUBLE, perm() AS INTEGER, signd AS INTEGER, rc AS INTEGER) DECLARE FUNCTION det# (n%, mat() AS DOUBLE) DECLARE SUB dgl (num%, n%, x#, y#(), f#()) DECLARE FUNCTION irwp% (A#, b#, h#, ystart#(), n%, awpnumm AS INTEGER, epsawp#, epsrb#, fmax AS INTEGER, itmax%, actiter AS INTEGER) DECLARE SUB gausol (n%, lumat() AS DOUBLE, perm() AS INTEGER, b#(), x#(), rc AS INTEGER) DECLARE SUB awp (x#, xend#, bspn AS INTEGER, n%, y#(), epsabs#, epsrel#, h#, methode%, fmax AS INTEGER, aufrufe AS INTEGER, fehler AS INTEGER) DECLARE SUB ruku23 (x#, y#(), bspn AS INTEGER, n%, h#, y2#(), y3#()) DECLARE SUB engl45 (x#, y#(), bspn AS INTEGER, n%, h#, y4#(), y5#()) DECLARE SUB prdo45 (x#, y#(), bspn AS INTEGER, n%, h#, y4#(), y5#(), steif1 AS INTEGER, steifanz AS INTEGER, steif2 AS INTEGER) DECLARE FUNCTION DistMax# (vector1#(), vector2#(), n%) DECLARE FUNCTION XMax# (A#, b#) DECLARE FUNCTION XMin# (A#, b#) DECLARE FUNCTION XNormMax# (vektor#(), n%) DECLARE SUB zeig (s AS STRING, v#(), n%) DECLARE SUB randbed (ya#(), yb#(), R#()) DECLARE SUB ISWAP (ia%, ib%) DECLARE SUB XSWAP (A#, b#) '************************************************************************ '* Solve a boundary value problem for a first order DE system of the * '* form: * '* * '* Y1' = F1(X,Y1,Y2,...,YN) * '* Y2' = F2(X,Y1,Y2,...,YN) * '* ... * '* YN' = FN(X,Y1,Y2,...,YN) * '* * '* via the shooting method by determining an approximation for the * '* initial value Y(A). * '* -------------------------------------------------------------------- * '* SAMPLE RUNS; * '* * '* TEST EXAMPLE (WITH METHOD #1): * '* Method: Runge-Kutta embedding formula of 4/5th order. * '* ============= * '* ANALYZED SET OF DIFFERENTIAL EQUATIONS: * '* --------------------------------------- * '* Y'( 1 ) = Y(2) * '* Y'( 2 ) = -Y(1)**3 * '* WITH BOUNDARY CONDITIONS : VALUE OF Y(1) * '* AT POINT A = 0 : 0 * '* AT POINT B = 1 : 0 * '* * '* AT LOCATION A = 0 THE FOLLOWING VALUES ARE PROVIDED: * '* * '* Y( 1 ) = 0 * '* Y( 2 ) = 12 * '* * '* REQUIRED PARAMETER: * '* ------------------: * '* -PRECISION FOR THE IVP = 1.0000000000E-07 * '* -PRECISION FOR THE BVP = 1.0000000000E-08 * '* -STEP WIDTH FOR THE IVP= 1.0000000000E-02 * '* -MAX. NUMBER OF F-EVALUATIONS FOR SOLUTION OF THE IVP'S = 20000 * '* -MAX. NUMBER OF NEWTON-ITERATION STEPS = 1000 * '* * '* SOLUTION OF THE PROBLEM: * '* ------------------------ * '* THE START VALUE OF A SOLUTION OF THE BVP IS APPROXIMATED * '* IN POINT A = 0.000E+00 AS FOLLOWS: * '* * '* Y( 1 ) = 0 * '* Y( 2 ) = 9.722981024221486 * '* * '* THIS REQUIRES 4 NEWTON-ITERATIONS! * '* * '* DETERMINATION COMPLETED AS PLANNED! * '* * '* * '* TEST EXAMPLE (WITH METHOD #2): * '* Method: Predictor-corrector method of order 4 by Adams-Bashforth- * '* Moulton. * '* ============= * '* ANALYZED SET OF DIFFERENTIAL EQUATIONS: * '* --------------------------------------- * '* Y'( 1 ) = Y(2) * '* Y'( 2 ) = -Y(1)**3 * '* WITH BOUNDARY CONDITIONS : VALUE OF Y(1) * '* AT POINT A = 0 : 0 * '* AT POINT B = 1 : 0 * '* * '* AT LOCATION A = 0 THE FOLLOWING VALUES ARE PROVIDED: * '* * '* Y( 1 ) = 0 * '* Y( 2 ) = 12 * '* * '* REQUIRED PARAMETER: * '* ------------------: * '* -PRECISION FOR THE IVP = .0000001 * '* -PRECISION FOR THE BVP = .00000001 * '* -STEP WIDTH FOR THE IVP= .01 * '* -MAX. NUMBER OF F-EVALUATIONS FOR SOLUTION OF THE IVP'S = 20000 * '* -MAX. NUMBER OF NEWTON-ITERATION STEPS = 1000 * '* * '* SOLUTION OF THE PROBLEM: * '* ------------------------ * '* THE START VALUE OF A SOLUTION OF THE BVP IS APPROXIMATED * '* IN POINT A = 0.000E+00 AS FOLLOWS: * '* * '* Y( 1 ) = 0 * '* Y( 2 ) = 9.722980952342329 * '* * '* THIS REQUIRES 4 NEWTON-ITERATIONS! * '* * '* DETERMINATION COMPLETED AS PLANNED! * '* * '* * '* TEST EXAMPLE (WITH METHOD #3): * '* Method: Extrapolation method of Bulirsch-Stoer. * '* ============= * '* ANALYZED SET OF DIFFERENTIAL EQUATIONS: * '* --------------------------------------- * '* Y'( 1 ) = Y(2) * '* Y'( 2 ) = -Y(1)**3 * '* WITH BOUNDARY CONDITIONS : VALUE OF Y(1) * '* AT POINT A = 0 : 0 * '* AT POINT B = 1 : 0 * '* * '* AT LOCATION A = 0 THE FOLLOWING VALUES ARE PROVIDED: * '* * '* Y( 1 ) = 0 * '* Y( 2 ) = 12 * '* * '* REQUIRED PARAMETER: * '* ------------------: * '* -PRECISION FOR THE IVP = .0000001 * '* -PRECISION FOR THE BVP = .00000001 * '* -STEP WIDTH FOR THE IVP= .01 * '* -MAX. NUMBER OF F-EVALUATIONS FOR SOLUTION OF THE IVP'S = 20000 * '* -MAX. NUMBER OF NEWTON-ITERATION STEPS = 1000 * '* * '* SOLUTION OF THE PROBLEM: * '* ------------------------ * '* THE START VALUE OF A SOLUTION OF THE BVP IS APPROXIMATED * '* IN POINT A = 0 AS FOLLOWS: * '* * '* Y( 1 ) = 0 * '* Y( 2 ) = 9.722981027825767 * '* * '* THIS REQUIRES 4 NEWTON-ITERATIONS! * '* * '* DETERMINATION COMPLETED AS PLANNED! * '* * '* -------------------------------------------------------------------- * '* REF.: "Numerical Algorithms with C, By Gisela Engeln-Muellges * '* and Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11]. * '* * '* Quick Basic Release 1.1 By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '* * '* Release 1.1 (12/25/06): corrected bug in ipraekorr (added xsave). * '************************************************************************ ' Initialize data; test example can be changed by adjusting the ' subroutines dgl, randbed, and the variable n in main program. 'Program Test_rwp DEFDBL A-H, O-Z DEFINT I-N OPTION BASE 0 DIM TMethod(3) AS STRING DIM fehler AS INTEGER TMethod(0) = "Runge-Kutta embedding formula of 4/5th order." TMethod(1) = "Predictor-corrector method of order 4 by Adams-Bashforth-Moulton." TMethod(2) = "Extrapolation method of Bulirsch-Stoer." n = 2 'Number of equations DIM ystri(n) AS STRING DIM yanf(n), yanf0(n) yanf0(0) = 0# yanf0(1) = 12# ystri(0) = "Y(2)" 'alpha-numeric form ystri(1) = "-Y(1)^3" 'of the DE A = 0# b = 1# hstep = .01# wa = 0# wb = 0# ifmax = 20000 itmax = 1000 epsawp = .0000001# epsrb = .00000001# 'output of the test examples for three varying methods PRINT FOR method = 0 TO 2 'yanf = yanf0 FOR i = 0 TO n - 1 yanf(i) = yanf0(i) NEXT i CLS PRINT "* TEST EXAMPLE (WITH METHOD #"; method + 1; ")" PRINT "* Method: "; TMethod(method) PRINT "=============" PRINT "ANALYZED SET OF DIFFERENTIAL EQUATIONS:" PRINT "---------------------------------------" FOR j = 0 TO n - 1 PRINT "* Y'("; j + 1; ") = "; ystri(j) NEXT j PRINT "WITH BOUNDARY CONDITIONS : VALUE OF Y(1)" PRINT "* AT POINT A = "; A; ": "; wa PRINT "* AT POINT B = "; b; ": "; wb PRINT "* AT LOCATION A = "; A; ", THE FOLLOWING VALUES ARE PROVIDED:" FOR j = 0 TO n - 1 PRINT "* Y("; j + 1; ") = "; yanf(j) NEXT j PRINT "REQUIRED PARAMETERS:" PRINT "------------------- " PRINT "* -PRECISION FOR THE IVP = "; epsawp PRINT "* -PRECISION FOR THE BVP = "; epsrb PRINT "* -STEP WIDTH FOR THE IVP= "; hstep PRINT "* -MAX. NUMBER OF F-EVALUATIONS FOR SOLUTION OF THE IVP'S = "; ifmax PRINT "* -MAX. NUMBER OF NEWTON-ITERATION STEPS = "; itmax 'call function irwp fehler = irwp(A, b, hstep, yanf(), n, method, epsawp, epsrb, ifmax, itmax, iter) 'print results IF fehler <> 0 THEN PRINT "* "; IF fehler = 1 THEN PRINT "ERRORBOUND(S) TOO SMALL" IF fehler = 2 THEN PRINT "ERROR: B <= A" IF fehler = 3 THEN PRINT "ERROR: STEP SIZE H <= 0" IF fehler = 4 THEN PRINT "ERROR: NUMBER N OF DEQ'S NOT CORRECT" IF fehler = 5 THEN PRINT "ERROR: WRONG IVP-PROGRAM CHOICE" IF fehler = 6 THEN PRINT "IFMAX TO SMALL FOR SOLUTION OF THE IVP-PROGRAM" IF fehler = 7 THEN PRINT "AFTER ITMAX STEPS PRECISION WAS NOT REACHED!" IF fehler = 8 THEN PRINT "ERROR: JACOBI-MATRIX IS SINGULAR!" ELSE 'no error? 'Output of solution PRINT "*" PRINT "SOLUTION OF THE PROBLEM:" PRINT "----------------------- " PRINT "THE START VALUE OF A SOLUTION OF THE BVP IS APPROXIMATED" PRINT "* IN POINT A = "; A; " AS FOLLOWS:" FOR j = 0 TO n - 1 PRINT "* Y("; j + 1; ") = "; yanf(j) NEXT j PRINT "*" PRINT "* THIS REQUIRES "; iter; " NEWTON-ITERATIONS!" PRINT "DETERMINATION COMPLETED AS PLANNED! " END IF INPUT R$ NEXT method END 'of main program SUB abmschritt (x0, n, h, methode, aufrufe AS INTEGER, f(), tmp(), kisum() AS DOUBLE, y1(), y2(), diff()) '************************************************************************ '* Perform one step of the Adams-Bashforth-Moulton method * '* * '* Input parameters: * '* ================= * '* x0 initial x-value * '* n number of DEs in system * '* h step size * '* hilf Structure with the vectors f[1]..f[4] for the starting values * '* * '* Output parameter: * '* ================= * '* x0 final x-value of the integration * '* methode equal to 1; hence the error estimation in ipraekorr is * '* performed with the factor for Adams-Bashforth-Moulton. * '* aufrufe current number of calls of dgl * '* hilf Structure with the following output : * '* y1: approximate solution used for the error estimation * '* (from predictor step) * '* y2: approximate solution from corrector step * '* f(4): new node for starting values * '* f(0),...,f(3) are altered and tmp is used as aux. storage * '* * '* global names used: * '* ================== * '* Init0vector, incvector, addvector. * '************************************************************************ DIM prae(4), tmp2(n) DIM korr(4) AS DOUBLE ' ---- the coefficients for Adams-Bashforth-Moulton ---------------- prae(0) = -9#: prae(1) = 37# prae(2) = -59#: prae(3) = 55# 'Predictor korr(0) = 1#: korr(1) = -5# korr(2) = 19#: korr(3) = 9# 'corrector tmp1 = f(0) 'the starting values are expected FOR j = 0 TO 3 'reside in f[0], ..., f[3] , but f(j) = f(j + 1) 'are stored in f[1]..f[4], we NEXT j 'rotate their pointer. f(4) = tmp1 init0vector y1(), n 'one predictor step FOR j = 0 TO 3 FOR i = 0 TO n - 1 tmp2(i) = f(j + i) NEXT i incvector y1(), tmp2(), prae(j), n NEXT j addvector y1(), y2(), y1(), h / 24#, n x0 = x0 + h 'move on by h; compute rhs dgl 0, n, x0, y1(), tmp2() 'of DE at (x0,y1) : this 'yields the start for the 'corrector. FOR i = 0 TO n - 1 f(4 + i) = tmp2(i) NEXT i init0vector tmp(), n 'one corrector step FOR j = 0 TO 3 FOR i = 0 TO n - 1 tmp2(i) = f(j + 1 + i) NEXT i incvector tmp(), tmp2(), korr(j), n NEXT j incvector y2(), tmp(), h / 24#, n 'y2: new approximate 'solution. dgl 0, n, x0, y2(), tmp2() 'f(4) comtains the node for 'the starting values FOR i = 0 TO n - 1 'derived from the corrector f(4 + i) = tmp2(i) NEXT i methode = 1 aufrufe = aufrufe + 2 '2 calls of dgl for one A-B-M step END SUB SUB addvector (summe(), summand1(), summand2(), factor, n) '************************************************************************ '* add factor * summand2 to the (0:n-1) vector summand1, store result in* '* vector summe. * '* * '* global names used: * '* ================== * '* None * '************************************************************************ FOR i = 0 TO n - 1 summe(i) = summand1(i) + factor * summand2(i) NEXT i END SUB '************************************************************************ '* * '* Solve an ordinary system of first order differential equations using * '* -------------------------------------------------------------------- * '* automatic step size control * '* ---------------------------- * '* * '* Programming language: ANSI C * '* Author: Klaus Niederdrenk (FORTRAN) * '* Adaptation: Juergen Dietel, Computer Center, RWTH Aachen * '* Source: existing C, Pascal, QuickBASIC and FORTRAN * '* codes * '* Date: 6.2.1992, 10.2.1995 * '* * '* Quick Basic Release By J-P Moreau, Paris. * '* -------------------------------------------------------------------- * '* REF.: "Numerical Algorithms with C, by Gisela Engeln-Muellges * '* and Frank Uhlig, Springer-Verlag, 1996". * '************************************************************************ ' 1st order DESs with automatic step size control ........................... SUB awp (x, xend, bspn AS INTEGER, n, y(), epsabs, epsrel, h, methode, fmax AS INTEGER, aufrufe AS INTEGER, fehler AS INTEGER) ' double x, initial/final x value .............. ' xend desired end point .................. ' integer bspn, # example ' n number of DEs ...................... ' double y(0:n-1), initial/final y value .............. ' epsabs, absolute error bound ............... ' epsrel, relative error bound ............... ' h initial/final step size ............ ' integer methode, desired method (3, 6, 7) ........... ' fmax, maximal # of calls of dgl ......... ' aufrufe, actual # of calls of dgl .......... ' fehler error code ......................... '************************************************************************ '* Compute the solution y of a system of first order ordinary * '* differential equations y' = f(x,y) at xend from the given * '* initial data (x0, y0). * '* We use automatic step size control internally so that the error of * '* y (absolutely or relatively) lies within the given error bounds * '* epsabs and epsrel. * '* * '* Input parameters: * '* ================= * '* x initial value for x * '* y initial values for y(0:n-1) * '* bspn # example * '* n number of differential equations * '* xend end of integration; xend > x0 * '* h initial step size * '* epsabs absolute error bound; >= 0; if = 0 we only check the * '* relative error. * '* epsrel relative error bound; >= 0; if = 0 we check only the * '* absolute eror. * '* fmax max number of evaluations of right hand side in dgl * '* methode chooses the method * '* = 3: Runge-Kutta method of 2nd/3rd order * '* = 6: England formula of 4th/5th order * '* = 7: Formula of Prince-Dormand of 4th/5th order * '* * '* Output parameters: * '* ================== * '* x final x-value for iteration. If fehler = 0 we usually have * '* x = xend. * '* y final y-values for the solution at x * '* h final step size used; leave for subsequent calls * '* aufrufe actual number of calls of dgl * '* * '* Return value (fehler): * '* ===================== * '* = 0: all ok * '* = 1: both error bounds chosen too small for the given mach. constant * '* = 2: xend <= x0 * '* = 3: h <= 0 * '* = 4: n <= 0 * '* = 5: more right hand side calls than allowed: aufrufe > fmax, * '* x and h contain the current values when stop occured. * '* = 6: improper input for embedding formula * '* = 7: lack of available memory (not used here) * '* = 8: Computations completed, but the Prince Dormand formula stiff- * '* ness test indicates possible stiffness. * '* = 9: Computations completed, but both Prince Dormand formula stiff- * '* ness tests indicate possible stiffness. Use method for stiff * '* systems instead ' * '* =10: aufrufe > fmax, see error code 5; AND the Prince Dormand formula* '* indicates stiffness; retry using a stiff DE solver ' * '* * '************************************************************************ DIM MACHEPS AS DOUBLE DIM mach2 AS DOUBLE DIM mach AS DOUBLE 'machine constant dependent variable which 'avoids using too little steps near xend. DIM xendh AS DOUBLE '|xend| - MACH2, carrying same sign as xend DIM ymax AS DOUBLE 'Maximum norm of newest approximation of max 'order. DIM hhilf AS DOUBLE 'aux storage for the latest value of h 'produced by step size control. It is saved 'here in order to avoid to return a `h' that 'resulted from an arbitrary reduction at the 'end of the interval. DIM diff AS DOUBLE 'distance of the two approximations from the 'embedding formula. DIM s AS DOUBLE 'indicates acceptance level for results from 'embeding formula. 'approximate solution of low order DIM ybad(n) 'ditto of high order DIM ygood(n) DIM amEnde AS INTEGER 'flag that shows if the end of the interval 'can be reached with the actual step size. DIM fertig AS INTEGER 'flag indicating end of iterations. DIM steif1 AS INTEGER 'Flag, that is set in prdo45() if its 'stiffness test (dominant eigenvalue) 'indicates so. Otherwise no changes. DIM steifanz AS INTEGER 'counter for number of successive successes 'of stiffness test of Shampine and Hiebert in 'prdo45(). DIM steif2 AS INTEGER 'Flag, set in prdo45(), when the stiffness 'test of Shampine and Hiebert wa successful 'three times in a row; otherwise no changes. 'initialize some variables fehler = 0 MACHEPS = 2.25D-16 mach2 = 100 * MACHEPS mach = MACHEPS ^ .75 amEnde = 0 fertig = 0 steif1 = 0 steif2 = 0 steifanz = 0 aufrufe = 1 ymax = XNormMax(y(), n) IF xend >= 0# THEN xendh = xend * (1# - mach2) ELSE xendh = xend * (1# + mach2) END IF ' ----------------------- check inputs ---------------------- IF epsabs <= mach2 * ymax AND epsrel <= mach2 THEN fehler = 1 EXIT SUB END IF IF xendh < x THEN fehler = 2 EXIT SUB END IF IF h < mach2 * ABS(x) THEN fehler = 3 EXIT SUB END IF IF n <= 0 THEN fehler = 4 EXIT SUB END IF IF methode <> 3 AND methode <> 6 AND methode <> 7 THEN fehler = 6 EXIT SUB END IF ' ********************************************************************** ' * * ' * I t e r a t i o n s * ' * * ' ********************************************************************** IF x + h > xendh THEN 'almost at end point ? hhilf = h 'A shortened step might be h = xend - x 'enough. amEnde = 1 END IF WHILE fertig = 0 'solve DE system by integrating from 'x0 to xend by suitable steps. 'choose method IF methode = 3 THEN CALL ruku23(x, y(), bspn, n, h, ybad(), ygood()) ELSEIF methode = 6 THEN CALL engl45(x, y(), bspn, n, h, ybad(), ygood()) ELSEIF methode = 7 THEN CALL prdo45(x, y(), bspn, n, h, ybad(), ygood(), steif1, steifanz, steif2) END IF aufrufe = aufrufe + methode diff = DistMax(ybad(), ygood(), n) IF (diff < mach2) THEN 'compute s s = 2# ELSE ymax = XNormMax(ygood(), n) s = SQR(h * (epsabs + epsrel * ymax) / diff) IF methode <> 3 THEN s = SQR(s) END IF IF s > 1# THEN 'integration acceptable ? FOR i = 0 TO n - 1 'accept highest order solution y(i) = ygood(i) 'move x NEXT i x = x + h IF amEnde <> 0 THEN 'at end of interval ? fertig = 1 'stop iteration IF methode = 7 THEN IF steif1 > 0 OR steif2 > 0 THEN fehler = 8 IF steif1 > 0 AND steif2 > 0 THEN fehler = 9 END IF ELSEIF aufrufe > fmax THEN 'too many calls of dgl ? hhilf = h 'save actual step size fehler = 5 'report error and stop fertig = 1 IF methode = 7 AND (steif1 > 0 OR steif2 > 0) THEN fehler = 10 ELSE 'Integration was successful 'not at the interval end ? h = h * XMin(2#, .98 * s) 'increase step size for next 'step properly, at most by 'factor two. Value `0.98*s' is 'used in order to avoid that 'the theoretical value s is 'exceeded by accidental 'rounding errors. IF x + h > xendh THEN 'nearly reached xend ? hhilf = h '=> One further step with h = xend - x 'reduced step size might be amEnde = 1 'enough. IF h < mach * ABS(xend) THEN 'very close to xend ? fertig = 1 'finish iteration. END IF END IF END IF ELSE 'step unsuccessful ? 'before repeating this step h = h * XMax(.5, .98 * s) 'reduce step size properly, at 'most by factor 1/2 (for factor amEnde = 0 '0.98: see above). END IF WEND h = hhilf 'return the latest step size computed by step 'size control and error code to the caller. END SUB SUB copyvector (A(), b(), n) FOR i = 0 TO n - 1 A(i) = b(i) NEXT i END SUB 'Determinant .............................. FUNCTION det (n, mat() AS DOUBLE) ' n, 'Dimension of the matrix .......... ' mat() as Double '(0:n-1,0:n-1) matrix here a vector '*====================================================================* '* * '* det computes the determinant of an n x n real matrix mat * '* * '*====================================================================* '* * '* Input parameter: * '* ================ * '* n integer; (n > 0) * '* Dimension of mat * '* mat (n x n matrix (here a n*n vector). * '* * '* Return value: * '* ============= * '* Double Determinant of mat. * '* If the return value = 0, then the matrix is singular * '* or the storage is insufficient * '* * '*====================================================================* '* * '* Functions in use : * '* =================== * '* * '* int gaudec (): LU decomposition of mat * '* * '*====================================================================* DIM MACHEPS AS DOUBLE DIM MAXROOT AS DOUBLE DIM rc AS INTEGER DIM signd AS INTEGER DIM perm(n) AS INTEGER DIM lu(n * n) AS DOUBLE MACHEPS = 2.25D-16 MAXROOT = 1D+16 IF n < 1 THEN det = 0# EXIT FUNCTION END IF gaudec n, mat(), lu(), perm(), signd, rc 'decompose IF rc <> 0 OR signd = 0 THEN det = 0# EXIT FUNCTION END IF tmpdet = 1# * signd FOR i = 0 TO n - 1 IF ABS(tmpdet) < MACHEPS THEN det = 0# EXIT FUNCTION ELSEIF ABS(tmpdet) > MAXROOT OR ABS(lu(i * n + i)) > MAXROOT THEN det = MAXROOT EXIT FUNCTION ELSE tmpdet = tmpdet * lu(i * n + i) 'compute det END IF NEXT i det = tmpdet END FUNCTION 'System of DEs to integrate '========================== SUB dgl (num, n, x, y(), f()) 'In this example, num,n,x are not used. f(0) = y(1) f(1) = -y(0) ^ 3 END SUB 'Maximum norm of a difference vector ........ FUNCTION DistMax (vector1(), vector2(), n) '************************************************************************ '* Compute the maximum norm of the difference of two [0..n-1] vectors * '* * '* global Name1 used: * '* ================= * '* None * '************************************************************************ ' double abstand reference value for computation of distance ' double hilf distance of two vector elements abstand = 0# FOR i = n - 1 TO 0 STEP -1 hilf = ABS(vector1(i) - vector2(i)) IF hilf > abstand THEN abstand = hilf NEXT i DistMax = abstand END FUNCTION ' England's Einbettungs formulas of 4th and 5th degree ................ SUB engl45 (x, y(), bspn AS INTEGER, n, h, y4(), y5()) ' double x starting point of integration ........ ' double y(0:n-1) initial value at x ................... ' integer bspn, # example ' n number of differential equations ..... ' double h step size ............................ ' double y4(0:n-1), 4th order approximation for y at x + h ' y5(0:n-1) 5th order approximation for y at x + h ' auxiliary vectors DIM yhilf(n) DIM k1(n) AS DOUBLE DIM k2(n) AS DOUBLE DIM k3(n) AS DOUBLE DIM k4(n) AS DOUBLE DIM k5(n) AS DOUBLE DIM k6(n) AS DOUBLE '************************************************************************ '* Compute 4th and 5th order approximates y4, y5 at x + h starting with * '* a solution y at x by using the England embedding formulas on the * '* first order system of n differential equations y' = f(x,y) , as * '* supplied by dgl(). * '* * '* Input parameters: * '* ================= * '* x initial x-value * '* y y-values at x, type pVEC * '* n number of differential equations * '* dgl function that evaluates the right hand side of the system * '* y' = f(x,y) * '* h step size * '* * '* yhilf, k1..K6: auxiliary vectors * '* * '* Output parameters: * '* ================== * '* y4 4th order approximation for y at x + h (pVEC) * '* y5 5th order approximation for y at x + h (pVEC) * '* * '************************************************************************ CALL dgl(bspn, n, x, y(), k1()) FOR i = 0 TO n - 1 yhilf(i) = y(i) + .5 * h * k1(i) NEXT i CALL dgl(bspn, n, x + .5 * h, yhilf(), k2()) FOR i = 0 TO n - 1 yhilf(i) = y(i) + (.25 * h * (k1(i) + k2(i))) NEXT i CALL dgl(bspn, n, x + .5 * h, yhilf(), k3()) FOR i = 0 TO n - 1 yhilf(i) = y(i) + h * (-k2(i) + 2# * k3(i)) NEXT i CALL dgl(bspn, n, x + h, yhilf(), k4()) FOR i = 0 TO n - 1 yhilf(i) = y(i) + h / 27# * (7# * k1(i) + 10# * k2(i) + k4(i)) NEXT i CALL dgl(bspn, n, x + (2# / 3#) * h, yhilf(), k5()) FOR i = 0 TO n - 1 yhilf(i) = y(i) + h / 625# * (28# * k1(i) - 125# * k2(i) + 546# * k3(i) + 54# * k4(i) - 378# * k5(i)) NEXT i CALL dgl(bspn, n, x + h / 5#, yhilf(), k6()) FOR i = 0 TO n - 1 y4(i) = y(i) + h / 6# * (k1(i) + 4# * k3(i) + k4(i)) y5(i) = y(i) + h / 336# * (14# * k1(i) + 35# * k4(i) + 162# * k5(i) + 125# * k6(i)) NEXT i END SUB SUB extrapol (row AS INTEGER, fhilf(), y(), n, epsabs, epsrel, x, h, h0, ahead AS INTEGER, index, extmax AS INTEGER, bufol() AS INTEGER) ' row: integer ' fhilf(12,2) ' y(n) ' ahead: boolean ' bufol(12): integer (bulirsch sequence) '************************************************************************ '* Perform one extrapolation step for function ibulstoe * '* * '* Input parameters: * '* ================= * '* row pointer to the arrays bufol and fhilf * '* fhilf pointer to Matrix for extrapolation values (here vector) * '* y y-value of solution at x * '* n number of equations in the system * '* epsabs absolute error bound * '* epsrel relative error bound * '* x x-value * '* h local step sizes * '* h0 * '* * '* Output parameters: * '* ================== * '* ahead *ahead = 1: Step is acceptable * '* *ahead != 1: Step not acceptable * '* index pointer to the arrays bufol and fhilf * '* x final x-value * '* h final step size * '* * '* global names used: * '* ================== * '* extmax, bufol, XMin, XMax, Copyvector. * '************************************************************************ DIM column AS INTEGER 'column of extrapolation scheme ' diffmax, maximal difference of two columns in scheme ' fhilfi1, column in extrapolation scheme ' fhilfi2, column to the right of fhilfi1 ' bufolr, element in Bulirsch sequence ' bufoli, idem ' ymax, maximal value in a column ' help aux variable ymax = 0# maxcol = Min(row + 1, extmax) ahead = 0 'False FOR column = 2 TO maxcol IF ahead <> 0 THEN EXIT SUB index = Min(11 - column, row - column + 3) diffmax = 0# FOR i = 0 TO n - 1 fhilfi1 = fhilf((index - 1) * n + i) fhilfi2 = fhilf((index - 2) * n + i) bufolr = 1# * bufol(row) bufoli = 1# * bufol(index - 2) fhilf((index - 2) * n + i) = fhilfi1 + (fhilfi1 - fhilfi2) / ((bufolr / bufoli) ^ 2 - 1#) fhilfi1 = fhilf((index - 1) * n + i) fhilfi2 = fhilf((index - 2) * n + i) ymax = XMax(ymax, ABS(fhilfi2)) diffmax = XMax(diffmax, ABS(fhilfi2 - fhilfi1)) NEXT i IF diffmax < epsrel * ymax + epsabs THEN 'Step acceptable ? x = x + h0 'Copyvector(y, fhilf(index - 2, n) FOR i = 0 TO n - 1 y(i) = fhilf((index - 2) * n + i) NEXT i help = 1# * (column - extmax) 'Step size for next step h = .9# * h * (.6 ^ help) row = -1 ahead = 1 'True END IF NEXT column END SUB 'Gauss decomposition ....................... SUB gaudec (n, mat() AS DOUBLE, lumat() AS DOUBLE, perm() AS INTEGER, signd AS INTEGER, rc AS INTEGER) ' n, size of matrix .................. ' mat() as Double, Input matrix .................... ' lumat() as Double, matrix decomposition ............ ' perm() as integer, row interchanges ................ ' signd as integer, sign of perm .................... ' rc as integer error code ...................... '*====================================================================* '* * '* gaudec decomposes a nonsingular n x n matrix into a product of a * '* unit lower and an upper triangular matrix. Both triangular factors* '* are stored in lumat (minus the unit diagonal, of course). * '* * '* ------------------------------------------------------------------ * '* * '* Input parameter: * '* ================ * '* n integer; (n > 0) * '* Dimension of mat and lumat, * '* size of b , x and perm. * '* mat type pMAT; * '* pointer to original system matrix in vector form. * '* * '* Output parameters: * '* ================== * '* lumat type pMAT; * '* pointer to LU factorization * '* perm type pIVEC; * '* pointer to row permutation vector for lumat * '* signd TYPE integer; * '* sign of perm. The determinant of mat can be computed * '* as the product of the diagonal entreis of lumat times* '* signd. * '* * '* Return value (rc): * '* ================= * '* = 0 all ok * '* = 1 n < 1 or invalid input * '* = 2 lack of memory * '* = 3 Matrix is singular * '* = 4 Matrix numerically singular * '* * '*====================================================================* DIM d(n) 'scaling vector for pivoting DIM MACHEPS AS DOUBLE MACHEPS = 2.25D-16 IF n < 1 THEN rc = 1 'Invalid parameters EXIT SUB END IF 'copy mat to lumat FOR i = 0 TO n - 1 FOR j = 0 TO n - 1 lumat(i * n + j) = mat(i * n + j) NEXT j NEXT i FOR i = 0 TO n - 1 perm(i) = i 'Initialize perm zmax = 0# FOR j = 0 TO n - 1 'find row maxima tmp = ABS(lumat(i * n + j)) IF tmp > zmax THEN zmax = tmp NEXT j IF zmax = 0# THEN 'mat is singular rc = 3 EXIT SUB END IF d(i) = 1# / zmax NEXT i signd = 1 'initialize sign of perm FOR i = 0 TO n - 1 piv = ABS(lumat(i * n + i)) * d(i) j0 = i 'Search for pivot element FOR j = i + 1 TO n - 1 tmp = ABS(lumat(j * n + i)) * d(j) IF piv < tmp THEN piv = tmp 'Mark pivot element and j0 = j 'its location. END IF NEXT j IF piv < MACHEPS THEN 'If piv is too small, mat is signd = 0 'nearly singular. rc = 4 EXIT SUB END IF IF j0 <> i THEN signd = -signd 'update signd ISWAP perm(j0), perm(i) 'swap pivot entries XSWAP d(j0), d(i) 'swap scaling vector FOR j = 0 TO n - 1 XSWAP lumat(j0 * n + j), lumat(i * n + j)'swap j0-th and i-th rows of lumat NEXT j END IF FOR j = i + 1 TO n - 1 'Gauss elimination IF lumat(i * n + i) <> 0# THEN lumat(j * n + i) = lumat(j * n + i) / lumat(i * n + i) tmp = lumat(j * n + i) FOR m = i + 1 TO n - 1 lumat(j * n + m) = lumat(j * n + m) - tmp * lumat(i * n + m) NEXT m END IF NEXT j NEXT i rc = 0'all ok END SUB '***************************************************************** '* Gauss algorithm for solving linear equations * '* ------------------------------------------------------------- * '* REF.: "Numerical Algorithms with C, by Gisela Engeln-Muellges * '* and Frank Uhlig, Springer-Verlag, 1996". * '* * '* TPW Release By J-P Moreau, Paris. * '* ------------------------------------------------------------- * '* REF.: "Numerical Algorithms with C, by Gisela Engeln-Muellges * '* and Frank Uhlig, Springer-Verlag, 1996". * '***************************************************************** 'Gauss solution ............................ SUB gausol (n, lumat() AS DOUBLE, perm() AS INTEGER, b(), x(), rc AS INTEGER) ' n, size of matrix .................. ' lumat() as double, decomposed matrix (LU) .......... ' perm() as integer, row permutation vector .......... ' b(), Right hand side ................. ' x(), solution ........................ ' rc as integer error code ...................... '*====================================================================* '* * '* gausol finds the solution x of the linear system lumat * x = b * '* for the product matrix lumat, that describes an LU decomposition, * '* as produced by gaudec. * '* * '*====================================================================* '* * '* Input parameters: * '* ================ * '* n integer; (n > 0) * '* Dimension of lumat, * '* lumat (0:n-1,0:n-1) matrix, here a (0:n*n-1) vector (1). * '* LU factorization matrix, as produced from gaudec * '* perm row permutation vector for lumat (integer). * '* b Right hand side of the system (size 0:n-1). * '* * '* Output parameter: * '* ================ * '* x Solution (0:n-1) vector. * '* * '* Return value (rc): * '* ================= * '* = 0 all ok * '* = 1 n < 1 or other invalid input parameter * '* = 3 improper LU decomposition ( zero diagonal entry) * '* * '* (1) element lumat(i,j) is here lumat(i*n+j). * '*====================================================================* DIM MACHEPS AS DOUBLE MACHEPS = 2.25D-16 IF n < 1 THEN rc = 1 'Invalid input parameter EXIT SUB END IF FOR k = 0 TO n - 1 'update b x(k) = b(perm(k)) FOR j = 0 TO k - 1 x(k) = x(k) - lumat(k * n + j) * x(j) NEXT j NEXT k FOR k = n - 1 TO 0 STEP -1 'back substitute sum = 0# FOR j = k + 1 TO n - 1 sum = sum + lumat(k * n + j) * x(j) NEXT j IF ABS(lumat(k * n + k)) < MACHEPS THEN rc = 3 EXIT SUB END IF x(k) = (x(k) - sum) / lumat(k * n + k) NEXT k rc = 0 'all ok END SUB 'Gauss algorithm for solving linear equations SUB gauss (mode, n, mat() AS DOUBLE, lumat() AS DOUBLE, perm() AS INTEGER, b(), x(), signd AS INTEGER, code AS INTEGER) ' mode, Modus: 0, 1, 2, 3 ............... ' n, Dimension of matrix ............. ' mat() as double, Input matrix (here a vector) .... ' lumat() as double, LU decomposition (here a vector). ' perm() as integer, row remutation vector ........... ' b(), right hand side ................. ' x(), solution of the system .......... ' signd as integer, sign of the permutation ......... ' code as integer '*====================================================================* '* * '* The procedure gauss solves a linear system : mat * x = b. * '* Here mat is the nonsingular system matrix, b the right hand side * '* of the system and x the solution vector. * '* * '* gauss uses the Gauss algorithm and computes a triangular factori- * '* zation of mat and scaled column pivot search. (Crout method with * '* row swaps). * '* * '* ------------------------------------------------------------------ * '* * '* Application: * '* ============ * '* Solve general linear system with a nonsingular coefficient * '* matrix. * '* * '* ------------------------------------------------------------------ * '* * '* Control parameter: * '* ================== * '* mode integer; * '* calling modus for gauss: * '* = 0 Find factorization and solve linear system * '* = 1 Find factorization only. * '* = 2 Solve linear system only; the factorization is * '* already available in lumat. This saves work when * '* solving a linear system repeatedly for several right * '* hand sides and the same system matrix such as when * '* inverting the matrix. * '* = 3 as under 2, additionally we improve the solution * '* via iterative refinement (not available here). * '* * '* Input parameters: * '* ================ * '* n integer; (n > 0) * '* Dimension of mat and lumat, * '* size of the vector b, the right hand side, the * '* solution x and the permutation vector perm. * '* mat (n,n) matrix of the linear system (here a n*n vector)* '* lumat (for mode = 2, 3) LU factors of mat (n*n vector). * '* lumat can be stored in the space of mat. * '* perm (for mode = 2, 3) * '* Permutation vector, of the row interchanges in lumat * '* b (for mode = 0, 2, 3) * '* Right hand side of the system. * '* signd type Integer; (for mode = 2, 3) * '* sign of the permutation in perm; the determinant of * '* mat can be computed as the product of the diagonal * '* entries of lumat times signd. * '* * '* Output parameters: * '* ================== * '* lumat n*n vector; (for mode = 0, 1) * '* LU factorization of mat. * '* perm n*n vector; (for mode = 0, 1) * '* row permutation vektor * '* x (0:n-1) vector (for mode = 0, 2, 3) * '* Solution vector. * '* signd type Integer; (for mode = 0, 1) * '* sign of perm. * '* * '* Return value (code): * '* =================== * '* =-1 Max. number (MAXITER) of iterative refinements * '* reached (MAXITER) while mod = 3 * '* = 0 all ok * '* = 1 n < 1 or other invalid input * '* = 2 lack of memory * '* = 3 Matrix singular * '* = 4 Matrix numerically singular * '* = 5 incorrect call * '* * '* ------------------------------------------------------------------ * '* * '* Subroutines used: * '* ================ * '* * '* gaudec: determines LU decomposition * '* gausol: solves the linear system * '* * '*====================================================================} DIM rc AS INTEGER IF n < 1 THEN code = 1 EXIT SUB END IF IF mode = 0 THEN 'Find factorization and solve system ............. gaudec n, mat(), lumat(), perm(), signd, rc IF rc = 0 THEN gausol n, lumat(), perm(), b(), x(), rc code = rc EXIT SUB ELSE code = rc EXIT SUB END IF END IF IF mode = 1 THEN 'Find factorization only ......................... gaudec n, mat(), lumat(), perm(), signd, rc code = rc EXIT SUB END IF IF mode = 2 THEN 'Solve only ...................................... gausol n, lumat(), perm(), b(), x(), rc code = rc EXIT SUB END IF IF mode = 3 THEN 'solve and then use iterative refinement ......... PRINT " fgauss: gausoli not implemented." code = 5 EXIT SUB END IF code = 5 'Wrong call END SUB 'Extrapolation method for 1st order DE systems. FUNCTION ibulstoe (x, xend, n, y(), epsabs, epsrel, h, hmax, neu, fmax AS INTEGER, aufrufe AS INTEGER) ' x: Double initial x-value/ final x-value ..... ' xend: Double desired end point .................. ' n: integer number of DEs ...................... ' y(n) initial y-value/ final y-value ..... ' epsabs, absolute error bound ............... ' epsrel: Double relative error bound ............... ' h: Double initial/final step size ............ ' hmax: Double maximal step size .................. ' neu: boolean use an outside given x ? ........... ' fmax: integer maximal # of calls of dgl ......... ' aufrufe: integer actual # of calls of dgl ........... '************************************************************************ '* Compute the solution of an intial value problem for a first order * '* system of ordinary differential equations * '* * '* y' = f(x,y) with initial value y(x0) = y0 * '* * '* by using the extrapolation method of Bulirsch-Stoer-Gragg. * '* The maximal extrapolation order is set internally depending on the * '* machine constant; the step size control follows the ideas in : * '* HALL, G.; WATT, J.M: Modern Numerical Methods for Ordinary * '* Differential Equations, Oxford 1976, * '* Clarendon Press, [HALL76] * '* REMARK : * '* The dynamic allocations from this function are only set free for a * '* new run or with improper input (fehler = TRUE or Return value >= 2). * '* Hence we advise the user to end work with this function by one call * '* with an improper value such as n = -2 in order to free up all storage* '* * '* Input parameters: * '* ================= * '* x x-value * '* y [0..n-1] vector with initial y-value at x * '* dgl pointer to a function which evaluates the right hand side * '* odf the DE system y' = f(x,y) * '* n number of equations in the system * '* xend desired final value for x, xend < x is allowed * '* h step size for next step * '* hmax maximal step size; hmax > 0 * '* epsabs\ error bounds for absolute and relative errors, both non * '* epsrel/ negative. * '* We test for the mixed error: * '* |lokaler Fehler| <= |y| * epsrel + epsabs. * '* If epsrel = 0, we test the absolute error; if epsabs = 0 we * '* test the relative error. epsrel will be set to 10 * machine * '* constant if it has been chosen smaller than this. * '* neu Flag, indicating whether this function is called to compute * '* on after one partial run that was cut short due to too many * '* evaluations: * '* neu = FALSE: continue with old values * '* neu = TRUE: start new * '* fmax upper bound for the number of allowed function evaluations * '* of the right hand side via dgl() * '* * '* Output parameters: * '* ================== * '* x final x-value of integration; usually x = xend * '* y [0..n-1] y-value vector at x * '* h final step size; should be used unchanged in the next call * '* if neu = TRUE; should be redefined for neu = TRUE. * '* aufrufe number of calls of dgl() * '* * '* Return value : * '* ============== * '* = 0: all ok; after a new setting for xend, the function can be used * '* again for the same problem, or all inputs may be changed * '* = 1: method has not reached xend after FMAX right hand side * '* evaluations; try another call with same parameters or change * '* error bounds. * '* = 2: step size less than 4 * machine constant. * '* For a future call, increase step size and error bounds. * '* = 3: epsabs or epsrel negative, or both zero. * '* = 4: xend = x * '* = 5: hmax negative. * '* = 6: n <= 0. * '* * '* global names used: * '* ================== * '* extrapol, MACHEPS, XMin, Min, XMax, copyvector, ABS. * '************************************************************************ DIM MACHEPS AS DOUBLE ' Smallest machine-dependant real DIM bufol(12) AS INTEGER 'The Bulirsch sequence DIM fhilf(12 * n) ' Index 0..7: Matrix for the extrapolation scheme ' Index 8..11: aux vectors ' x0 as Double x-value on call of dgl DIM row AS INTEGER ' row of extrapolation scheme ' absh, magnitude of step size ' absh0, magnitude of aux step size h0 ' h0, auxiliary step size ' hilf: Double auxiliary variable DIM count AS INTEGER ' Loop variable for mid-point rule ' index: integer Index for the vectors fhilf and bufol DIM ahead AS INTEGER ' Flag indicating acceptance of step DIM extmax AS INTEGER ' Maximum order of extrapolation DIM tmp(n), tmp1(n) ' Auxiliary vectors of size n for calls of dgl MACHEPS = 2.225D-16 'Initialize bufol array bufol(0) = 2: bufol(1) = 4: bufol(2) = 6: bufol(3) = 8: bufol(4) = 12 bufol(5) = 16: bufol(6) = 24: bufol(7) = 32: bufol(8) = 48 bufol(9) = 64: bufol(10) = 96: bufol(11) = 128 FOR i = 0 TO 11 FOR j = 0 TO n - 1 fhilf(i * n + j) = 0# NEXT j NEXT i h0 = 0# 'check input IF epsabs < 0# OR epsrel < 0# OR epsabs + epsrel <= 0# THEN ibulstoe = 3 EXIT FUNCTION ELSEIF xend = x THEN ibulstoe = 4 EXIT FUNCTION ELSEIF hmax <= 0# THEN ibulstoe = 5 EXIT FUNCTION END IF IF n <= 0 THEN ibulstoe = 6 EXIT FUNCTION END IF IF neu = 0 THEN 'new call or repeat call due to excessive evaluations ? x = fhilf(9 * n) 'copyvector(y, fhilf(8), n) FOR i = 0 TO n - 1 y(i) = fhilf(8 * n + i) NEXT i END ELSE row = -1 END IF aufrufe = 0 ahead = 1 'TRUE extmax = INT(-LOG(MACHEPS) / LOG(2#) / 7.5#) epsrel = XMax(epsrel, 10# * MACHEPS) 'main integration loop 15 IF neu <> 0 THEN 'new calculation IF ahead <> 0 THEN 'new step ? absh = ABS(h) absh = XMin(absh, hmax) h0 = XMin(absh, ABS(xend - x)) IF xend < x THEN h0 = -h0 absh0 = ABS(h0) IF absh0 <= 4# * MACHEPS * ABS(x) THEN ibulstoe = 0 EXIT FUNCTION 'normal exit END IF ahead = 0 'FALSE END IF 20 row = row + 1 'find step size for extrapolation h = h0 / bufol(row) x0 = x 'Euler step; save initial values 'copyvector(fhilf(8), y, n) FOR i = 0 TO n - 1 fhilf(8 * n + i) = y(i) NEXT i FOR i = 0 TO n - 1 tmp(i) = fhilf(8 * n + i) NEXT i dgl 0, n, x0, tmp(), tmp1() FOR i = 0 TO n - 1 fhilf(11 * n + i) = tmp1(i) NEXT i FOR i = 0 TO n - 1 fhilf(9 * n + i) = fhilf(8 * n + i) + h * fhilf(11 * n + i) NEXT i x0 = x0 + h FOR i = 0 TO n - 1 tmp(i) = fhilf(9 * n + i) NEXT i dgl 0, n, x0, tmp(), tmp1() FOR i = 0 TO n - 1 fhilf(11 * n + i) = tmp1(i) NEXT i 'use mid-point rule FOR count = 1 TO bufol(row) - 1 FOR i = 0 TO n - 1 fhilf(10 * n + i) = fhilf(8 * n + i) + 2# * h * fhilf(11 * n + i) NEXT i x0 = x0 + h FOR i = 0 TO n - 1 tmp(i) = fhilf(10 * n + i) NEXT i dgl 0, n, x0, tmp(), tmp1() FOR i = 0 TO n - 1 fhilf(11 * n + i) = tmp1(i) NEXT i FOR j = 8 TO 10 'store for next step 'copyvector(fhilf(j), fhilf(j + 1), n) FOR i = 0 TO n - 1 fhilf(j * n + i) = fhilf((j + 1) * n + i) NEXT i NEXT j NEXT count FOR i = 0 TO n - 1 tmp(i) = fhilf(9 * n + i) NEXT i dgl 0, n, x0, tmp(), tmp1() FOR i = 0 TO n - 1 fhilf(11 * n + i) = tmp1(i) NEXT i FOR i = 0 TO n - 1 'stabilize with trapezoidal rule fhilf(row * n + i) = .5# * (fhilf(9 * n + i) + fhilf(8 * n + i) + h * fhilf(11 * n + i)) NEXT i aufrufe = aufrufe + bufol(row) + 2 IF row = 0 THEN GOTO 20 'Extrapolation extrapol row, fhilf(), y(), n, epsabs, epsrel, x, h, h0, ahead, index, extmax, bufol() IF aufrufe >= fmax THEN ' too many calls ? ' => stop fhilf(9 * n) = x 'copyvector(fhilf(8), y, n) FOR i = 0 TO n - 1 fhilf(8 * n + i) = y(i) NEXT i x = x0 'copyvector(y, fhilf(index - 2), n) FOR i = 0 TO n - 1 y(i) = fhilf((index - 2) * n + i) NEXT i ibulstoe = 1 EXIT FUNCTION END IF END IF 'neu <> 0 IF ahead = 0 OR neu = 0 THEN 'do we need to repeat step ? neu = 1 'set flag for a new call IF row >= Min(7, extmax - 1) THEN ' store differently since the extrapolation scheme has at most} ' 11 rows, but we allow only 8 extrapolations. FOR j = 0 TO 6 'copyvector(fhilf(j), fhilf(j+1), n) FOR i = 0 TO n - 1 fhilf(j * n + i) = fhilf((j + 1) * n + i) NEXT i NEXT j END IF ' accuracy could not be reached from the whole extrapolation ' table; repeat step with smaller step size. IF row >= extmax + 2 THEN hilf = 1# * (row - extmax + 1) h0 = .9# * h * (.6 ^ hilf) IF ABS(h0) <= 4# * MACHEPS * ABS(x0) THEN 'step too small ibulstoe = 2 EXIT FUNCTION END IF row = -1 END IF END IF GOTO 15 'end of main loop END FUNCTION 'IBULSTOE SUB incvector (ziel(), quelle(), factor, n) '************************************************************************ '* add factor * quelle to the vector ziel * '* * '* global names used: * '* ================== * '* None * '************************************************************************ FOR i = 0 TO n - 1 ziel(i) = ziel(i) + factor * quelle(i) NEXT i END SUB SUB init0vector (vektor(), n) '************************************************************************ '* initialize the [0..n-1] vector vektor to be the zero vector * '* * '* global names used: * '* ================== * '* None * '************************************************************************ FOR i = 0 TO n - 1 vektor(i) = 0# NEXT i END SUB 'Predictor-corrector method for 1st order DE systems FUNCTION ipraekorr (x, y(), n, xend, h, epsabs, epsrel, fmax AS INTEGER, aufrufe AS INTEGER, hmax, neu) ' x initial/final x-value ........ ' y() initial value/ solution ...... ' n number of DEs ................ ' xend desired final x-value ........ ' h starting/final step size ..... ' epsabs absolute error bound ......... ' epsrel relative error bound ......... ' fmax maximal # of calls for dgl ... ' aufrufe actual # of calls of dgl ..... ' hmax maximal step size ............ ' neu delete old data ? ............ '************************************************************************ '* Solve a system of ordinary differential equations * '* * '* y' = f(x,y) with the initial condition y(x0) = y0 * '* * '* using the predictor-corrector method of Adams-Bashforth-Moulton. * '* The needed starting values are obtained by an initial run of Runge- * '* Kutta with the same order as the A-B-M method. * '* We work with automatic step size control with doubling/halfing of * '* the step size according to the subsequent error estimates. * '* * '* REMARK : * '* This procedure frees certain buffers only at the beginning of a * '* new run or when memory is not sufficient, or for improper input. * '* Hence we advise to terminate a connected run of this function with * '* one call with improper input such as n < 0 to clear all memory. * '* * '* Input parameters: * '* ================= * '* x x-value where a corresponding y-value is known * '* y [0..n-1] vector with the initial y-data * '* dgl pointer to a function f which evaluates the right hand side * '* of the DE system y' = f(x,y) * '* n number of DEs * '* xend x-value for which we want to find the solution; may be less * '* than x * '* h step size for next step; this is usually determined inside * '* this function * '* epsabs\ error bounds for absolute and relative errors; each >= 0 * '* epsrel/ We apply a mixed test : * '* |local error| <= |y| * epsrel + epsabs. * '* For epsrel = 0, we test for absolute error only. * '* For epsabs = 0, we test for relative error. * '* Each error bound is internally set to 10 * machine constant * '* if it should be prescribed to be too small. * '* fmax upper bound for the allowed number of evaluations of the * '* right hand side dgl() od the DE system * '* hmax maximal step size; must be positive * '* neu If TRUE, we start off using Runge-Kutta, even if the earlier* '* integration counld be continued with another A-B-M step * '* * '* Output parameters: * '* ================== * '* x final x-value of integration (if run was successful: x=xend)* '* y [0..n-1] approximate solution vector at x * '* h final step size; should be used for start of next call. * '* If changed, please reset flag as well. * '* aufrufe counter for calls of dgl() * '* * '* Return value : * '* ============== * '* error code : * '* = 0: no error, xend was reached * '* = 1: after fmax calls of dgl() we have not reached xend: * '* a new call with unchanged parameters might be successful * '* (or try raising the error bounds) * '* = 2: the step size is less than 8 * machine constant. * '* Before subsequent calls increase h and the error bounds. * '* = 3: epsabs or epsrel negative, or both equal to zero. * '* = 4: xend = x * '* = 5: hmax negative * '* = 6: n <= 0 * '* = 7: lack of available memory * '* * '* global names used: * '* ================== * '* rkstart, abmschritt, initpraeko, MACHEPS, XMax, Xnormmax, copyvector * '************************************************************************ DIM MACHEPS AS DOUBLE 'Machine-dependant smallest real number DIM CHECK AS INTEGER 'check new step size in rkstart ' Vector with factors for error estimation: ' guess(0) for Runge-Kutta, guess(1) for Adams-Bashforth-Moulton DIM guess(2) 'local auxiliary vectors DIM f(6), tmp(n), tmp2(n), y1(n), y2(n), diff(n) DIM kisum(4) AS DOUBLE DIM rks AS INTEGER 'Result of irkstart DIM folgeaufruf AS INTEGER 'Flag indicating we can continue 'an earlier call with A-B-M method. MACHEPS = 2.25D-16 CHECK = 1 folgeaufruf = 0 guess(0) = 1# / 80# guess(1) = -19# / 270# ' ------------------- Check input data ------------------------------- IF epsabs < 0# OR epsrel < 0# OR epsabs + epsrel <= 0# THEN folgeaufruf = 0 ipraekorr = 3 EXIT FUNCTION END IF IF xend = x THEN folgeaufruf = 0 ipraekorr = 4 EXIT FUNCTION END IF IF hmax <= 0# THEN folgeaufruf = 0 ipraekorr = 5 EXIT FUNCTION END IF IF n <= 0 THEN folgeaufruf = 0 ipraekorr = 6 EXIT FUNCTION END IF ' -------------- Prepare integration loop --------------------------- epsrel = XMax(epsrel, 10# * MACHEPS) epsabs = XMax(epsabs, 10# * MACHEPS) aufrufe = 1 IF neu <> 0 THEN folgeaufruf = 0 'Force start with R-K step ? IF folgeaufruf <> 0 THEN 'Use values of previous call ? abmschritt x0, n, h, methode, aufrufe, f(), tmp(), kisum(), y1(), y2(), diff() ELSE 'very first call ? dgl 0, n, x, y(), tmp2() 'store beginning of starting FOR i = 0 TO n - 1 'values in f(1)..f(n) f(1 + i) = tmp2(i) NEXT i aufrufe = aufrufe + 1 hsave = h 'save starting step size xsave = x rks = irkstart(x, x0, y(), n, xend, h, hmax, CHECK, methode, aufrufe, f(), tmp(), kisum(), y1(), y2(), diff()) x = xsave IF rks <> 0 THEN h = hsave folgeaufruf = 0 ipraekorr = 0 EXIT FUNCTION END IF ' ---------- starting values available in ---------- ' ---------- f(1), f(2), f(3) and f(4) ---------- END IF ' ---------------------- Integration loop -------------------------- 10 IF aufrufe > fmax THEN 'excessive function calls ? x = x0 copyvector y(), y2(), n 'newest approximation for y() folgeaufruf = 1 ipraekorr = 1 'stop and report excessive calls EXIT FUNCTION END IF FOR i = 0 TO n - 1 'errro estimation diff(i) = guess(methode) * (y2(i) - y1(i)) NEXT i diffnorm = XNormMax(diff(), n) ynorm = XNormMax(y2(), n) IF diffnorm >= epsrel * ynorm + epsabs THEN 'error too large ? h = h * .5# 'halve the step size 'and repeat IF ABS(h) <= 8# * MACHEPS * ABS(x0) THEN 'step size too small ? folgeaufruf = 0 ipraekorr = 2 EXIT FUNCTION END IF xsave = x i = irkstart(x, x0, y(), n, xend, h, hmax, 0, methode, aufrufe, f(), tmp(), kisum(), y1(), y2(), diff()) x = xsave ELSE ' error not excessive ? ' step was successful, continue x = x0 ' on with the previous step size FOR i = 0 TO n - 1 ' add estimated error onto new y2(i) = y2(i) + diff(i) ' approximation NEXT i copyvector y(), y2(), n dgl 0, n, x0, y(), tmp2() ' correct last node for the next FOR i = 0 TO n - 1 ' A-B-M step. f(4 + i) = tmp2(i) NEXT i aufrufe = aufrufe + 1 'accuracy excessive ? IF diffnorm <= .02 * (epsrel * ynorm + epsabs) THEN h = h + h 'double step size hsave = h IF (h > 0# AND x0 >= xend) OR (h < ZERO AND x0 <= xend) THEN folgeaufruf = 0 ipraekorr = 0 'all done ! EXIT FUNCTION 'normal exit END IF ' Continue integration with doubled step size. ' First find a new set of starting values via Runge-Kutta. 'copyvector f(1), f(4), n FOR i = 0 TO n - 1 f(1 + i) = f(4 + i) NEXT i xsave = x IF irkstart(x, x0, y(), n, xend, h, hmax, CHECK, methode, aufrufe, f(), tmp(), kisum(), y1(), y2(), diff()) <> 0 THEN h = hsave folgeaufruf = 0 ipraekorr = 0 EXIT FUNCTION END IF x = xsave ELSE 'last step successful ? IF (h > 0# AND x0 >= xend) OR (h < 0# AND x0 <= xend) THEN folgeaufruf = 0 ipraekorr = 0 EXIT FUNCTION END IF IF (h > 0# AND x0 + h >= xend) OR (h < 0# AND x0 + h <= xend) THEN 'to xend ? ' The next step would lead beyond xend. Hence we do not use ' xend - x0 as the step size, and start anew using R-K. h = xend - x0 hsave = h 'copyvector f(1), f(4), n old starting values to be the FOR i = 0 TO n - 1 ' first value of the new one. f(1 + i) = f(4 + i) NEXT i xsave = x IF irkstart(x, x0, y(), n, xend, h, hmax, CHECK, methode, aufrufe, f(), tmp(), kisum(), y1(), y2(), diff()) <> 0 THEN h = hsave folgeaufruf = 0 ipraekorr = 0 EXIT FUNCTION END IF x = xsave ELSE 'perform one A-B-M step abmschritt x0, n, h, methode, aufrufe, f(), tmp(), kisum(), y1(), y2(), diff() END IF END IF END IF GOTO 10 END FUNCTION 'ipraekorr FUNCTION irkstart (x, x0, y(), n, xend, h, hmax, newstep, methode, aufrufe AS INTEGER, f(), tmp(), kisum() AS DOUBLE, y1(), y2(), diff()) '************************************************************************ '* Find the starting values using the Runge-Kutta method as needed in * '* ipraekorr for the Adams-Bashforth-Moulton method * '* * '* Input parameters: * '* ================= * '* x x-value for start of integration * '* y() Initial value for the solution at x * '* n number of DE equations * '* xend x-value where solution is wanted; xend may be less than x * '* h Step size * '* hmax maximal step size, must be positive * '* newstep <> 0: check new step size for properness * '* = 0: do not check new step size * '* * '* Output parameters: * '* ================== * '* x0 x-value until which we have integrated * '* h step size for next integration * '* methode value 0; hence the error estimate in praekorr is performed * '* for the factor for Runge-Kutta values. * '* aufrufe actual number of calls of dgl * '* hilf Record, which contains info of the solution: * '* y1: approximate solution needed for error estimation* '* (Step size 3*h) * '* y2: actual approximate solution * '* f[2]..f[4]: starting values * '* The entries in kisum and tmp represent auxiliary values. * '* * '* Return value: * '* ============= * '* = 0: all ok * '* = 1: new step size too small relative to machine constant * '* * '* Note: * '* ===== * '* x is modified during the process and must be saved and restored by * '* the calling subroutine, if necessary. * '* * '* global names used: * '* ================== * '* rkschritt, incvector, addvector, MACHEPS, Min, copyvector. * '************************************************************************ DIM MACHEPS AS DOUBLE DIM tmp2(n) MACHEPS = 2.25D-16 IF newstep <> 0 THEN 'check new step size for properness h = XMin(h, hmax) h = XMin(ABS(h), ABS(xend - x) / 3#) IF xend <= x THEN h = -h IF ABS(h) <= 8# * MACHEPS * ABS(x) THEN irkstart = 1 EXIT FUNCTION END IF END IF x0 = x 'save initial value copyvector y2(), y(), n 'y2 <- y FOR j = 2 TO 4 'three steps with steps h 'accumulate kisum rkschritt x, y2(), n, f(), tmp(), kisum(), y1(), y2(), diff(), h x = x + h incvector y2(), kisum(), h, n 'y2 = y2 + h*kisum dgl 0, n, x, y2(), tmp2() 'copy remainder of starting FOR i = 0 TO n - 1 'values in f(2)..f(4). f(j + i) = tmp2(i) NEXT i NEXT j ' after three steps with size h, we now perform one step with 3*h ' result put into y1. 'compute kisum rkschritt x0, y(), n, f(), tmp(), kisum(), y1(), y2(), diff(), 3# * h addvector y1(), y(), kisum(), 3# * h, n x0 = x0 + 3# * h methode = 0 aufrufe = aufrufe + 13 '13 calls of dgl for the starting values irkstart = 0 END FUNCTION ' ------------------------------------------------------------------- ' 'Shooting method for boundary value problem of 1st order' FUNCTION irwp (A, b, h, ystart(), n, awpnumm AS INTEGER, epsawp, epsrb, fmax AS INTEGER, itmax, actiter AS INTEGER) ' a, left end point ..................... ' b, right end point .................... ' h, starting step size ................. ' ystart(), initial approximation or solution .. ' initial value problem y(a) ........ ' n, number of differntial equations .... ' awpnumm as integer, Number of desired IVP solver ....... ' epsawp, error bound for initial value ' problem ............................ ' epsrb, error bound for boundary value ' problem ............................ ' fmax as integer, maximal number of calls of dgl ..... ' itmax, maximal number of Newton iterations ' actiter AS INTEGER actual number of Newton steps ...... ' return value error code ......................... '************************************************************************ '* This function solves the first order boundary value problem * '* * '* y' = F(x,y), a <= x <= b, R(y(a),y(b)) = 0. * '* * '* It uses the shooting method starting with an approximate ystart for * '* the associated initial value y(a), from which it constructs an * '* approximate solution using as initial value solver the function awp. * '* The nonlinear system arising in the shooting method is solved using * '* Newton's method. * '* * '* Input parameters: * '* ================= * '* a left end point * '* b right end point (b > a) * '* h suitable step size for solving the associated initial value * '* problem for the shooting method approximately * '* ystart (0:n-1) vector, initial approximation for y(a) * '* n number of differential equations * '* awpnumm label for the desired IVP solver in the shooting method: * '* = 0: Runge-Kutta embedding formula of 4/5th order; England * '* formula from awp. * '* = 1: Predictor-corrector method of order 4 by Adams- * '* Bashforth-Moulton (from ipraekorr). * '* = 2: Runge-Kutta embedding formula of 7/8th order (from * '* ibulstoe). * '* epsawp desired accuracy for the solution of the associated initial * '* value problem * '* epsrb desired accuracy for which the approximation ystart for * '* the initial value y(a) of a solution of the boundary value * '* problem should satisfy the boundary condition R. * '* fmax upper bound for number of calls of dgl() in the system of * '* differential equations when solving the associated initial * '* value problem * '* itmax upper bound for number of Newton iterations when solving * '* the nonlinear systems in the shooting method * '* * '* Output parameters: * '* ================== * '* ystart (0..n-1) approximation for the initial value y(a) of a * '* solution y of the boundary problem * '* actiter number of Newton iterations performed * '* * '* Return value : * '* ============== * '* error code * '* = 0: all is ok * '* = 1: at least one of the accuracy bounds epsawp, epsrb too small * '* = 2: b <= a * '* = 3: h <= 0 * '* = 4: n <= 0 * '* = 5: improper input for awpnumm * '* = 6: maximal number of allowed function evaluations exceeded * '* = 7: actiter > itmax. number of allowd Newton steps exceeded without * '* finding a suitable value for ystart. * '* = 8: The Jacobi matrix for Newton iterations is singular. * '* > 9: error in one of the IVP solvers at first node * '* > 19: error in one of the IVP solvers at second node * '* * '************************************************************************ DIM Debug AS INTEGER ' if True, print debugging infos DIM yk(n) ' (0..n-1) vector with value of solution ' at right end point DIM yaj(n) ' (0..n-1) vector of modified ystart used ' for forming the Jacobi matrix DIM R(n) ' (0..n-1) vector with left boundary ' condition, which serves as the right hand ' side for the linear system in the Newton ' method. DIM rj(n) ' (0..n-1) vector with left end boundary ' condition of the modified initial value ' problem. DIM d(n) ' (0..n-1) vector, the solution of the ' linear system. DIM amat(n * n) ' (0..n-1,0..n-1) array of the Jacobi ' matrix for the Newton step (here a vector) DIM lumat(n * n) AS DOUBLE ' LU Decomposition of amat (here a vector) ' xk, left end point for awp ' hk, desired step size for awp ' epsabs, absolute accuracy desired in awp ' epsrel, relative accuracy in awp ' rmax, Maximum norm of the left boundary condition ' delta Aux variable for Jacobi matrix DIM mach AS DOUBLE ' accuracy bounds depending on computer DIM mach2 AS DOUBLE ' DIM aufrufe AS INTEGER ' number of function evaluations in awp DIM pivot(n) AS INTEGER ' Permutation vector for Gaussian elimin. ' jacobi, counter for columns of Jacobi matrix ' mark error code from gauss DIM signdet AS INTEGER DIM fehler AS INTEGER ' error code in awp CONST SINGU = 1 ' Name for the return value of gauss, ' which indicates that the matrix is ' singular CONST iengl45 = 6 ' awp shall use the England embedding ' formula of 4th and 5th order CONST VONVORNE = 1 CONST NOTSAVE = 0 CONST HALF = .5# DIM MACHEPS AS DOUBLE MACHEPS = 2.25D-16 mach = MACHEPS ^ .75 mach2 = 100# * MACHEPS Debug = 0: irwp = 0 IF (epsrb < mach2) THEN irwp = 1 'check input IF (b <= A) THEN irwp = 2 IF (h <= mach2 * ABS(A)) THEN irwp = 3 IF (n <= 0) THEN irwp = 4 IF (awpnumm < 1 OR awpnumm > 5) THEN irwp = 5 actiter = 0 epsabs = HALF * epsawp epsrel = epsabs IF (epsrel < .0000000001#) THEN epsrel2 = .0000000001# ELSE epsrel2 = epsrel END IF '********************************************************************** '* If ystart proves to be a sufficiently good approximation for y(a) * '* in the following loop, we report the success and check the * '* remaining input data in awp. * '********************************************************************** 'main loop WHILE (1 = 1) IF (Debug <> 0) THEN zeig "Y_start:", ystart(), n copyvector yk(), ystart(), n xk = A hk = h IF awpnumm = 0 THEN CALL awp(xk, b, 0, n, yk(), epsabs, epsrel, hk, iengl45, fmax, aufrufe, fehler) IF awpnumm = 1 THEN fehler = ipraekorr(xk, yk(), n, b, hk, epsabs, epsrel, fmax, aufrufe, 1#, VONVORNE) IF awpnumm = 2 THEN fehler = ibulstoe(xk, b, n, yk(), epsabs, epsrel, hk, 1#, VONVORNE, fmax, aufrufe) IF (fehler <> 0) THEN irwp = 10 + fehler 'error in IVP solver ? EXIT FUNCTION END IF IF (Debug <> 0) THEN zeig "Yk:", yk(), n randbed ystart(), yk(), R() rmax = 0# FOR i = 0 TO n - 1 rmax = XMax(rmax, ABS(R(i))) NEXT i IF (rmax < epsrb) THEN 'boundary condition satisfied ? irwp = 0 'report success (normal exit) EXIT FUNCTION END IF actiter = actiter + 1 IF actiter > itmax THEN irwp = 6 'after itmax Newton steps ? EXIT FUNCTION 'report failure END IF 'find a better approximation ystart for y(a) by using the 'Newton method with a Jacobi matrix that is approximated by one- 'sides difference quotients. FOR jacobi = 0 TO n - 1 copyvector yk(), ystart(), n copyvector yaj(), ystart(), n IF (ABS(yk(jacobi)) < mach2) THEN yk(jacobi) = yk(jacobi) + mach delta = 1# / mach ELSE yk(jacobi) = yk(jacobi) * (1# + mach) delta = 1# / (mach * yk(jacobi)) END IF yaj(jacobi) = yk(jacobi) xk = A hk = h IF (Debug <> 0) THEN zeig "Yaj:", yaj(), n IF awpnumm = 0 THEN CALL awp(xk, b, 0, n, yk(), epsabs, epsrel, hk, iengl45, fmax, aufrufe, fehler) IF awpnumm = 1 THEN fehler = ipraekorr(xk, yk(), n, b, hk, epsabs, epsrel, fmax, aufrufe, 1#, VONVORNE) IF awpnumm = 2 THEN fehler = ibulstoe(xk, b, n, yk(), epsabs, epsrel, hk, 1#, VONVORNE, fmax, aufrufe) IF (fehler <> 0) THEN 'error in IVP solver ? irwp = 20 + fehler EXIT FUNCTION END IF IF (Debug <> 0) THEN zeig "Yk:", yk(), n randbed yaj(), yk(), rj() FOR i = 0 TO n - 1 amat(i * n + jacobi) = (rj(i) - R(i)) * delta NEXT i NEXT jacobi gauss 0, n, amat(), lumat(), pivot(), R(), d(), signdet, mark IF (mark = SINGU) THEN 'Jacobi matrix singular ? irwp = 8 'return error EXIT FUNCTION END IF IF (Debug <> 0) THEN zeig "D:", d(), n FOR i = 0 TO n - 1 'correct ystart ystart(i) = ystart(i) - d(i) NEXT i WEND 'while True END FUNCTION 'IRWP SUB ISWAP (ia, ib) itmp = ib: ib = ia: ia = itmp END SUB FUNCTION Max (ia, ib) IF ia > ib THEN Max = ia ELSE Max = ib END IF END FUNCTION FUNCTION Min (ia, ib) IF ia < ib THEN Min = ia ELSE Min = ib END IF END FUNCTION ' embedding formulas of Prince-Dormand of 4./5. order ......................... SUB prdo45 (x, y(), bspn AS INTEGER, n, h, y4(), y5(), steif1 AS INTEGER, steifanz AS INTEGER, steif2 AS INTEGER) ' double x, starting point of integration ..... ' y(0:n-1) initial value at x ................ ' integer bspn, # example ......................... ' n number of DEs ..................... ' double h, step size ......................... ' y4(0:n-1), solution of 4th order at x+h ...... ' y5(0:n-1) solution of 5th order at x+h ...... ' auxiliary flags ' integer steif1, steifanz,steif2 ' auxiliary vectors DIM yhilf(n) DIM k1(n) AS DOUBLE, k2(n) AS DOUBLE DIM k3(n) AS DOUBLE, k4(n) AS DOUBLE DIM k5(n) AS DOUBLE, k6(n) AS DOUBLE DIM k7(n) AS DOUBLE DIM g6(n) AS DOUBLE, g7(n) AS DOUBLE '************************************************************************ '* Compute 4th and 5th order approximates y4, y5 at x + h starting with * '* a solution y at x by using the Prince-Dormand embedding formulas on * '* the first order system of n differential equations y' = f(x,y) , as * '* supplied by dgl(). * '* Simultaneously we perform two tests for stiffness whose results are * '* stored in steif1 and steif2. * '* * '* Input parameters: * '* ================= * '* x initial x-value * '* y y-values at x (pVEC) * '* n number of differential equations * '* dgl function that evaluates the right hand side of the system * '* y' = f(x,y) * '* h step size * '* * '* yhilf, k1..k7,g6,g7: auxiliary vectors. * '* * '* Output parameters: * '* ================== * '* y4 4th order approximation for y at x + h * '* y5 5th order approximation for y at x + h * '* * '***********************************************************************} DIM steifa AS INTEGER 'Flag which is set if the second test for stiffness 'Shampine und Hiebert) is positive; otherwise the 'flag is erased. CALL dgl(bspn, n, x, y(), k1()) FOR i = 0 TO n - 1 yhilf(i) = y(i) + .2 * h * k1(i) NEXT i CALL dgl(bspn, n, x + .2 * h, yhilf(), k2()) FOR i = 0 TO n - 1 yhilf(i) = y(i) + .075 * h * (k1(i) + 3# * k2(i)) NEXT i CALL dgl(bspn, n, x + .3 * h, yhilf(), k3()) FOR i = 0 TO n - 1 yhilf(i) = y(i) + h / 45# * (44# * k1(i) - 168# * k2(i) + 160# * k3(i)) NEXT i CALL dgl(bspn, n, x + .8 * h, yhilf(), k4()) FOR i = 0 TO n - 1 yhilf(i) = y(i) + h / 6561# * (19372# * k1(i) - 76080# * k2(i) + 64448# * k3(i) - 1908# * k4(i)) NEXT i CALL dgl(bspn, n, x + (8# / 9#) * h, yhilf(), k5()) FOR i = 0 TO n - 1 g6(i) = y(i) + h / 167904# * (477901# * k1(i) - 1806240# * k2(i) + 1495424# * k3(i) + 46746# * k4(i) - 45927# * k5(i)) NEXT i CALL dgl(bspn, n, x + h, g6(), k6()) FOR i = 0 TO n - 1 g7(i) = y(i) + h / 142464# * (12985# * k1(i) + 64000# * k3(i) + 92750# * k4(i) - 45927# * k5(i) + 18656# * k6(i)) NEXT i CALL dgl(bspn, n, x + h, g7(), k7()) FOR i = 0 TO n - 1 y5(i) = g7(i) y4(i) = y(i) + h / 21369600# * (1921409# * k1(i) + 969088# * k3(i) + 13122270# * k4(i) - 5802111# * k5(i) + 1902912# * k6(i) + 534240# * k7(i)) NEXT i ' Test for stiffness via dominant eigenvalue IF DistMax(k7(), k6(), n) > 3.3 * DistMax(g7(), g6(), n) THEN steif1 = 1 ' one step in steffness test of Shampine & Hiebert FOR i = 0 TO n - 1 g6(i) = h * (2.2 * k2(i) + .13 * k4(i) + .144 * k5(i)) g7(i) = h * (2.134 * k1(i) + .24 * k3(i) + .1 * k6(i)) NEXT i IF DistMax(g6(), g7(), n) < DistMax(y4(), y5(), n) THEN steifa = 1 ELSE steifa = 0 END IF IF (steifa > 0) THEN steifanz = steifanz + 1 IF steifanz >= 3 THEN steif2 = 1 ELSE steifanz = 0 END IF END SUB SUB PrintHilf (f(), tmp(), kisum() AS DOUBLE, y1(), y2(), diff(), n) PRINT " Record hilf" printvec " Vector f:", 5, f() printvec " Vector tmp:", n, tmp() printvec " Vector ki_sum:", 4, kisum() printvec " Vector y1:", n, y1() printvec " Vector y2:", n, y2() printvec " Vector diff:", n, diff() END SUB 'print a REAL square matrix with Name1 (debug only) SUB PrintMat (Name1 AS STRING, n, xmat()) ' Name1 Matrix caption ' n Size of matrix ' double xmat(0:n-1,0:n-1) stored in a vector(n*n) ' matrix to be printed PRINT " "; Name1 FOR i = 0 TO n - 1 FOR j = 0 TO n - 1 PRINT " "; xmat(n * i + j); NEXT j PRINT NEXT i END SUB 'debug only SUB printvec (Name1 AS STRING, n, v()) PRINT Name1 FOR i = 0 TO n - 1 PRINT " "; v(i); NEXT i PRINT END SUB 'Boundary Conditions '=================== SUB randbed (ya(), yb(), R()) R(0) = ya(0) R(1) = yb(0) END SUB SUB rkschritt (x0, y0(), n, f(), tmp(), kisum() AS DOUBLE, y1(), y2(), diff(), h) '************************************************************************ '* perform one Runge-Kutta integration * '* * '* Input parameters: * '* ================= * '* x0 x-value from wchich to integrate * '* y0 initial value of solution at x0 * '* n number of DEs * '* h step size for current Runge-Kutta integration * '* * '* Output parameter: * '* ================= * '* hilf stores the result in kisum. The entries in y1 and tmp have no * '* meaning. * '* * '* global names used: * '* ================== * '* incvector, addvector. * '************************************************************************ DIM A(4), a1(4) CONST HALF = .5# CONST ONE = 1# CONST ZERO = 0# ' ---- Coefficients of the classical Runge-Kutta method ------------- A(0) = ONE / 6#: A(1) = ONE / 3# A(2) = ONE / 3#: A(3) = ONE / 6# a1(0) = ZERO: a1(1) = HALF a1(2) = HALF: a1(3) = ONE 'use also as the b(j,s) dgl 0, n, x0, y0(), tmp() 'k(0) <- right hand side at (x0,y0) FOR i = 0 TO n - 1 'initialize kisum with kisum(i) = A(0) * tmp(i) 'first term NEXT i FOR j = 1 TO 3 'compute k(1)..k(3), add to kisum 'multiplied by factor addvector y1(), y0(), tmp(), a1(j) * h, n 'y1 <- y0 + 'a1(j) * h * k(j - 1) dgl 0, n, x0 + a1(j) * h, y1(), tmp() 'compute k(j) incvector kisum(), tmp(), A(j), n 'kisum = kisum + A(j) * k(j) NEXT j END SUB ' Runge-Kutta embedding formulas of 2nd, 3rd degree .................... SUB ruku23 (x, y(), bspn AS INTEGER, n, h, y2(), y3()) '************************************************************************ '* Compute 2nd and 3rd order approximates y2, y3 at x + h starting with * '* a solution y at x by using Runge-Kutta embedding formulas on the * '* first order system of n differential equations y' = f(x,y) , as * '* supplied by dgl(). * '* * '* Input parameters: * '* ================= * '* x x-value of left end point * '* y y-values at x * '* bspn # example * '* n number of differential equations * '* dgl function that evaluates the right hand side of the system * '* y' = f(x,y) * '* h step size * '* * '* yhilf,k1,k2,k3: auxiliary vectors defined in module awp. * '* * '* Output parameters: * '* ================== * '* y2 2nd order approximation for y at x + h * '* y3 3rd order approximation for y at x + h * '* * '************************************************************************ DIM yhilf(n) DIM k1(n) AS DOUBLE DIM k2(n) AS DOUBLE DIM k3(n) AS DOUBLE CALL dgl(bspn, n, x, y(), k1()) FOR i = 0 TO n - 1 yhilf(i) = y(i) + h * k1(i) NEXT i CALL dgl(bspn, n, x + h, yhilf(), k2()) FOR i = 0 TO n - 1 yhilf(i) = y(i) + .25 * h * (k1(i) + k2(i)) NEXT i CALL dgl(bspn, n, x + .5 * h, yhilf(), k3()) FOR i = 0 TO n - 1 y2(i) = y(i) + .5 * h * (k1(i) + k2(i)) y3(i) = y(i) + h / 6# * (k1(i) + k2(i) + 4# * k3(i)) NEXT i END SUB FUNCTION XMax (A, b) IF A >= b THEN XMax = A ELSE XMax = b END IF END FUNCTION FUNCTION XMin (A, b) IF A <= b THEN XMin = A ELSE XMin = b END IF END FUNCTION ' Find the maximum norm of a REAL vector ............................ FUNCTION XNormMax (vektor(), n) ' double vektor(0:n-1) vector ................. ' integer n length of vector ....... ' ************************************************************************ ' * Return the maximum norm of a [0..n-1] vector v * ' * * ' ************************************************************************ DIM norm AS DOUBLE ' local max. ' double betrag ' magnitude of a component norm = 0# FOR i = 0 TO n - 1 betrag = ABS(vektor(i)) IF betrag > norm THEN norm = betrag NEXT i XNormMax = norm END FUNCTION SUB XSWAP (A, b) tmp = b: b = A: A = tmp END SUB 'print a vector of length n and its name (debug only) SUB zeig (s AS STRING, v(), n) PRINT s FOR i = 0 TO n - 1 PRINT " "; v(i); NEXT i PRINT END SUB