Monday, September 11, 2017

Adams-Bashforth Method

'********************************************************************
'*    Solve Y' = F(X,Y) with Initial Condition Y(X0)=Y0 using the   *
'*    Adams-Bashforth Method                                        *
'* ---------------------------------------------------------------- *
'* REFERENCE: "Méthode de calcul numérique- Tome 2 - Programmes en  *
'*             Basic et en Pascal By Claude Nowakowski, Edition du  *
'*             P.S.I., 1984" [4].                                   *
'* ---------------------------------------------------------------- *
'* SAMPLE RUN:                                                      *
'* (Solve Y' = -Y + X/((1+X)*(1+X))  with Y(0)=1 for X0=0 up to     *
'*  X1=1.0, exact solution is Y = 1 / (1+X).                        *
'*                                                                  *
'*      X           Y        Y exact      Error                     *
'* -----------------------------------------------                  *
'*   0.000000    1.000000    1.000000   0                           *
'*   0.000000    1.000000    1.000000   0.000000                    *
'*   0.050000    0.952381    0.952381   0.000000                    *
'*   0.100000    0.909091    0.909091   0.000000                    *
'*   0.150000    0.869525    0.869565   0.000040                    *
'*   0.200000    0.833265    0.833333   0.000068                    *
'*   0.250000    0.799910    0.800000   0.000090                    *
'*   0.300000    0.769125    0.769231   0.000106                    *
'*   0.350000    0.740623    0.740741   0.000117                    *
'*   0.400000    0.714160    0.714286   0.000125                    *
'*   0.450000    0.689525    0.689655   0.000131                    *
'*   0.500000    0.666533    0.666667   0.000134                    *
'*   0.550000    0.645026    0.645161   0.000135                    *
'*   0.600000    0.624865    0.625000   0.000135                    *
'*   0.650000    0.605926    0.606061   0.000134                    *
'*   0.700000    0.588103    0.588235   0.000133                    *
'*   0.750000    0.571298    0.571429   0.000131                    *
'*   0.800000    0.555428    0.555556   0.000128                    *
'*   0.850000    0.540416    0.540541   0.000125                    *
'*   0.900000    0.526194    0.526316   0.000121                    *
'*   0.950000    0.512703    0.512821   0.000118                    *
'*   1.000000    0.499886    0.500000   0.000114                    *
'* -----------------------------------------------                  *
'*                                                                  *
'*                              Basic Release By J-P Moreau, Paris. *
'*                                       (www.jpmoreau.fr)          *
'********************************************************************
DefDbl A-H, O-Z
DefInt I-N

Option Base 0

  H = 0.05 'integration step

  Dim B(4), X(4), Y(4)

  'Initial conditions
  X(0) = 0#      'starting X
  X1 = 1#        'ending X
  Y(0) = 1#      'initial Y
  
  'write header
  Cls
  Print
  Print "       X           Y        Y exact      Error    "
  Print "  ------------------------------------------------"
  'write initial line
  F$ = "#####.######"
  Print USING; F$; X(0);
  Print USING; F$; Y(0);
  Print USING; F$; Y(0);
  Print "    0"
  'use Runge-Kutta to start
  For K = 0 To 1
    GoSub 1000 'call RK4
    XX = X(K + 1): GoSub 600: ER = Abs(FX - Y(K + 1))
    Print USING; F$; X(K + 1);
    Print USING; F$; Y(K + 1);
    Print USING; F$; FX;
    Print USING; F$; ER
  Next K
  'main integration loop
  While X(2) < X1
    For I = 1 To 3
      XX = X(3 - I): YY = Y(3 - I): GoSub 500
      B(I) = F
    Next I
    X(3) = X(2) + H
    Y(3) = Y(2) + H * (23# * B(1) - 16# * B(2) + 5# * B(3)) / 12#
    XX = X(3): GoSub 600: ER = Abs(Y(3) - FX)
    Print USING; F$; X(3);
    Print USING; F$; Y(3);
    Print USING; F$; FX;
    Print USING; F$; ER
    For K = 0 To 2
      X(K) = X(K + 1): Y(K) = Y(K + 1)
    Next K
  Wend
  Print "  ------------------------------------------------"

End

'User defined function Y'=F(XX,YY)
500  F = -YY + XX / ((1# + XX) * (1# + XX))
     Return

'Exact solution Y=FX(XX)
600  FX = 1# / (1# + XX)
     Return

'Runge-Kutta method to calculate first points only
1000 XX = X(K): YY = Y(K): GoSub 500: C1 = F
     XX = X(K) + H / 2#: YY = Y(K) + H / 2# * C1: GoSub 500: C2 = F
     XX = X(K) + H / 2#: YY = Y(K) + H / 2# * C2: GoSub 500: C3 = F
     XX = X(K) + H: YY = Y(K) + H * C3: GoSub 500: C4 = F
     X(K + 1) = X(K) + H
     Y(K + 1) = Y(K) + H * (C1 + 2 * C2 + 2 * C3 + C4) / 6#
     Return

'end of file adambash.bas



 EXPLANATION FILE OF PROGRAM ADAMBASH
====================================


  Ordinary Differential Equations Y' = F(x,y)
  -------------------------------------------


  Linked Steps Method
  -------------------

    We have seen that Euler, Runge-Kutta... methods can be put under the general
  form:

                y    = y  + h phi(x , y , h)
                 n+1    n          n   n

  and each point y    of the solution is only determined from previous point, y
                  n+1                                                          n

  Calculations are made for y1, y2,...,yn-1, yn, yn+1... These methods are cal-
  led "with separate steps": each step is independant from the previous one. To
  obtain a given accuracy, each step has to be divided into intermediate steps.

  Let us take an example:

  For the Runge-Kutta method of order 4, we use the formulas:

                k1 = h f(xn, yn)

                k2 = h f(xn+h/2, yn+k1/2)

                k3 = h f(xn+h/2, yn+k2/2)

                k4 = h f(xn+h, yn+k3)

                y    = y  + (1/6)(k1+2k2+2k3+k4)
                 n+1    n

  Another method consists in using the previous calculations to improve speed;
  the y   point is evaluated from y , y    and y    points.
       n+1                         n   n-1      n-2

  The general formula is:

                y    = a y + a y   + ... + a y 
                 n+1    0 n   1 n-1         k n-k

                       + h (b  f   + b f + b f    + ... + b f
                             -1 n+1   0 n   1 n-1          k n-k

  This can be written as:

                        k              k
                y    = Sum a y    + h Sum  b  f
                 n+1   j=0  j n-j     j=-1  j  n-j


    These methods are called "with linked steps" and f(x,y) evaluations are only
  made at points x0, x1, x2,..., xn.

  If b   = 0, the process is explicit; y    is directly obtained by applying the
      -1                                n+1
  formula.

  If b   <> 0, the process is implicit: we must solve an equation of the form
      -1

  y    = phi(y   ) to obtain y   .
   n+1        n+1             n+1

    For the three first points, y1, y2, y3, we have no previous points to calcu-
  late them: these methods cannot start by themselves like "with separate steps"
  methods.

    The truncation error Tn is estimated by using the Taylor formula for y(x   )
                                                                            n+j
                                 (jh)²                (jh)^p  (p)
    y(x   ) = y(x ) + jhy'(x ) + ----- y"(x ) + ... + ------ y   (x )
       n+j       n          n      2!      n            p!         n

                    + h^p eps(h)

    and also for f(x   , y(x   )):
                    n+j     n+j


    f(x   , y(x   )) = y'(x +jh)
       n+j     n+j         n                         p-1
                                 h               (jh)      (p)
                     = y'(x ) + j y"(x ) + ... + -------- y   (x ) + h^p eps(h)
                           n          n           (p-1)!        n

    Example: let us estimate the truncation error for the implicit formula:

                     y    = y  + (h/2) (f    + f )
                      n+1    n           n+1    n

                     T    = y(x   ) - y
                      n+1      n+1     n+1
                                                                    3
    y(x   ) = y(x ) + hy'(x ) + (h^2/2) y"(x ) + (h^3/3) y"'(x ) + h  eps(h)
       n+1       n         n                n                 n

                                                    2               3
    (h/2) f    = (h/2) y'(x + h) = (h/2) y'(x ) + (h /2) y"(x ) + (h /4)
           n+1             n                 n               n
                              3
                   y"'(x ) + h  eps(h)
                        n

    (h/2) f  = (h/2) y'(x )
           n             n
                                                                  
    Hence   y(x   ) - y    = y(x   ) - y  + (h/2) (f    + f ) = 
               n+1     n+1      n+1     n           n+1    n

                               3                 3
                           = (h /12) y"'(x )  + h  eps(h)
                                          n

    Here the method is of order two.


    The methods with separate steps can be integrated by the formula

                                    xn+2h
                y(x   ) - y(x   ) =  Sum  f(t,y(t)) dt
                   n+1       n-M     xn

    and by applying the Simpson's formula

                 b             b - a            a + h
                Sum (f(t) dt = ----- {f(x) + 4f(-----) + f(b)}
                 a               6                2

    Hence       y    - y  = (h/3)[f + 4f   + f   ]
                 n+2    n          n    n+1   n+2

    Here we have an implicit process.


    In a more general way, knowing y , y   , y   , we can calculate f , f   ,
                                    n   n-1   n-2                    n   n-1

  f    and approximate y' = f(x,y) by an interpolation polynomial at points  
   n-2

  x , x   , x   :
   n   n-1   n-2                      x
                                       n+1
                        y    = y    + Sum  P(x) dx 
                         n+1    n-M   x
                                       n-M

   where y    is an approximation of y(x   ) and y    is an approximation of
          n-M                           n-M       n+1

   y(x   ).
      n+1

  In the case of an implicit method, y(x   ) can be approximated by the formula:
                                        n+1
                                     x
                         p            n+1
                        y   = y    + Sum  P(x) dx
                         n+1   n-M   x
                                      n-M

                          p             p
  This allows evaluating f   = f(x   , y   ) and we can resume the interpolation
                          n+1     n+1   n+1

  step with a new polynomial P*(x). y    is calculated by the correction formula
                                     n+1

                                     x
                         c            n+1
                        y   = y    + Sum  P*(x) dx
                         n+1   n-M   x
                                      n-M

  As the points are equally spaced and indices are decreasing, xn, xn-1, xn-2...
  we can use the Newton formula with back diferences:

             Div (f ) = f - f   
                   n     n   n-1

                2
             Div (f ) = f - 2f    +f
                   n     n    n-1   n-2

             --------------------------
                k           |k|      |k|               k
             Div (f ) = f - | | f    | | f   +...+ (-1) f
                   n     n  |1|  n-1 |2|  n-2            n-k

  Hence              x                             2
         p            n+1      Div(fn)          Div (fn)
        y   = y    + Sum  [f + ------- (x-x ) + -------- (x-x )(x-x   ) + ...
         n+1   n-M   x      n    h         n      2h^2       n     n-1
                      n-M

                        k
                     Div (fn)
                   + -------- (x-x )(x-x   )...(x-x     )] dx
                      k! h^k      n     n-1        n-k+1

        Let us put u = (x-x )/h, then  du = dx/dh and
                           n                            2
                    p             1                  Div (fn)
                   y    = y    + Sum [f + Div(f )u + -------- u(u+1) + ...
                    n+1    n-M   -M    n       n        2

                        k
                     Div (fn)
                   + -------- u(u+1)(u+2)...(u+k-1) h du
                        k!

       After integration:

        p               _     _           _     2             _    k
       y    = y    + h [P f + P Div(f ) + P  Div (f ) + ... + P Div (f )]
        n+1    n-M       0 n   1     n     2       n           k      n

              _    1   1
       where  P  = -- Sum u(u+1)(u+2)...(u+j-1) du 
               j   j! -M

                j           |j|      |j|               j
       and   Div (f ) = f - | | f    | | f   +...+ (-1) f
                   n     n  |1|  n-1 |2|  n-2            n-j

                            |j|      j!
                    with    | | = ---------   (Newton's coefficients) 
                            |m|   m! (j-m)!

       Fianally   p
                 y    = y
                  n+1    n-M + h [P f  + P f    + ... + P f   ]
                                   0 n    1 n-1          k n-k

       Example: M=3, k=2

                     x                              2
        p             n+1       Div(fn)          Div (fn)
       y    = y    + Sum  [f  + ------- (x-x ) + -------- (x-x )(x-x   )] dx
        n+1    n-3   x      n     n         n      2h²        n     n-1
                      n-3

       Let us put u = (x-xn)/h, then
                                             2
        p             1                   Div (fn)
       y    = y    + Sum [f + Div(f ) u + -------- u (u+1)] h du
        n+1    n-3   -3    n       n         2

       After integration:

        p                                        2
       y    = y    + h [4f - 4 Div(f ) + (8/3)Div (f )
        n+1    n-3        n         n               n

  This leads to the Milne's formula (explicit process):

        p
       y    = y    + (4h/3) [2f - f   + 2f   ]
        n+1    n-3             n   n-1    n-2

  We do in a similar way for the implicit process. For example, with M=1 and
  k=2, the Milne's corrector formula is:

        c
       y    = y    + (h/3) [f    + 4f  + f   ]
        n+1    n-1           n+1     n    n-1

  This corresponds to the numerical integration of

       x
        n+1
       Sum  f(x) dx   by the Simpson's method. 
       x
        n-1

  For the Adams-Moulton's implicit formulas, we have:

              p          h
       k=1   y    = y  + - [3f - f   ]
              n+1    n   2    n   n-1

              c          h
             y    = y  + - [f   + f ]
              n+1    n   2   n+1   n

              p          h
       k=2   y    = y  + -- [23f - 16f   + 5f   ]
              n+1    n   12     n     n-1    n-2

              c          h
             y    = y  + -- [5f   + 8f  - f   ]
              n+1    n   12    n+1    n    n-1

              p          h
       k=3   y    = y  + -- [55f - 59f   + 17f    - 9f   ]
              n+1    n   24     n     n-1     n-2     n-3

              c          h
             y    = y  + -- [9f   + 19f  - 5f    +f   ]
              n+1    n   24    n+1     n     n-1   n-2

              p           h
       k=4   y    = y  + --- [1901f - 2984f   + 2616f    - 1274f   + 251f   ]
              n+1    n   720       n       n-1       n-2        n-3      n-4

              c           h
             y    = y  + --- [251f   + 646f  - 264f    +106f   - 19f   ]
              n+1    n   720      n+1      n       n-1      n-2     n-3


    In these formulas, the main part of what is left aside in the integration
  corresponds to the truncation error.
                                                 3
                  c                             h   "'
    Example: for y    = y  + (h/2) (f   + f ) - -- y  (ksi)
                  n+1    n           n+1   n    12

             with x  <= ksi <= x
                   n            n+1

    Another way to obtain an explicit or implicit process is using a recursive
  formula:

            y    = a y  + a y    + ... + a y
             n+1    0 n    1 n-1          k n-k

                                       + h [b  f    + b f  + ... + b y   ]
                                             -1 n+1    0 n          k n-k

    For y, the polynomial of greatest degree is:

           y(x) = 1,  y'(x) = 0
           y(x) = x,  y'(x) = 1
           y(x) = x², y'(x) = 2x
           ---------------------
                   m            m-1
           y(x) = x , y'(x) = mx

    So
           1 = a0 + a1 + ... + ak
           x = (x-h)a  + (x-2h)a  + ... + ha   + b  +b + b + ... +b
                     0          1           k-1   -1  0   1        k
           
           x² = (x-h)²a  + (x-2h)²a  + ... + ha
                       0           1           k-1

                       + 2 [xb  +(x-h)b + ... + hb
                              -1       0          k-1
           --------------------------------------------------
            l        l           l
           x  = (x-h) a  + (x-2h) a  + ... + ha
                       0           1           k-1

                             l-1          l-1
                       + l [x   b   +(x-h)   b + ... + hb
                                 -1           0          k-1

  Example: Nystroem's explicit formula:

                                        3
                       y    = y    + h Sum b  f
                        n+1    n-1     j=0  j  n-j

                            3         2
                       y = x , y' = 3x

                      3                3
         ==>  y    = x ,  y    = (x-2h)
               n+1         n-1
                                        2                2                2
                       f  = y'  = 3(x-h) , f    = 3(x-2h) , f    = 3(x-3h)
                        n     n             n-1              n-2

               3         2              2           2           2
         ==>  x  = (x-2h)  + 3h [b (x-h)  + b (x-2h)  + b (x-3h)
                                  0          1           2

  By developing and identifying:

                   b0 + b1 + b2 = 2

                   b0 + 2b1 +3b2 = 2

                   b0 + 4b1 + 9b2 = 8/3

     We find:      b0 = 7, b1 = -2, b2 = 1.


  So we have the Adams's formulas:

     Explicit, order 2:  y    = y  + h/2 (3f  - f   )
                          n+1    n          n    n-1

               order 3:  y    = y  + h/12 (23f  - 16f    + 5f   )
                          n+1    n            n      n-1     n-2

     Implicit, order 2:  y    = y  + h/2 (f   + f )
                          n+1    n         n+1   n

  and Nystroem's formulas:

     Explicit, order 2:  y    = y   + 2h f
                          n+1    n-1      n

               order 3:  y    = y   + h/3 (7f  - 2f    + f   )
                          n+1    n+1         n     n-1    n-2

     Implicit, order 3:  y    = y   + h/12 (f   + 4f -f   )
                          n+1    n-1         n+1    n  n-1

  From [BIBLI 04].
------------------------------------------------
End of file adambash.txt



4 comments: