Prediction-correction method

'***********************************************************
'*        SOLVING DIFFERENTIAL EQUATIONS OF ORDER 1        *
'*             (Prediction-correction method)              *
'* ------------------------------------------------------- *
'* Reference:                                              *
'*   "Analyse en Turbo Pascal versions 5.5 et 6.0 By Marc  *
'*    DUCAMP et Alain REVERCHON - Eyrolles, Paris 1991"    *
'*    [BIBLI 03].                                          *
'*                                                         *
'*                    Basic version by J-P Moreau, Paris   *
'*                            (www.jpmoreau.fr)            *
'* ------------------------------------------------------- *
'* SAMPLE RUN:                                             *
'* (Integrate y'=(4xy+4x*sqrt(y))/(1+xý) from x=0 to x=10  *
'*  with initial condition y(0)=1)                         *
'*                                                         *
'*  Input x begin: 0                                       *
'*  Input x end  : 10                                      *
'*  Input y begin: 1                                       *
'*  Input number of points: 10                             *
'*  Input finesse: 10                                      *
'*                                                         *
'*         X               Y                               *
'*  -----------------------------                          *
'*     1.000000         8.999984                           *
'*     2.000000        80.999881                           *
'*     3.000000       360.999496                           *
'*     4.000000      1088.998512                           *
'*     5.000000      2600.996484                           *
'*     6.000000      5328.992838                           *
'*     7.000000      9800.986875                           *
'*     8.000000     16640.977767                           *
'*     9.000000     26568.964559                           *
'*    10.000000     40400.946171                           *
'*  -----------------------------                          *
'*                                                         *                                                    *
'***********************************************************
DefInt I-N
DefDbl A-H, O-Z

  MAXDATA = 100

  Dim tx(MAXDATA), ty(MAXDATA)
  'integer ifi, n

  Cls
  Print
  input " Input x begin: ", xi
  input " Input x end  : ", xf
  input " Input y begin: ", yi
  input " Input number of points: ", n
  input " Input finesse: ", ifi

  GoSub 3000  'call equadiff_pc(tx,ty,xi,xf,yi,n,ifi)

' print results
  f$ = "#####.######     #####.######"
  Print
  Print "       X                Y     "
  Print " -----------------------------"
  For i = 1 To n
    Print using; f$; tx(i); ty(i)
  Next i
  Print " -----------------------------"
  Print

End 'of main program

1000  'Function y'=(4xy+4x*sqrt(y))/(1+xý)
  If yy > 0 Then
    f = (4 * xx * yy + 4 * xx * Sqr(yy)) / (1 + xx * xx)
  Else
    f = 0
  End If
Return

2000 'classical Runge-Kutta method of order 4
  xx = X: yy = y: GoSub 1000: a = h * f
  xx = X + h / 2: yy = y + a / 2: GoSub 1000: B = h * f
  yy = y + B / 2: GoSub 1000: c = h * f
  X = X + h: xx = X: yy = y + c: GoSub 1000: d = h * f
  y = y + (a + B + B + c + c + d) / 6
Return

  '**********************************************************
  '*             Prediction-correction method               *
  '* ------------------------------------------------------ *
  '* INPUTS:                                                *
  '*           xi      begin x value                        *
  '*           xf      end x value                          *
  '*           yi      begin y value (at x=xi)              *
  '*           n       number of points to calculate        *
  '*           ifi     finesse (number of intermediate      *
  '*                   points (for example 10)              *
  '* OUTPUTS:                                               *
  '*           tx      table of x values (m values)         *
  '*           ty      table of y values (m values)         *
  '*                                                        *
  '* DESCRIPTION:                                           *
  '* The algorithm has the following steps:                 *
  '*     1. calculate y2, y3, y4 using a classical Runge-   *
  '*        Kutta method of order 4.                        *
  '*     2. from point (x4, y4), first estimate y(n+1) by   *
  '*        formula:                                        *
  '*        y(n+1)=y(n) + h/24(55y'(n)-59y'(n-1)+37y'(n-2)  *
  '*               -9y'(n-3)                                *
  '*        then continue with formula:                     *
  '*        y(n+1)=y(n) + h/24(9y'(n+1)+19y'(n)-5y'(n-1)    *
  '*               +y'(n-2),                                *
  '*        noting that y'(n+1)=f(x(n+1),y(n+1)) with the   *
  '*        estimated value of y(n+1), until convergence is *
  '*        obtained.                                       *
  '**********************************************************
3000 'Subroutine equadiff_pc(tx, ty, xi, xf, yi, n, ifi)
  Dim p(3)
  z = yi
  m = n
  If (m > MAXDATA) Or (ifi < 1) Then Return
  h = (xf - xi) / ifi / m
  xx = xi: yy = yi: GoSub 1000: p(3) = f
  tx(0) = xi: ty(0) = yi
  k = 0
  For i = 1 To m
    ni = (i - 1) * ifi - 1
    For j = 1 To ifi
      X = xi + h * (ni + j)
      k = k + 1
      If (k < 4) Then
        'call runge_kutta(h,x,z)
        y = z: GoSub 2000: z = y
    xx = X: yy = z: GoSub 1000: p(3 - k) = f
      Else
    X = X + h
        w = z + (h / 24) * (55 * p(0) - 59 * p(1) + 37 * p(2) - 9 * p(3))
        'continue
10  y = w
          xx = X: yy = y: GoSub 1000
          w = z + (h / 24) * (9 * f + 19 * p(0) - 5 * p(1) + p(2))
        If Abs(y - w) > 0.0000000001 Then GoTo 10
    z = w: p(3) = p(2): p(2) = p(1)
    p(1) = p(0): xx = X: yy = z: GoSub 1000: p(0) = f
      End If
    Next j 'j loop
    tx(i) = X: ty(i) = z
  Next i 'i loop
Return

' end of file teqdifpc.bas