Attribute VB_Name = "Module4" DefDbl A-H, O-Z DefInt I-N Option Base 0 '************************************************************************ '* 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]. * '* * '* Visual Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '************************************************************************ ' Add to project: awp.bas, fgauss.bas, gear.bas, t_dgls.bas ' --------------------------------------------------------- Sub Exec_Mgear() ' 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(4) As String 'text of equations ' -------------------- read input -------------------- Form1.Print ' Input example number (0 to 11): ' bspnummer = 0 for example 0 bspnummer = 3 If bspnummer <> 0 And bspnummer <> 3 Then Form1.Print Form1.Print " Example not registered." Exit Sub End If ' input absolute and relative errors ' epsabs = 0.000001: epsrel = 0.00000001 (Ex. 0) epsabs = 0.0000000001: epsrel = 0.0000000001 ' input x starting value x0 = 0# Call settxt(bspnummer, n, dgltxt) 'read text of equations 'and parameter n ' allocate memory for vectors y0, yex ReDim y0(n) ReDim yex(n) ' input initial values y0(i) ' y0(0) = 0.5: y0(1) = 0.5 (Ex. 0) For i = 0 To n - 1 y0(i) = 1# Next i ' input initial step size ' h = 0.0001 (Ex. 0) h = 0.001 ' input ending x value xend = 1.5 ' input maximum number of calls to right hand side fmax = 6000 ' ----------------- print input data ---------------------- Form1.Print " ==============================================" Form1.Print " Solve a first order ordinary system of DEs " Form1.Print " using the implicit method of Gear of 4th order" Form1.Print " ==============================================" Form1.Print Form1.Print " System of DEs:" Form1.Print " ------------- " For i = 0 To n - 1 Form1.Print " "; dgltxt(i) Next i Form1.Print " Example # "; bspnummer Form1.Print " Number of DEs = "; n Form1.Print " x0 = "; x0 Form1.Print " xend = "; xend Form1.Print " epsabs = "; epsabs Form1.Print " epsrel = "; epsrel Form1.Print " fmax = "; fmax Form1.Print " h = "; h For i = 0 To n - 1 Form1.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 Form1.Print " Gear4: error n° "; 10 + fehler Exit Sub End If ' ---------------------- print results --------------------- Form1.Print Form1.Print " Output data:" Form1.Print " -----------" Form1.Print " error code from gear4: "; fehler Form1.Print " final local step size: "; h Form1.Print " number of calls of right hand side: "; aufrufe Form1.Print " Integration stopped at x = "; x0 Form1.Print For i = 0 To n - 1 Form1.Print " approximate solution y"; i + 1; "(x) = "; y0(i) Next i Form1.Print End Sub ' -------------------------- END mgear.bas -------------------------