'******************************************************************** '* Solve Y' = F(X,Y) with Initial Condition Y(X0)=Y0 using the * '* Euler-Romberg Method * '* ---------------------------------------------------------------- * '* SAMPLE RUN: * '* (Solve Y' = X*X*Y with Y(0)=1 for X0=0 up to Xn=1.1, exact * '* solution is Y = Exp(X^3/3) ). * '* * '* X Y Y Error Number of * '* estimated exact subdivisions * '* ---------------------------------------------------------- * '* 0.1 1.000333 1.000333 0.00000000 4 * '* 0.2 1.002670 1.002670 0.00000001 4 * '* 0.3 1.009041 1.009041 0.00000006 4 * '* 0.4 1.021562 1.021562 0.00000014 4 * '* 0.5 1.042547 1.042547 0.00000027 4 * '* 0.6 1.074654 1.074655 0.00000086 4 * '* 0.7 1.121125 1.121126 0.00000107 4 * '* 0.8 1.186094 1.186095 0.00000126 4 * '* 0.9 1.275067 1.275069 0.00000133 4 * '* 1.0 1.395611 1.395612 0.00000114 4 * '* 1.1 1.558410 1.558412 0.00000047 4 * '* ---------------------------------------------------------- * '* * '* Ref.: "Methodes de calcul numerique By Claude Nowakowski, Tome 2 * '* PSI Editions, France, 1981" [BIBLI 04]. * '* * '******************************************************************** 'Program EulerRomberg DEFINT I-N DEFDBL A-H, O-Z OPTION BASE 0 NMAX = 100 H = .1 'initial integration step ER = .000001# 'desired precision LA = 10 'maximum number of subdivisions NC = 10 'number of calculations NC = (Xn+1 - X1)/H DIM X(NMAX), Y(NMAX), T(20) 'Initial conditions X(0) = 0#: Y(0) = 1# 'write header CLS PRINT PRINT " X Y Y Error Number of " PRINT " estimated exact subdivisions " PRINT "---------------------------------------------------------" 'main integration loop FOR N = 0 TO NC XC = X(N): YC = Y(N) XX = XC: YY = YC: GOSUB 1000 T(1) = Y(N) + H * F L = 1: LM = 2: ET=1# WHILE L < LA AND ET >= ER XC = X(N): YC = Y(N) FOR J = 1 TO LM XC = XC + H / LM XX = XC: YY = YC: GOSUB 1000 YC = YC + H / LM * F NEXT J T(L + 1) = YC: M = 1: K = L: MM = 2: ET = 1# IF K > 1 THEN WHILE ET >= ER AND K > 1 T(K) = (MM * T(K + 1) - T(K)) / (MM - 1) ET = ABS(T(K) - T(K - 1)) M = M + 1: K = K - 1: MM = MM * 2 WEND END IF IF K = 1 THEN L = L + 1: LM = LM * 2 END IF WEND X(N + 1) = X(N) + H: Y(N + 1) = T(K) XX = X(N + 1): GOSUB 2000: YEX = FX EF = ABS(YEX - Y(N + 1)) PRINT USING "###.#"; X(N + 1); PRINT USING "#####.######"; Y(N + 1); PRINT USING "#####.######"; YEX; PRINT USING "#####.########"; EF; PRINT " "; L NEXT N PRINT "---------------------------------------------------------" END 'of main program 'Y' = F(X,Y) 1000 'Function F(XX,YY) F = XX * XX * YY RETURN 'Exact solution FX(XX) 2000 'Function FX(XX) FX = EXP(XX ^ 3 / 3) RETURN 'end of file eulromb.bas