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