Attribute VB_Name = "Module3" '************************************************************************ '* * '* 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 * '* * '* Visual Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '* -------------------------------------------------------------------- * '* REF.: "Numerical Algorithms with C, By Gisela Engeln-Muellges * '* and Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11]. * '************************************************************************ DefDbl A-H, O-Z DefInt I-N Option Base 0 Dim hilf() Dim zj(), zjp1() Dim f(), ykp1() Dim con(), iperm() 'print a REAL square matrix with name (debug only) Sub PrintMat(Name As String, n, xmat()) ' Name Matrix caption ' n Size of matrix ' double xmat(0:n-1,0:n-1) stored in a vector(n*n) ' matrix to be printed Form1.Print " "; Name For i = 0 To n - 1 For j = 0 To n - 1 Form1.Print " "; xmat(n * i + j); Next j Form1.Print Next i End Sub '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, 'MACH_EPS^0.75; used instead of MACH_EPS ' 'to check whether we have reached xend in ' 'order to avoid minuscule step sizes ' eps2, '100 * MACH_EPS; for comparison with zero ' hs 'optimal step size for Jacobi matrix ' 'approximation ReDim hilf(n) '(0..n-1)-vector ReDim zj(5 * n) '(0:4,0:n-1)-matrix stored in a vector(5*n) ReDim zjp1(5 * n) '(0:4,0:n-1)-matrix stored in a vector(5*n) ReDim f(n), ykp1(n) '(0..n-1)-vectors ReDim fs(n * n), fsg(n * n) '(0..n-1,0..n-1)-matrices stored in vectors(n*n) ReDim con(n) '(0..n-1)-vector ReDim 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 aufrufe_awp As Integer Dim signdet As Integer 'sign of determinant in Gauss ReDim dum(n), dum1(n) 'dummy vectors for gauss Dim MACH_EPS 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 MACH_EPS = 1.2E-16 'for IBM PC eps1 = MACH_EPS ^ 0.75 eps2 = 100# * MACH_EPS hs = 10# * Sqr(MACH_EPS) 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 = norm_max(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, aufrufe_awp, fehler aufrufe = aufrufe + aufrufe_awp 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 = dist_max(hilf, con, n) ymax = norm_max(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 ' -------------------------- END gear.bas -------------------------