DECLARE SUB copyvector (a#(), b#(), n%) DECLARE SUB rkschritt (num%, x0#, y0#(), n%, f#(), tmp#(), kisum() AS DOUBLE, y1#(), y2#(), diff#(), h#) DECLARE FUNCTION XMax# (a#, b#) DECLARE FUNCTION XMin# (a#, b#) DECLARE FUNCTION XNormMax# (vektor#(), n%) DECLARE SUB abmschritt (num%, x0#, n%, h#, methode%, aufrufe AS INTEGER, f#(), tmp#(), kisum() AS DOUBLE, y1#(), y2#(), diff#()) DECLARE SUB init0vector (vektor#(), n%) DECLARE SUB incvector (ziel#(), quelle#(), factor#, n%) DECLARE SUB addvector (summe#(), summand1#(), summand2#(), factor#, n%) DECLARE SUB dgl (num%, n%, x#, Y#(), f#()) DECLARE FUNCTION ipraekorr% (num%, x#, Y#(), n%, xend#, h#, epsabs#, epsrel#, fmax AS INTEGER, aufrufe AS INTEGER, hmax#, neu%) DECLARE FUNCTION irkstart% (num%, x#, x0#, Y#(), n%, xend#, h#, hmax#, newstep%, methode%, aufrufe AS INTEGER, f#(), tmp#(), kisum() AS DOUBLE, y1#(), y2#(), diff#()) '******************************************************************* '* Solve an ordinary system of first order differential equations * '* using the predictor-corrector method of Adams-Bashforth-Moulton * '* --------------------------------------------------------------- * '* SAMPLE RUNS: * '* Example #1 * '* * '* # example (0 to 3) ..........: 0 * '* Number of equations in system: 2 * '* Starting x ..................: 0 * '* Ending x ....................: 2 * '* Absolute precision ..........: 1D-8 * '* Relative precision ..........: 1D-10 * '* Starting integration step ...: 0.1 * '* Y(0) = 0 * '* Y(1) = 0 * '* Max. number of function calls: 1000 * '* * '* D.E. System to solve: * '* y1' = y1 * y2 + cos(x) - 0.5 * sin(2.0*x) * '* y2' = y1 * y1 + y2 * y2 - (1 + sin(x)) * '* * '* x= 2 * '* Y(0) = 9.193162470825378D-02 * '* Y(1) = -1.363855037177237 * '* * '* Number of function calls: 566 * '* Final step size ........: 8.333333333333304D-03 * '* Error code .............: 0 * '* * '* Example #2 * '* * '* # example (0 to 3) ..........: 1 * '* Number of equations in system: 1 * '* Starting x ..................: 0 * '* Ending x ....................: 1 * '* Absolute precision ..........: 1D-8 * '* Relative precision ..........: 1D-10 * '* Starting integration step ...: 0.05 * '* Y(0) = 1 * '* Max. number of function calls: 500 * '* * '* D.E. System to solve: * '* y' = -y + x/((1+x)*(1+x)) * '* * '* x= 1 * '* Y(0) = .499999991411362 * '* * '* Number of function calls: 136 * '* Final step size ........: 3.333333333333333D-02 * '* Error code .............: 0 * '* * '* Example #3 * '* * '* # example (0 to 3) ..........: 3 * '* Number of equations in system: 5 * '* Starting x ..................: 0 * '* Ending x ....................: 2 * '* Absolute precision ..........: 1D-8 * '* Relative precision ..........: 1D-10 * '* Starting integration step ...: 0.1 * '* Y(0) = 1 * '* Y(1) = 1 * '* Y(2) = 1 * '* Y(3) = 1 * '* Y(4) = 1 * '* Max. number of function calls: 2000 * '* * '* D.E. System to solve: * '* y1' = y2 * '* y2' = y3 * '* y3' = y4 * '* y4' = y5 * '* y5' = (45 * y3 * y4 * y5 - 40 * y4 * y4 * y4) / (9 * y3 * y3) * '* * '* x= 2 * '* Y 1 = 6.708203932218345 * '* Y 2 = 5.341640785050235 * '* Y 3 = 2.414953413220991 * '* Y 4 = -1.448972049484087 * '* Y 5 = -1.448972043348845 * '* * '* Number of function calls: 1075 * '* Final step size ........: 1.041666666666593D-02 * '* Error code .............: 0 * '* * '* --------------------------------------------------------------- * '* Reference: [BIBLI 11]. * '* * '* Quick Basic Release 1.2 By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '* * '* Release 1.1 (12/15/06): added two more examples. * '* Release 1.2 (12/25/06): corrected bug in Sub ipraekorr * '* (added xsave). * '******************************************************************* 'Program Test_ABMOU DEFINT I-N DEFDBL A-H, O-Z OPTION BASE 0 DIM erreur AS INTEGER DIM bspn AS INTEGER DIM fmax AS INTEGER CLS PRINT INPUT " # Example (0 to 3) ......: ", bspn 'see sub dgl INPUT " Number of DEs ...........: ", n 'number of equations in system DIM Y(n) INPUT " Starting x value ........: ", xbegin INPUT " Ending x value ..........: ", xend INPUT " Absolute precision ......: ", epsabs INPUT " Relative precision ......: ", epsrel INPUT " Starting integration step: ", h FOR i = 0 TO n - 1 'initial value of Y PRINT " Y("; i; ") = "; : INPUT "", Y(i) NEXT i INPUT " Max. number of dgl calls : ", fmax 'call integration procedure with step size control erreur = ipraekorr(bspn, xbegin, Y(), n, xend, h, epsabs, epsrel, fmax, ncalls, 1#, 1) 'write results CLS PRINT PRINT " x= "; xend FOR i = 0 TO n - 1 PRINT " y"; i + 1; " = "; Y(i) NEXT i PRINT PRINT " Number of function calls: "; ncalls PRINT " Final step size ........: "; h PRINT " error code .............: "; erreur END 'of main program SUB abmschritt (num, 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: * '* ================= * '* num #example in dgl * '* 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 erreur 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 erreur 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 num, 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 num, 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 SUB copyvector (a(), b(), n) FOR i = 0 TO n - 1 a(i) = b(i) NEXT i END SUB 'System of DEs to integrate '========================== SUB dgl (bspn AS INTEGER, n, x, Y(), f()) IF bspn = 0 THEN f(0) = Y(0) * Y(1) + COS(x) - .5# * SIN(2# * x) f(1) = Y(0) * Y(0) + Y(1) * Y(1) - (1# + SIN(x)) ELSEIF bspn = 1 THEN f(0) = -Y(0) + x / ((1# + x) ^ 2) ELSEIF bspn = 2 THEN f(0) = Y(1) f(1) = -Y(0) ^ 3 ELSEIF bspn = 3 THEN f(0) = Y(1) f(1) = Y(2) f(2) = Y(3) f(3) = Y(4) f(4) = (45# * Y(2) * Y(3) * Y(4) - 40# * Y(3) ^ 3) / (9# * Y(2) ^ 2) END IF END SUB 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 (num, x, Y(), n, xend, h, epsabs, epsrel, fmax AS INTEGER, aufrufe AS INTEGER, hmax, neu) ' num # example of DEs ' 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 erreur bound ......... ' epsrel relative erreur 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: * '* ================= * '* num #example in dgl * '* x x-value where a corresponding y-value is known * '* y [0..n-1] vector with the initial y-data * '* 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 of 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 * '* * '* 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 irkstart ' Vector with factors for erreur estimation: ' guess(0) for Runge-Kutta, guess(1) for Adams-Bashforth-Moulton DIM guess(2) 'local auxiliary vectors DIM f(5 * n), 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 num, x0, n, h, methode, aufrufe, f(), tmp(), kisum(), y1(), y2(), diff() ELSE 'very first call ? dgl num, 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(num, 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 'erreur 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(num, x, x0, Y(), n, xend, h, hmax, 0, methode, aufrufe, f(), tmp(), kisum(), y1(), y2(), diff()) x = xsave ELSE ' erreur not excessive ? ' step was successful, continue x = x0 ' on with the previous step size FOR i = 0 TO n - 1 ' add estimated erreur onto new y2(i) = y2(i) + diff(i) ' approximation NEXT i copyvector Y(), y2(), n dgl num, 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(num, 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(num, 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 num, x0, n, h, methode, aufrufe, f(), tmp(), kisum(), y1(), y2(), diff() END IF END IF END IF GOTO 10 END FUNCTION FUNCTION irkstart (num, 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: * '* ================= * '* num # example in dgl * '* 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 * '* * '* 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 num, x, y2(), n, f(), tmp(), kisum(), y1(), y2(), diff(), h x = x + h incvector y2(), kisum(), h, n 'y2 = y2 + h*kisum dgl num, 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 num, 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 SUB rkschritt (num, x0, y0(), n, f(), tmp(), kisum() AS DOUBLE, y1(), y2(), diff(), h) '************************************************************************ '* perform one Runge-Kutta integration * '* * '* Input parameters: * '* ================= * '* num # example in dgl * '* 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 num, 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 num, 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 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