Euler-Romberg Method

'********************************************************************
'*    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