Friday, September 22, 2017

The Runge-Kutta Method

'********************************************************************
'*             Differential equations of order 1                    *
'*             by Runge-Kutta method of order 4                     *
'*                                                                  *
'*                             Basic version by J-P Moreau, Paris   *
'*                                     (www.jpmoreau.fr)            *
'* ---------------------------------------------------------------- *
'* Reference: "Analyse en Turbo Pascal versions 5.5 et 6.0 By Marc  *
'*             DUCAMP et Alain REVERCHON - Eyrolles, Paris 1991"    *
'*             [BIBLI 03].                                          *
'* ---------------------------------------------------------------- *
'* SAMPLE RUN:                                                      *
'*                                                                  *
'* Example: integrate y'=4x*(y+sqrt(y))/1+x^2 from x=0 to x=1       *
'*                                                                  *
'*        DIFFERENTIAL EQUATION WITH 1 VARIABLE OF ORDER 1          *
'*                     of type y' = f(x,y)                          *
'*                                                                  *
'*   begin value x   : ? 0                                          *
'*   end value x     : ? 1                                          *
'*   y value at x0   : ? 1                                          *
'*   number of points: ? 11                                         *
'*   finesse         : ? 10                                         *
'*                                                                  *
'*       X            Y                                             *
'* -------------------------                                        *
'*   0.000000     1.000000                                          *
'*   0.100000     1.040400                                          *
'*   0.200000     1.166400                                          *
'*   0.300000     1.392400                                          *
'*   0.400000     1.742400                                          *
'*   0.500000     2.250000                                          *
'*   0.600000     2.958400                                          *
'*   0.700000     3.920400                                          *
'*   0.800000     5.198400                                          *
'*   0.900000     6.864400                                          *
'*   1.000000     9.000000                                          *
'*                                                                  *
'********************************************************************
DefInt I-N
DefDbl A-H, O-Z

'ifi:INTEGER

  Cls
  Print
  Print "    DIFFERENTIAL EQUATION WITH 1 VARIABLE OF ORDER 1"
  Print "             of type y' = f(x,y)"
  Print
  print "    begin value x   : "; : input xi
  print "    end value x     : "; : input xf
  print "    y value at x0   : "; : input yi
  print "    number of points: "; : input m
  print "    finesse         : "; : input ifi

  Dim t(m + 1)

  'equadiff1(fp,t,xi,xf,yi,m,ifi)
  GoSub 2000

End


1000 'Example: y'=4x(y+rac(y))/(1+x²)
    fp = 4# * xx * (yy + Sqr(yy)) / (1# + xx * xx)
Return

'***************************************************************************
'*        SOLVING DIFFERENTIAL EQUATIONS WITH 1 VARIABLE OF ORDER 1        *
'*                       of type y' = f(x,y)                               *
'* ----------------------------------------------------------------------- *
'*  INPUTS:                                                                *
'*    fp        Given equation to integrate (see test program)             *
'*    xi, xf    Begin, end values of variable x                            *
'*    yi        Begin Value of y at x=xi                                   *
'*    m         Number of points to calculate                              *
'*    fi        finesse (number of intermediate points)                    *
'*                                                                         *
'*  OUTPUTS:                                                               *
'*    t       : real vector storing m results for function y               *
'***************************************************************************
2000 'Subroutine  Equadiff1
  If ifi < 1 Then Return
  h = (xf - xi) / ifi / (m - 1)
  Y = yi
  t(1) = yi

  For i = 1 To m
    ni = (i - 1) * ifi - 1
    For j = 1 To ifi
      X = xi + h * (ni + j)
      xx = X: yy = Y
      GoSub 1000 'calculate current fp(x,y)
      a = h * fp
      xx = X + h / 2#: yy = Y + a / 2#
      GoSub 1000 'calculate fp(x+h/2,y+a/2)
      B = h * fp
      yy = Y + B / 2#
      GoSub 1000 'calculate fp(x+h/2,y+b/2)
      c = h * fp
      X = X + h
      xx = X: yy = Y + c
      GoSub 1000 'calculate fp(x,y+c)
      d = h * fp
      Y = Y + (a + B + B + c + c + d) / 6
    Next j
    t(i + 1) = Y
  Next i
  'Affiche(t,m,xi,xf)
  GoSub 3000
Return

3000 'Subroutine Affiche
  h = (xf - xi) / (m - 1)
  X = xi - h
  Cls
  Print
  Print "      X         Y     "
  Print "----------------------"
  For i = 1 To m
    X = X + h
    Print using; " ##.######  ##.######"; X; t(i)
  Next i
Return

'End of file teqdif1.bas



EXPLANATION FILE OF PROGRAM TEQUDIF1
====================================

  
  The Runge-Kutta Method
  ----------------------


  We present here the Runge-Kutta method of order 4 to integrate an ODE

  of order 1:       Y' = F(X, Y)

  The development of Y around x  coincidates with its Taylor development
                               n
  of order 4:

       y    = y  + h y' + (h^2/2) y" + (h^3/6) y"' + (h^4/24) y""
        n+1    n      n            n            n              n

  Other orders may also be used for the Runge-Kutta method, for instance

  the Euler's method is a Runge-Kutta of order one.

  Order 4 is often used because it is a good compromise between speed and
  
  accuracy.

  From a current point (x , y ), the next point is determined by:
                         n   n

   K  = h f(x , y )
         0        n   n

 k  = h f(x + h/2, y + K /2)
         1        n        n   0

 k  = h f(x + h/2, y + K /2)
         2        n        n   1

 k  = h f(x + h, y + K )
         3        n      n   2

  and   y   = y  + (1/6)(K + 2 K + K )
         n+1   n          0     1   3

  It can be proved that the error at each step is about h^5.

  The only drawback is the time of calculation: at each step, four function

  evaluations are necessary.

  From [BIBLI 03].

-------------------------------------------------------------
End of file rkutta.txt.

Wednesday, September 20, 2017

Runge-Kutta method of order 4

'********************************************************************
'*        Differential equations with p variables of order 1        *
'*             by Runge-Kutta method of order 4                     *
'* ---------------------------------------------------------------- *
'* Reference: "Analyse en Turbo Pascal versions 5.5 et 6.0 by Marc  *
'*             DUCAMP et Alain REVERCHON - Eyrolles, Paris 1991"    *
'*                                                                  *
'*                        Basic version 1.1 By J-P Moreau, Paris    *
'*                                  (www.jpmoreau.fr)               *
'* ---------------------------------------------------------------- *
'* SAMPLE RUN:                                                      *
'*                                                                  *
'* Example #1: integrate system of equations from x=0 to x=3:       *
'*              y1' = y2 + y3 - 3*y1                                *
'*              y2' = y1 + y3 - 3*y2                                *
'*              y3' = y1 + y2 - 3*y3                                *
'* with initial conditions: y1(0)=1, y2(0)=2 and y3(0)=-1           *
'*                                                                  *
'*                                                                  *
'*        DIFFERENTIAL EQUATION WITH P VARIABLE OF ORDER 1          *
'*              of type yi' = f(y1,y2,...,yn), i=1..n               *
'*                                                                  *
'*   number of variables: 3                                         *
'*   begin value x      : 0                                         *
'*   end value x        : 3                                         *
'*   y1 value at x0     : 1                                         *
'*   y2 value at x0     : 2                                         *
'*   y3 value at x0     : -1                                        *
'*   number of points   : 7                                         *
'*   finesse            : 30                                        *
'*                                                                  *
'*       X        Y1         Y2         Y3                          *
'* --------------------------------------------------------         *
'*   0.000000   1.000000   2.000000  -1.000000                      *
'*   0.500000   0.449466   0.584801   0.178795                      *
'*   1.000000   0.251358   0.269674   0.214727                      *
'*   1.500000   0.149580   0.152058   0.144622                      *
'*   2.000000   0.090335   0.090671   0.089664                      *
'*   2.500000   0.054738   0.054784   0.054648                      *
'*   3.000000   0.033193   0.033200   0.033181                      *
'* --------------------------------------------------------         *
'*                                                                  *
'* Example #2: integrate system of equations from x=0 to PI*Sqrt(2) *
'*              y1' = y2                                            *
'*              y2' = -4*y1 - 3*y3                                  *
'*              y3' = y4                                            *
'*              y4' = -8*y1 - 2*y3                                  *
'* with initial conditions: y1(0)=3, y2(0)=0, y3(0)=4, y4(0)=0      *
'*                                                                  *
'*                                                                  *
'*        DIFFERENTIAL EQUATION WITH P VARIABLE OF ORDER 1          *
'*              of type yi' = f(y1,y2,...,yn), i=1..n               *
'*                                                                  *
'*   number of variables: 4                                         *
'*   begin value x      : 0                                         *
'*   end value x        : 4.442883                                  *
'*   y1 value at x0     : 3                                         *
'*   y2 value at x0     : 0                                         *
'*   y3 value at x0     : 4                                         *
'*   y4 value at x0     : 0                                         *
'*   number of points   : 9                                         *
'*   finesse            : 30                                        *
'*                                                                  *
'*       X        Y1         Y2         Y3         Y4               *
'* --------------------------------------------------------         *
'*   0.000000   3.000000   0.000000   4.000000   0.000000           *
'*   0.555360   0.000000  -8.485281   0.000000 -11.313708           *
'*   1.110721  -3.000000  -0.000001  -4.000000  -0.000002           *
'*   1.666081  -0.000001   8.485281  -0.000001  11.313708           *
'*   2.221442   3.000000   0.000003   4.000000   0.000003           *
'*   2.776802   0.000001  -8.485281   0.000002 -11.313708           *
'*   3.332162  -3.000000  -0.000004  -4.000000  -0.000005           *
'*   3.887523  -0.000002   8.485281  -0.000002  11.313708           *
'*   4.442883   3.000000   0.000005   4.000000   0.000007           *
'* --------------------------------------------------------         *
'*                                                                  *
'* Release 1.1: Added example #2 (Sept. 2008).                      *
'********************************************************************
DefInt I-N
DefDbl A-H, O-Z

'ifi,ip : INTEGER
Dim y(10), yi(10), t1(50), t2(50), t3(50), t4(50)

  Cls
  Print
  Print "    DIFFERENTIAL EQUATIONS WITH P VARIABLES OF ORDER 1"
  Print "        of type yi' = f(y1,y2,...,yn), i=1..n"
  Print
  PRINT "    number of variables : "; : INPUT ip
  Print
  PRINT "    begin value x       : "; : INPUT xi
  PRINT "    end value x         : "; : INPUT xf
  For i = 0 To ip - 1
    PRINT "    y"; i + 1; " value at x0    : "; : INPUT yi(i)
  Next i
  PRINT "    number of points    : "; : INPUT m
  PRINT "    finesse             : "; : INPUT ifi

  'call subroutine eqdifp
  GoSub 2000

End


'Example #1: y1'=y2+y3-3y1, y2'=y1+y3-3y2, y3'=y1+y2-3y3
'1000 'FUNCTION fp
'  IF k = 0 THEN
'   fp = y(1) + y(2) - 3# * y(0)
'  ELSEIF k = 1 THEN
'   fp = y(0) + y(2) - 3# * y(1)
'  ELSEIF k = 2 THEN
'   fp = y(0) + y(1) - 3# * y(2)
'  ELSE
'   fp = 0#
'  END IF
'RETURN

'Example #2: y1'=y2, y2'=-4y1-3y3, y'3=y4, y'4=-8y1-2y3
1000 'FUNCTION fp
  If k = 0 Then
    fp = y(1)
  ElseIf k = 1 Then
    fp = -4# * y(0) - 3# * y(2)
  ElseIf k = 2 Then
    fp = y(3)
  ElseIf k = 3 Then
    fp = -8# * y(0) - 2# * y(2)
  Else
    fp = 0#
  End If
Return

'****************************************************************************
'*         SOLVING DIFFERENTIAL SYSTEMS WITH P VARIABLES OF ORDER 1         *
'*                 of type yi' = f(y1,y2,...,yn), i=1..n                    *
'* ------------------------------------------------------------------------ *
'*  INPUTS:                                                                 *
'*    m         number of points to calculate                               *
'*    xi, xf    begin, end values of variable x                             *
'*    yi        table of begin values of functions at xi                    *
'*    ip        number of independant variables                             *
'*    ifi       finesse (number of intermediary points)                     *
'*                                                                          *
'*  OUTPUTS:                                                                *
'*    t1,t2     real vectors storing the results for first two functions,   *
'*              y1 and y2.                                                  *
'* ------------------------------------------------------------------------ *
'*  EXAMPLE:    y1'=y2+y3-3y1, y2'=y1+y3-3y2, y3'=y1+y2-3y3                 *
'*              Exact solution :  y1 = 1/3 (exp(-4x)  + 2 exp(-x))          *
'*                                y2 = 1/3 (4exp(-4x) + 2 exp(-x))          *
'*                                y3 = 1/3 (-5exp(-4x)+ 2 exp(-x))          *
'****************************************************************************
2000 'subroutine Eqdifp
  Dim ta(10), tb(10), tc(10), td(10), z(10)
  h = (xf - xi) / ifi / (m - 1)
   ip = ip - 1
   t1(1) = yi(0)
   t2(1) = yi(1)
   t3(1) = yi(2)
   t4(1) = yi(3)
   For k = 0 To ip
     y(k) = yi(k): z(k) = yi(k)
   Next k
   For i = 1 To m
     ni = (i - 1) * ifi - 1
     For j = 1 To ifi
       X = xi + h * (ni + j)
       For k = 0 To ip
         y(k) = z(k)
       Next k
       For k = 0 To ip
         GoSub 1000
         ta(k) = h * fp
       Next k
       For k = 0 To ip
         y(k) = z(k) + ta(k) / 2#
       Next k
       X = X + h / 2#
       For k = 0 To ip
         GoSub 1000
         tb(k) = h * fp
       Next k
       For k = 0 To ip
         y(k) = z(k) + tb(k) / 2#
       Next k
       For k = 0 To ip
         GoSub 1000
         tc(k) = h * fp
       Next k
       For k = 0 To ip
         y(k) = z(k) + tc(k)
       Next k
       X = X + h / 2#
       For k = 0 To ip
         GoSub 1000
         td(k) = h * fp
       Next k
       For k = 0 To ip
         z(k) = z(k) + (ta(k) + 2# * tb(k) + 2# * tc(k) + td(k)) / 6#
       Next k
     Next j
     t1(i + 1) = z(0)
     t2(i + 1) = z(1)
     t3(i + 1) = z(2)
     t4(i + 1) = z(3)
   Next i
   'Display(t1,t2,t3,t4,m,p,xi,xf)
   GoSub 3000
Return

3000 'subroutine Display
  h = (xf - xi) / (m - 1)
  X = xi - h
  Print
  Print "      X";
  For i = 1 To ip + 1
    Print "       Y"; i;
  Next i
  Print
  Print "--------------------------------------------------------"
  For i = 1 To m
    X = X + h
    If ip = 1 Then
      Print USING; " ##.######  ##.######  ##.######"; X; t1(i); t2(i)
    ElseIf ip = 2 Then
      Print USING; " ##.######  ##.######  ##.######  ##.######"; X; t1(i); t2(i); t3(i)
    ElseIf ip > 2 Then
      Print USING; " ##.######  ##.######  ##.######  ##.######  ###.######"; X; t1(i); t2(i); t3(i); t4(i)
    End If
  Next i
  Print "--------------------------------------------------------"
Return

'End of file teqdifp.bas

Monday, September 18, 2017

Adams-Moulton Prediction-Correction Method

'********************************************************************
'*   Solve Y' = F(X,Y) with initial conditions using the Adams-     *
'*   Moulton Prediction-Correction Method                           *
'* ---------------------------------------------------------------- *
'* SAMPLE RUN                                                       *
'* (Integrate Y' = -Y + X/(1+X)^2 from X=0 to X=1 with initial      *
'*  condition Y(0) = 1 )                                            *
'*                                                                  *
'*     X          Y         Y True   |Y-Y True|                     *
'* ---------------------------------------------                    *
'*  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.869569   0.869565    0.000004   1                  *
'*  0.200000   0.833340   0.833333    0.000006   1                  *
'*  0.250000   0.800008   0.800000    0.000009   1                  *
'*  0.300000   0.769241   0.769231    0.000010   1                  *
'*  0.350000   0.740752   0.740741    0.000011   1                  *
'*  0.400000   0.714298   0.714286    0.000012   1                  *
'*  0.450000   0.689668   0.689655    0.000012   1                  *
'*  0.500000   0.666679   0.666667    0.000013   1                  *
'*  0.550000   0.645174   0.645161    0.000013   1                  *
'*  0.600000   0.625013   0.625000    0.000013   1                  *
'*  0.650000   0.606073   0.606061    0.000013   1                  *
'*  0.700000   0.588248   0.588235    0.000013   1                  *
'*  0.750000   0.571441   0.571429    0.000013   1                  *
'*  0.800000   0.555568   0.555556    0.000012   1                  *
'*  0.850000   0.540553   0.540541    0.000012   1                  *
'*  0.900000   0.526327   0.526316    0.000012   1                  *
'*  0.950000   0.512832   0.512821    0.000011   1                  *
'*  1.000000   0.500011   0.500000    0.000011   1                  *
'* ---------------------------------------------------------------- *
'* REFERENCE: "Méthode de calcul numérique- Tome 2 - Programmes en  *
'*             Basic et en Pascal By Claude Nowakowski, Edition du  *
'*             P.S.I., 1984".                                       *
'*                                                                  *
'*                      Quick Basic Release By J-P Moreau, Paris.   *
'*                                  (www.jpmoreau.fr)               *
'********************************************************************
'See explanation file adambash.txt
'---------------------------------
DefInt I-N
DefDbl A-H, O-Z

Option Base 0        'index begins from zero

Dim X(3), Y(3)

H = 0.05             'integration step
X(0) = 0#: Y(0) = 1# 'Initial conditions
EC = 0.000001        'Precision

F$ = " ##.######  ##.######  ##.######   ##.######  ##"

Cls
Print "     X          Y         Y True   |Y-Y True| "
Print " ---------------------------------------------"
Print USING; F$; X(0); Y(0); Y(0); X(0)

'Start with Runge-Kutta
For K = 0 To 1
  XX = X(K): YY = Y(K): GoSub 1000: C1 = C
  XX = X(K) + H / 2#: YY = Y(K) + H / 2# * C1: GoSub 1000: C2 = C
  YY = Y(K) + H / 2# * C2: GoSub 1000: C3 = C
  X(K + 1) = X(K) + H
  XX = X(K + 1): YY = Y(K) + H * C3: GoSub 1000: C4 = C
  Y(K + 1) = Y(K) + H * (C1 + 2 * C2 + 2 * C3 + C4) / 6#
  XX = X(K + 1): GoSub 2000: ER = Abs(F - Y(K + 1))

  Print USING; F$; X(K + 1); Y(K + 1); F; ER

Next K

100 K = 2
XX = X(K): YY = Y(K): GoSub 1000: C1 = C
XX = X(K - 1): YY = Y(K - 1): GoSub 1000: C2 = C
XX = X(K - 2): YY = Y(K - 2): GoSub 1000: C3 = C
X(K + 1) = X(K) + H
YP = Y(K) + H / 12# * (23 * C1 - 16 * C2 + 5 * C3)

L = 0
200 XX = X(K + 1): YY = YP: GoSub 1000: C1 = C
XX = X(K): YY = Y(K): GoSub 1000: C2 = C
XX = X(K - 1): YY = Y(K - 1): GoSub 1000: C3 = C
YC = Y(K) + H / 12# * (5 * C1 + 8 * C2 - C3)

'PRINT YC

If Abs(YP - YC) > EC Then
  YP = YC: L = L + 1: GoTo 200
End If

Y(K + 1) = YC: XX = X(K + 1): GoSub 2000: ER = Abs(F - Y(K + 1))

Print USING; F$; X(K + 1); Y(K + 1); F; ER; L


For K = 0 To 2
  X(K) = X(K + 1): Y(K) = Y(K + 1)
Next K

If X(2) < 1# Then
  GoTo 100
End If

End

1000 C = -YY + XX / ((1# + XX) ^ 2)
Return

2000 F = 1# / (1# + XX)
Return

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