'******************************************************************** '* 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 = 0.1 'initial integration step ER = 0.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

## Thursday, September 7, 2017

### Euler-Romberg Method

Subscribe to:
Post Comments (Atom)

## No comments:

## Post a Comment