DECLARE SUB SWAP1 (a#, b#) DECLARE FUNCTION det# (n%, xmat#()) DECLARE FUNCTION XNormMax# (vektor#(), n%) DECLARE SUB settxt (bspn AS INTEGER, n%, dgltxt$()) DECLARE SUB gear4 (x#, xend#, bspn AS INTEGER, n%, y#(), epsabs#, epsrel#, h#,
fmax AS INTEGER, aufrufe AS INTEGER, fehler AS INTEGER)
DECLARE FUNCTION DistMax# (vector1#(), vector2#(), n%) DECLARE SUB ruku23 (x#, y#(), bspn AS INTEGER, n%, h#, y2#(), y3#()) DECLARE SUB dgl (bspn AS INTEGER, n%, x#, y#(), f#()) 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 SUB awp (x#, xend#, bspn AS INTEGER, n%, y#(), epsabs#, epsrel#, h#,
methode%, fmax AS INTEGER, aufrufe AS INTEGER, fehler AS INTEGER)
DECLARE FUNCTION XMax# (a#, b#) DECLARE FUNCTION XMin# (a#, b#) DECLARE SUB ISWAP (ia%, ib%) DECLARE SUB gauss (mode%, n%, xmat#(), xlumat#(), iperm%(), b#(), x#(),
signd AS INTEGER, code AS INTEGER)
DECLARE SUB gaudec (n%, xmat#(), xlumat#(), iperm%(), signd AS INTEGER, rc AS INTEGER) DECLARE SUB gausol (n%, xlumat#(), iperm%(), b#(), x#(), rc AS INTEGER) '************************************************************************ '* Solve a first order Stiff System of Differential Equations using * '* the implicit Gear method of order 4. * '* -------------------------------------------------------------------- * '* Mode of operation: * '* ================== * '* This program can solve two of the 11 examples of file t_dgls using * '* the implicit Gear method of order 4 (see file gear.pas). * '* To test other systems of DEs, please proceed as explained in file * '* t_dgls.pas. * '* * '* Inputs: * '* ======= * '* bspnummer Number of DE system from t_dgls.pas * '* epsabs desired absolute error bound * '* epsrel desired relative error bound * '* x0 left edge of integration * '* y0[0] \ known approximation for the solution at x0 * '* .. . > * '* y0[n-1] / * '* h initial step size * '* xend right endpoint of integration * '* fmax maximal number of calls of the right hand side * '* * '* The size n of the DE system is passed on from t_dgls.bas. * '* -------------------------------------------------------------------- * '* SAMPLE RUN * '* * '* Example #1: * '* (Solve set of differential equations (n=2): * '* f[0] = y[0] * y[1] + COS(x) - HALF * SIN(TWO * x); * '* f[1] = y[0] * y[0] + y[1] * y[1] - (ONE + SIN(x)); * '* Find values of f(0), f(1) at x=1.5). * '* * '* Input example number (0 to 11): 0 * '* abs. epsilon: 1e-6 * '* rel. epsilon: 1e-8 * '* x0: 0 * '* y0[0]: 0.5 * '* y0[1]: 0.5 * '* initial step size h: 0.0001 * '* right edge xend: 1.5 * '* maximal number of calls of right hand side: 6000 * '* * '* Input data: * '* ----------- * '* Example # 0 * '* Number of DEs = 2 * '* x0 = 0 * '* xend = 1.5 * '* epsabs = 0.000001 * '* epsrel = 0.00000001 * '* fmax = 6000 * '* h = 0.0001 * '* y0(0) = 0.5 * '* y0(1) = 0.5 * '* * '* Output data: * '* ------------ * '* error code from gear4: 0 * '* final local step size: 6.06783655109067E-02 * '* number of calls of right hand side: 360 * '* Integration stopped at x = 1.5 * '* * '* approximate solution y1(x) = 1.23598612893281 * '* approximate solution y2(x) = -0.104949617987246 * '* * '* Example #2: * '* (Solve set of differential equations (n=5): * '* f[0] = y[1]; * '* f[1] = y[2]; * '* f[2] = y[3]; * '* f[3] = y[4]; * '* f[4] = ((REAL)45.0 * y[2] * y[3] * y[4] - * '* (REAL)40.0 * y[3] * y[3] * y[3]) / (NINE * y[2] * y[2]); * '* Find values of f(0), ..., f(4) at x=1.5). * '* * '* Input example number (0 to 11): 3 * '* abs. epsilon: 1e-10 * '* rel. epsilon: 1e-10 * '* x0: 0 * '* y0[0]: 1 * '* y0[1]: 1 * '* y0[2]: 1 * '* y0[3]: 1 * '* y0[4]: 1 * '* initial step size h: 0.001 * '* right edge xend: 1.5 * '* maximal number of calls of right hand side: 6000 * '* * '* Input data: * '* ----------- * '* Example # 3 * '* Number of DEs = 5 * '* x0 = 0 * '* xend = 1.5 * '* epsabs = 0.0000000001 * '* epsrel = 0.0000000001 * '* fmax = 6000 * '* h = 0.001 * '* y0[0] = 1 * '* y0[1] = 1 * '* y0[2] = 1 * '* y0[3] = 1 * '* y0[4] = 1 * '* * '* Output data: * '* ------------ * '* error code from gear4: 0 * '* final local step size: 4.86347661993806E-03 * '* number of calls of right hand side: 3423 * '* Integration stopped at x = 1.5 * '* * '* approximate solution y1(x) = 4.36396102990278 * '* approximate solution y2(x) = 4.00000000763431 * '* approximate solution y3(x) = 2.82842715661993 * '* approximate solution y4(x) = 4.86163232512205E-08 * '* approximate solution y5(x) = -3.7712362229557 * '* * '* -------------------------------------------------------------------- * '* REF.: "Numerical Algorithms with C, By Gisela Engeln-Muellges * '* and Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11]. * '* * '* Quick Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '************************************************************************ DefDbl A-H, O-Z DefInt I-N Option Base 0 ' epsabs, absolute error bound ' epsrel, relative error bound ' x0 left edge of integration interval ' y0 (0:n-1)-vector initial values ' yex(0:n-1)-vector exact solution (when given) ' h initial, final step size ' xend right edge of integration interval Dim fmax As Integer ' maximal number of calls of right side in gear4 Dim aufrufe As Integer ' actual number of function calls Dim bspnummer As Integer ' # example ' n ' number of DEs in system (see t_dlgs.bas) Dim fehler As Integer ' error code from subroutine gear4 Dim dgltxt$(5) ' text of equations ' -------------------- read input -------------------- Cls Print INPUT " Example number (0 to 11): ", bspnummer If bspnummer <> 0 And bspnummer <> 3 Then Print Print " Example not registered." End 'stop program End If ' input absolute and relative errors INPUT " epsabs = ", epsabs INPUT " epsrel = ", epsrel ' input x starting value INPUT " x0 = ", x0 Call settxt(bspnummer, n, dgltxt$()) 'read text of equations 'and parameter n ' allocate memory for vectors y0, yex Dim y0(n) Dim yex(n) ' input initial values y0(i) ' y0(0) = 0.5: y0(1) = 0.5 (Ex. 0) For i = 0 To n - 1 PRINT " y0("; i; ") = "; : INPUT "", y0(i) Next i ' input initial step size INPUT " h = ", h ' input ending x value INPUT " xend = ", xend ' input maximum number of calls to right hand side INPUT " fmax = ", fmax ' ----------------- print input data ---------------------- Cls Print " ==============================================" Print " Solve a first order ordinary system of DEs " Print " using the implicit method of Gear of 4th order" Print " ==============================================" Print Print " System of DEs:" Print " ------------- " For i = 0 To n - 1 Print " "; dgltxt$(i) Next i Print " Example # "; bspnummer Print " Number of DEs = "; n Print " x0 = "; x0 Print " xend = "; xend Print " epsabs = "; epsabs Print " epsrel = "; epsrel Print " fmax = "; fmax Print " h = "; h For i = 0 To n - 1 Print " y0("; i; ") = "; y0(i) Next i ' ------------------- Solve system of DEs ------------------------------ gear4 x0, xend, bspnummer, n, y0(), epsabs, epsrel, h, fmax, aufrufe, fehler If fehler <> 0 Then 'something went wrong Print " Gear4: error nø "; 10 + fehler End 'stop program End If ' ---------------------- print results --------------------- Print Print " Output data:" Print " -----------" Print " error code from gear4: "; fehler Print " final local step size: "; h Print " number of calls of right hand side: "; aufrufe Print " Integration stopped at x = "; x0 Print For i = 0 To n - 1 Print " approximate solution y"; i + 1; "(x) = "; y0(i) Next i Print End 'of main program '************************************************************************ '* * '* 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 * '* dgl function that evaluates the right hand side of the system * '* y' = f(x,y) (see t_dgls) (Removed from list of parameters) * '* 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 MACH1 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 = 1.2E-16 MACH2 = 100 * MachEps MACH1 = MachEps ^ 0.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#, 0.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 < MACH1 * 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(0.5, 0.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 Function det(n, xmat()) ' Determinant ....................................................... ' integer n Dimension of the matrix ......... ' double xmat(0:n-1,0:n-1) matrix .......................... '*====================================================================* '* * '* det computes the determinant of an n x n real matrix xmat * '* * '*====================================================================* '* * '* Input parameter: * '* ================ * '* n integer; (n > 0) * '* Dimension of xmat * '* xmat matrix(0:n-1,0:n-1) * '* stored in a vector(n*n). * '* * '* Return value: * '* ============= * '* REAL Determinant of mat. * '* If the return value = 0, then the matrix is singular * '* or the storage is insufficient * '* * '*====================================================================* '* * '* subroutine used: * '* ================ * '* * '* gaudec (): LU decomposition of mat * '* * '*====================================================================* Dim rc As Integer, signd As Integer Dim iperm(n), xlu(n * n) Dim MachEps As Double Dim MAXROOT As Double If n < 1 Then 'n not valid det = 0# Exit Function End If MachEps = 1.2E-16 MAXROOT = 1E+16 Call gaudec(n, xmat(), xlu(), iperm(), 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(xlu(n * i + i)) > MAXROOT Then det = MAXROOT Exit Function Else tmpdet = tmpdet * xlu(n * i + i) 'compute det End If Next i det = tmpdet End Function '************************************************************************ '* * '* Test examples for methods to solve first order ordinary systems of * '* differential equations. * '* * '* We supply the right hand sides, their alpha-numeric description and * '* the exact analytic solution, if known. * '* * '* When running the main test program, the user can select among the * '* examples below. * '* * '* Quick Basic Release By J-P Moreau, Paris. * '* -------------------------------------------------------------------- * '* REF.: "Numerical Algorithms with C, by Gisela Engeln-Muellges * '* and Frank Uhlig, Springer-Verlag, 1996". * '************************************************************************ ' Note: here, only examples #0 and #3 are implemented in Quick Basic ' Sub dgl(bspn As Integer, n, x, y(), f()) If bspn = 0 Then f(0) = y(0) * y(1) + Cos(x) - 0.5 * Sin(2# * x) f(1) = y(0) ^ 2 + y(1) ^ 2 - (1# + Sin(x)) 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 '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) + 0.5 * h * k1(i) Next i Call dgl(bspn, n, x + 0.5 * h, yhilf(), k2()) For i = 0 To n - 1 yhilf(i) = y(i) + (0.25 * h * (k1(i) + k2(i))) Next i Call dgl(bspn, n, x + 0.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 ' Gauss decomposition ................................................ Sub gaudec(n, xmat(), xlumat(), iperm(), signd As Integer, rc As Integer) ' integer n size of matrix .................. ' double xmat(0:n-1,0:n-1), Input matrix .................... ' xlumat(0:n-1,0:n-1) matrix decomposition ............ ' integer iperm(0:n-1), row interchanges ................ ' signd, sign of perm .................... ' rc error code ...................... 'Note: in matrices xmat,xlumat, iperm, element [i,j] is replaced by ' vector element [n*i+j]. '*====================================================================* '* * '* 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. * '* xmat pointer to original system matrix in vector form. * '* * '* Output parameters: * '* ================== * '* xlumat pointer to LU factorization * '* iperm pointer to row permutation vector for xlumat * '* signd sign of perm. The determinant of xmat can be computed* '* as the product of the diagonal entries of xlumat * '* 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 MachEps As Double Dim d(n) As Double 'scaling vector for pivoting MachEps = 1.2E-16 If n < 1 Then rc = 1 'Invalid parameters Exit Sub End If 'copy xmat to xlumat For i = 0 To n - 1 For j = 0 To n - 1 xlumat(n * i + j) = xmat(n * i + j) Next j Next i For i = 0 To n - 1 iperm(i) = i 'Initialize iperm zmax = 0# For j = 0 To n - 1 'find row maxima tmp = Abs(xlumat(n * i + j)) If tmp > zmax Then zmax = tmp Next j If zmax = 0# Then 'xmat is singular rc = 3 Exit Sub End If d(i) = 1# / zmax Next i signd = 1 'initialize sign of iperm For i = 0 To n - 1 piv = Abs(xlumat(n * i + i)) * d(i) j0 = i 'Search for pivot element For j = i + 1 To n - 1 tmp = Abs(xlumat(n * j + 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 small, xmat is signd = 0 'nearly singular rc = 4 Exit Sub End If If j0 <> i Then signd = -signd 'update signd Call ISWAP(iperm(j0), iperm(i)) 'SWAP1 pivot entries Call SWAP1(d(j0), d(i)) 'SWAP1 scaling vector For j = 0 To n - 1 'SWAP1 j0-th and i-th rows of xlumat Call SWAP1(xlumat(n * j0 + j), xlumat(n * i + j)) Next j End If For j = i + 1 To n - 1 'Gauss elimination If xlumat(n * j + i) <> 0# Then xlumat(n * j + i) = xlumat(n * j + i) / xlumat(n * i + i) tmp = xlumat(n * j + i) For m = i + 1 To n - 1 xlumat(n * j + m) = xlumat(n * j + m) - tmp * xlumat(n * i + m) Next m End If Next j Next i rc = 0 'all ok End Sub ' Gauss solution ..................................................... Sub gausol(n, xlumat(), iperm(), b(), x(), rc As Integer) ' integer n size of matrix .................. ' real*8 xlumat(0:n-1,0:n-1) decomposed matrix (LU) .......... ' integer perm(0:n-1) row permutation vector .......... ' real*8 b(0:n-1), Right hand side ................. ' x(0:n-1) solution ........................ ' integer rc 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 xlumat. * ' xlumat Matrix(0:n-1,0:n-1) stored in a vector(n*n). * ' LU factorization, as produced from gaudec. * ' iperm integer vector(0:n-1) * ' row permutation vector for xlumat * ' b vector(0:n-1) * ' Right hand side of the system. * ' * ' Output parameter: * ' ================ * ' x vector(0:n-1) * ' Solution vector * ' * ' Return value (rc): * ' ================= * ' = 0 all ok * ' = 1 n < 1 or other invalid input parameter * ' = 3 improper LU decomposition ( zero diagonal entry) * ' * '====================================================================* Dim MachEps As Double MachEps = 1.2E-16 If n < 1 Then rc = 1 'Invalid input parameter Exit Sub End If MachEps = 1.2E-16 For k = 0 To n - 1 'update b x(k) = b(iperm(k)) For j = 0 To k - 1 x(k) = x(k) - xlumat(n * k + 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 + xlumat(n * k + j) * x(j) Next j If Abs(xlumat(n * k + k)) < MachEps Then rc = 3 Exit Sub End If x(k) = (x(k) - Sum) / xlumat(n * k + k) Next k rc = 0 'all ok End Sub ' Gauss algorithm for solving linear equations ....................... Sub gauss(mode, n, xmat(), xlumat(), iperm(), b(), x(), signd As Integer, code As Integer) ' integer mode Modus: 0, 1, 2, 3 ............... ' integer n Dimension of matrix ............. ' double xmat(0:n-1,0:n-1), Input matrix .................... ' xlumat(0:n-1,0:n-1) LU decomposition ................ ' integer iperm(0:n-1,0:n-1) row remutation vector ........... ' double b(0:n-1), right hand side ................. ' x(0:n-1) solution of the system .......... ' integer signd, sign of the permutation ......... ' code return error code ............... 'Note: in matrices xmat,xlumat, iperm, element [i,j] is replaced by ' vector element [n*i+j]. '*====================================================================* '* * '* 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 SWAP1s). * '* * '* ------------------------------------------------------------------ * '* * '* 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. * '* xmat matrix of the linear system. It is stored in vector * '* form. * '* xlumat (for mode = 2, 3) * '* LU factors of mat * '* xlumat can be stored in the space of xmat. * '* iperm (for mode = 2, 3) * '* Permutation vector, of the row interchangfes in * '* xlumat. * '* b Right hand side of the system. * '* signd (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: * '* ================== * '* xlumat (for mode = 0, 1) * '* LU factorization of xmat. * '* iperm (for mode = 0, 1) * '* row ermutation vector * '* x (for mode = 0, 2, 3) * '* solution vector(0:n-1). * '* signd (for mode = 0, 1) * '* sign of perm. * '* * '* Return value (code): * '* =================== * '* =-1 Max. number (MAXITER) of iterative refinements * '* reached (MAXITER) while mode = 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 Return End If 'Select mode If mode = 0 Then 'Find factorization and solve system ................... Call gaudec(n, xmat(), xlumat(), iperm(), signd, rc) If rc = 0 Then Call gausol(n, xlumat(), iperm(), b(), x(), rc) code = rc Else code = rc Exit Sub End If ElseIf mode = 1 Then 'Find factorization only ........................... Call gaudec(n, xmat(), xlumat(), iperm(), signd, rc) code = rc Exit Sub ElseIf mode = 2 Then 'Solve only ........................................ Call gausol(n, xlumat(), iperm(), b(), x(), rc) code = rc Exit Sub ElseIf 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 '************************************************************************ '* * '* Solve a first order system of DEs using the implicit Gear method * '* of order 4. * '* * '* Programming language: ANSI C * '* Author: Klaus Niederdrenk, 1.22.1996 (FORTRAN 77) * '* Adaptation: Juergen Dietel, Computing Center, RWTH Aachen * '* Source: FORTRAN 77 source code * '* Date: 2.26.1996 * '* * '* Quick Basic Release By J-P Moreau, Paris. * '* -------------------------------------------------------------------- * '* REF.: "Numerical Algorithms with C, by Gisela Engeln-Muellges * '* and Frank Uhlig, Springer-Verlag, 1996". * '************************************************************************ 'Gear method of 4th order for DESs of 1st order Sub gear4(x, xend, bspn As Integer, n, y(), epsabs, epsrel, h, fmax As Integer, aufrufe As Integer, fehler As Integer) ' double x, starting or end point ................ ' xend desired end point (> x) .............. ' integer bspn, # of example ......................... ' n number of DEs ........................ ' double y(0:n-1), initial value or solution ............ ' epsabs, absolute error bound ................. ' epsrel, relative error bound ................. ' h starting or final step size .......... ' integer fmax, maximal number of calls of dgl() ..... ' aufrufe, actual number of calls of dgl() ...... ' fehler error code ........................... '************************************************************************ '* Compute the value of the solution at xend, starting with the IVP. * '* We use the implicit method of Gear of 4th order with step size * '* control which is especially suited for stiff DEs. * '* The local step size control insures that the two error bounds are met* '* The number of function calls of the right hand side is limited by * '* fmax. This function can be used inside a loop to find solutions at * '* a specified point to arbitrary accuracy. * '* * '* Input parameters: * '* ================= * '* x initial value x0 * '* xend final value for the integration (xend > x0) * '* bspn # example * '* n number of DEs * '* dgl Function to compute the right hand side f(x0,y0) for the * '* system of DEs (removed from parameters). * '* y [0..n-1] solution vector y0 of the system of DEs at x0 * '* epsabs absolute error bound (>= 0); if zero, we only check for the * '* relative error. * '* epsrel relative error bound (>= 0); if zero, we only check for the * '* absolute error. * '* h given step size for first step * '* fmax maximal number of calls of the right hand side of the system* '* * '* Output parameters: * '* ================= * '* x final x-value of the integration; normally equal to xend * '* h final step size; keep for further calls * '* y [0..n-1]-vector, the solution of the system of DEs at x * '* aufrufe counter of calls of dgl() * '* * '* erroe code (fehler): * '* ==================== * '* = 0: all ok * '* = 1: Both error bounds too small * '* = 2: xend <= x0 * '* = 3: Step size h <= 0 * '* = 4: n <= 0 * '* = 5: # calls > fmax; we have not reached the desired xend; * '* repeat function call with actual values of x, y, h. * '* = 6: Jacobi matrix is singular; x, y, h are the last values * '* that could be computed with accuracy * '* = 7: ran out of memory (not used here) * '* = 8: error when calling gauss() for the second time; * '* should not occur. * '* * '************************************************************************ ' double eps1, 'MachEps^0.75; used instead of MachEps ' 'to check whether we have reached xend in ' 'order to avoid minuscule step sizes ' eps2, '100 * MachEps; for comparison with zero ' hs 'optimal step size for Jacobi matrix ' 'approximation Dim hilf(n) '(0..n-1)-vector Dim zj(5 * n) '(0:4,0:n-1)-matrix stored in a vector(5*n) Dim zjp1(5 * n) '(0:4,0:n-1)-matrix stored in a vector(5*n) Dim f(n), ykp1(n) '(0..n-1)-vectors Dim fs(n * n), fsg(n * n) '(0..n-1,0..n-1)-matrices stored in vectors(n*n) Dim con(n) '(0..n-1)-vector Dim iperm(n) '(0..n-1) permutation vector for Gauss elimination ' sg, sign of xend ' xe |xend| - eps2, carrying the sign of xend Dim amEnde As Integer 'Flag, indicating that we shall reach 'xend with the current step ' ymax, Maximum norm of the current 'approximate value of y ' dummy, aux storage for h ' auxiliary variables (double); ' xka, xke , hka, hk1, diff, eps, q, halt, quot1, quot2, quot3, quot4 Dim done As Integer 'Boolean (0 or 1) ' nochmal 'Boolean (0 or 1) Dim aufrufeawp As Integer Dim signdet As Integer 'sign of determinant in Gauss Dim dum(n), dum1(n) 'dummy vectors for gauss Dim MachEps As Double 'machine smallest real number 'auxiliary vectors 'Dim yhilf(n), k1(n), k2(n), k3(n), k4(n), k5(n), k6(n) ' ------------------------- Initialize ------------------------------ dummy = 0#: done = 0 MachEps = 1.2E-16 'for IBM PC eps1 = MachEps ^ 0.75 eps2 = 100# * MachEps hs = 10# * Sqr(MachEps) ONE = 1# If (xend >= 0#) Then sg = 1# Else sg = -1# End If xe = (1# - sg * eps2) * xend fehler = 0 aufrufe = 1 amEnde = 0 ' --------- check input parameters ------------------- ymax = XNormMax(y(), n) If epsabs <= eps2 * ymax And epsrel <= eps2 Then fehler = 1 ElseIf xe < x Then fehler = 2 ElseIf h < eps2 * Abs(x) Then fehler = 3 ElseIf n <= 0 Then fehler = 4 End If If fehler > 0 Then Exit Sub ' ------------ first integration step --------------- If x + h > xe Then h = xend - x dummy = h amEnde = 1 End If For i = 0 To n - 1 hilf(i) = y(i) Next i xka = x xke = xka hka = 0.25 * h hk1 = hka For k = 1 To 4 xke = xke + hka awp xka, xke, bspn, n, hilf(), epsabs, epsrel, hk1, 6, fmax - aufrufe, aufrufeawp, fehler aufrufe = aufrufe + aufrufeawp If fehler <> 0 Then Exit Sub For i = 0 To n - 1 zjp1(n * k + i) = hilf(i) Next i Next k dgl bspn, n, x, y(), f() aufrufe = aufrufe + 1 ' ---------- Compute first Gear-Nordsieck approximation ------------------- For i = 0 To n - 1 zj(i) = y(i) zj(i + n) = h * f(i) zj(i + 2 * n) = ONE / 24 * (35 * y(i) - 104 * zjp1(i + n) + 114 * zjp1(i + 2 * n) - 56 * zjp1(i + 3 * n) + 11 * zjp1(i + 4 * n)) zj(i + 3 * n) = ONE / 12 * (-5 * y(i) + 18 * zjp1(i + n) - 24 * zjp1(i + 2 * n) + 14 * zjp1(i + 3 * n) - 3 * zjp1(i + 4 * n)) zj(i + 4 * n) = ONE / 24 * (y(i) - 4 * zjp1(i + n) + 6 * zjp1(i + 2 * n) - 4 * zjp1(i + 3 * n) + zjp1(i + 4 * n)) Next i ' ------------------------ adjust step size -------------------------- While done = 0 ' --- use Newton method for an implicit approximation --- For i = 0 To n - 1 ykp1(i) = zj(i) + zj(i + n) + zj(i + 2 * n) + zj(i + 3 * n) + zj(i + 4 * n) Next i dgl bspn, n, x + h, ykp1(), f() For k = 0 To n - 1 'copy vector ykp1 in hilf For i = 0 To n - 1 hilf(i) = ykp1(i) Next i hilf(k) = hilf(k) - hs dgl bspn, n, x + h, hilf(), dum() For i = 0 To n - 1 fs(k * n + i) = dum(i) Next i For i = 0 To n - 1 fs(k * n + i) = -h * 0.48 * (f(i) - fs(k * n + i)) / hs Next i fs(k * n + k) = fs(k * n + k) + ONE Next k 'update number of calls to dgl aufrufe = aufrufe + n + 1 For i = 0 To n - 1 con(i) = ykp1(i) - 0.48 * (zj(i + n) + 2 * zj(i + 2 * n) + 3 * zj(i + 3 * n) + 4 * zj(i + 4 * n)) For k = 0 To n - 1 fsg(k * n + i) = fs(i * n + k) Next k Next i gauss 1, n, fsg(), fsg(), iperm(), dum(), dum1(), signdet, fehler If fehler > 0 Then 'error in gauss ? fehler = 6 Exit Sub End If For iter = 1 To 3 For i = 0 To n - 1 hilf(i) = -ykp1(i) For k = 0 To n - 1 hilf(i) = hilf(i) + fs(k * n + i) * ykp1(k) Next k hilf(i) = h * 0.48 * f(i) + hilf(i) + con(i) Next i gauss 2, n, fsg(), fsg(), iperm(), hilf(), ykp1(), signdet, fehler If fehler > 0 Then fehler = 8 Exit Sub End If dgl bspn, n, x + h, ykp1(), f() Next iter 'update number of calls to dgl aufrufe = aufrufe + 3 ' ---- Compute corresponding Gear-Nordsieck approximation ---- For i = 0 To n - 1 hilf(i) = h * f(i) - zj(i + n) - 2 * zj(i + 2 * n) - 3 * zj(i + 3 * n) - 4 * zj(i + 4 * n) Next i For i = 0 To n - 1 zjp1(i) = ykp1(i) zjp1(i + n) = h * f(i) zjp1(i + 2 * n) = zj(i + 2 * n) + 3# * zj(i + 3 * n) + 6# * zj(i + 4 * n) + 0.7 * hilf(i) zjp1(i + 3 * n) = zj(i + 3 * n) + 4# * zj(i + 4 * n) + 0.2 * hilf(i) zjp1(i + 4 * n) = zj(i + 4 * n) + 0.02 * hilf(i) Next i ' --- decide whether to accept last step --- ' copy vector zjp1(4) in hilf and zj(4) in con For i = 0 To n - 1 hilf(i) = zjp1(i + 4 * n) con(i) = zj(i + 4 * n) Next i diff = DistMax(hilf(), con(), n) ymax = XNormMax(ykp1(), n) eps = (epsabs + epsrel * ymax) / 6# q = Sqr(Sqr(eps / diff)) / 1.2 If (diff < eps) Then ' --- accept last step; prepare for next one --- x = x + h 'copy vector ykp1 in y For i = 0 To n - 1 y(i) = ykp1(i) Next i ' stop integration, if interval end xend has been reached ' or if there has been too many function dgl calls. nochmal = 0 While nochmal = 0 If amEnde <> 0 Then h = dummy Exit Sub ElseIf aufrufe > fmax Then fehler = 5 Exit Sub End If ' --- adjust step size for next step --- halt = h h = XMin(q, 2#) * h If x + h >= xe Then dummy = h h = xend - x amEnde = 1 ' --- close enough to xend => stop integration --- If h < eps1 * Abs(xend) Then nochmal = 1 End If If nochmal = 0 Then GoTo 10 Wend ' ------ compute Gear-Nordsieck approximation ----- ' ------ for the next step ----- 10 quot1 = h / halt quot2 = quot1 ^ 2 quot3 = quot2 * quot1 quot4 = quot3 * quot1 For i = 0 To n - 1 zj(i) = zjp1(i) zj(i + n) = quot1 * zjp1(i + n) zj(i + 2 * n) = quot2 * zjp1(i + 2 * n) zj(i + 3 * n) = quot3 * zjp1(i + 3 * n) zj(i + 4 * n) = quot4 * zjp1(i + 4 * n) Next i Else ' ------ repeat last step with smaller step size; ----- ' -------- adjust Gear-Nordsieck approximation --------- halt = h h = XMax(0.5, q) * h quot1 = h / halt quot2 = quot1 ^ 2 quot3 = quot2 * quot1 quot4 = quot3 * quot1 For i = 0 To n - 1 zj(i + n) = quot1 * zj(i + n) zj(i + 2 * n) = quot2 * zj(i + 2 * n) zj(i + 3 * n) = quot3 * zj(i + 3 * n) zj(i + 4 * n) = quot4 * zj(i + 4 * n) Next i amEnde = 0 End If Wend 'while done=0 End Sub 'gear4 Sub ISWAP(ia, ib) Dim tnp As Integer tmp = ib: ib = ia: ia = tmp End Sub ' Gauss for multiple right hand sides ................................ Sub mgauss(n, k, xmat(), rmat(), code As Integer) ' integer n, Dimension of system ............. ' k number of right hand sides ...... ' real*8 mat(0:n-1,0:n-1), original matrix ................. ' rmat(0:n-1,0:n-1) Right hand sides/solutions ...... ' integer code Error code ...................... '*====================================================================* '* * '* mgauss finds the solution matrix x for the linear system * '* mat * x = rmat with an n x n coefficient matrix mat and a * '* n x k matrix rmat of right hand sides. Here mat must be * '* nonsingular. * '* * '* ------------------------------------------------------------------ * '* * '* Input parameters: * '* ================ * '* n integer; (n > 0) * '* Dimension of xmat. * '* k integer k; (k > 0) * '* number of right hand sides * '* xmat matrix(0:n-1,0:n-1) stored in a vector(n*n) * '* n x n original system matrix * '* rmat matrix(0:n-1,0:n-1) stored in a vector(n*n) * '* Right hand sides * '* * '* Output parameter: * '* ================ * '* rmat solution matrix for the system. * '* The input right hand sides are lost. * '* * '* Return value (code): * '* =================== * '* = 0 all ok * '* = 1 n < 1 or k < 1 or invalid input parameter * '* = 2 lack of memory (not used here) * '* = 3 mat is numerically singular. * '* * '* ------------------------------------------------------------------ * '* * '* subroutine used: * '* ================ * '* * '* gaudec: LU decomposition of mat. * '* * '*====================================================================* Dim signd As Integer, rc As Integer Dim iperm(n) Dim xlu(n * n), x(n) Dim MachEps As Double If n < 1 Or k < 1 Then 'Invalid parameter code = 1 Exit Sub End If MachEps = 1.2E-16 Call gaudec(n, xmat(), xlu(), iperm(), signd, rc) 'compute factorization 'in matrix lu If rc <> 0 Or signd = 0 Then 'if not possible code = 3 'exit with code=3 Exit Sub End If For m = 0 To k - 1 'Loop over the right hand sides For i = 0 To n - 1 'Updating the b's x(i) = rmat(n * iperm(i) + m) For j = 0 To i - 1 x(i) = x(i) - xlu(n * i + j) * x(j) Next j Next i For i = n - 1 To 0 Step -1 'back substitution Sum = 0# For j = i + 1 To n - 1 Sum = Sum + xlu(n * i + j) * x(j) Next j If Abs(xlu(n * i + i)) < MachEps Then 'invalid LU decomposition code = 2 Exit Sub End If x(i) = (x(i) - Sum) / xlu(n * i + i) Next i For j = 0 To n - 1 'Save result rmat(n * j + m) = x(j) Next j Next m code = 0 End Sub ' 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) + 0.2 * h * k1(i) Next i Call dgl(bspn, n, x + 0.2 * h, yhilf(), k2()) For i = 0 To n - 1 yhilf(i) = y(i) + 0.075 * h * (k1(i) + 3# * k2(i)) Next i Call dgl(bspn, n, x + 0.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 + 0.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) + 0.13 * k4(i) + 0.144 * k5(i)) g7(i) = h * (2.134 * k1(i) + 0.24 * k3(i) + 0.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 '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 ' 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) + 0.25 * h * (k1(i) + k2(i)) Next i Call dgl(bspn, n, x + 0.5 * h, yhilf(), k3()) For i = 0 To n - 1 y2(i) = y(i) + 0.5 * h * (k1(i) + k2(i)) y3(i) = y(i) + h / 6# * (k1(i) + k2(i) + 4# * k3(i)) Next i End Sub '*************************************************************** '* This subroutine defines n, number of DEs and dgltxt(), the * '* text decription of current example. * '*************************************************************** 'NOTE: here, only examples #0, #3 are implemented. Sub settxt(bspn As Integer, n, dgltxt$()) If bspn = 0 Then n = 2 dgltxt$(0) = " y1' = y1 * y2 + cos(x) - 0.5 * sin(2.0*x)" dgltxt$(1) = " y2' = y1 * y1 + y2 * y2 - (1 + sin(x))" ElseIf bspn = 3 Then n = 5 dgltxt$(0) = " y1' = y2" dgltxt$(1) = " y2' = y3" dgltxt$(2) = " y3' = y4" dgltxt$(3) = " y4' = y5" dgltxt$(4) = " y5' = (45 * y3 * y4 * y5 - 40 * y4 * y4 * y4) / (9 * y3 * y3)" End If End Sub Sub SWAP1(a, b) tmp = b: b = a: a = tmp 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 ' -------------------------- END mgear1.bas -------------------------
No comments:
Post a Comment