'***********************************************************************
'* Integrate a System of Ordinary Differential Equations By the        *
'* Runge-Kutta-Fehlberg method (double precision)                      *
'* ------------------------------------------------------------------- *
'* REFERENCE:     H A Watts and L F Shampine,                          *
'*                Sandia Laboratories,                                 *
'*                Albuquerque, New Mexico.                             *
'* ------------------------------------------------------------------- *
'*                                                                     *
'* Basic Release 1.1 By J-P Moreau, Paris.                             *
'* FreeBASIC adaptation by J. Debord                                   *
'*                                                                     *
'* Release 1.1: added test #3.                                         *
'***********************************************************************

#INCLUDE "fbmath.bi"

SUB DiffEq1(T AS DOUBLE, Y() AS DOUBLE, Yp() AS DOUBLE)
' Differential equation for Test 1

  Yp(1) = .25 * Y(1) * (1 - Y(1) / 20)
END SUB

FUNCTION Yexact1(T AS DOUBLE) AS DOUBLE
' Exact solution of the ODE for Test 1

  Yexact1 = 20 / (1 + 19 * EXP(-.25 * T))
END FUNCTION

SUB DiffEq2(T AS DOUBLE, Y() AS DOUBLE, Yp() AS DOUBLE)
' Differential equations for Test 2

  Yp(1) = Y(2)
  Yp(2) = - Y(1)
END SUB

SUB Yexact2(T AS DOUBLE, BYREF Y1 AS DOUBLE, BYREF Y2 AS DOUBLE)
' Exact solution of the ODE's for Test 2

  Y1 = COS(T)
  Y2 = - SIN(T)
END SUB

SUB DiffEq3(T AS DOUBLE, Y() AS DOUBLE, Yp() AS DOUBLE)
' Differential equations for Test 3

  Yp(1) = Y(2)
  Yp(2) = Y(3)
  Yp(3) = Y(4)
  Yp(4) = Y(5)
  Yp(5) = (45 * Y(3) * Y(4) * Y(5) - 40 * Y(4) * Y(4) * Y(4)) / (9 * Y(3) * Y(3))
END SUB

SUB Test1

  CONST Neqn = 1                           ' Number of equations

  DIM AS DOUBLE  Y(1 TO Neqn) = {1}        ' Functions and initial cond.
  DIM AS DOUBLE  Yp(1 TO Neqn)             ' Derivatives

  DIM AS DOUBLE  Tstart = 0, Tstop = 20    ' Integration interval
  DIM AS INTEGER Nstep  = 5                ' Number of steps
  DIM AS DOUBLE  StepSize = _
                 (Tstop - Tstart) / Nstep  ' Step size

  DIM AS DOUBLE  AbsErr = 1D-6, _
                 RelErr = 1D-6             ' Abs. and relative errors

  DIM AS INTEGER Flag = 1                  ' Error flag

  DIM AS DOUBLE  T = Tstart
  DIM AS DOUBLE  Tout
  DIM AS INTEGER I
  DIM AS DOUBLE  Yc

  PRINT
  PRINT " TEST01"
  PRINT " Solve a scalar equation:"
  PRINT
  PRINT "  Y' = 0.25 * Y * ( 1 - Y / 20 )"
  PRINT

  PRINT " T            Y            Y exact            Error"
  PRINT

  Yc = Yexact1(T)

  PRINT USING "##"; T;
  PRINT USING "       #.####^^^^"; Y(1); Yc; Y(1) - Yc

  FOR I = 1 TO Nstep
    Tout = T + StepSize
    RKF45 @DiffEq1, Y(), Yp(), T, Tout, RelErr, AbsErr, Flag
    Yc = Yexact1(Tout)
    PRINT USING "##"; T;
    PRINT USING "       #.####^^^^"; Y(1); Yc; Y(1) - Yc
    T = Tout
  NEXT I
END SUB

SUB Test2

  CONST Neqn = 2                             ' Number of equations

  DIM AS DOUBLE  Y(1 TO Neqn) = {1, 0}       ' Functions and initial cond.
  DIM AS DOUBLE  Yp(1 TO Neqn)               ' Derivatives

  DIM AS DOUBLE  Tstart = 0, Tstop = 2 * Pi  ' Integration interval
  DIM AS INTEGER Nstep  = 12                 ' Number of steps

  DIM AS DOUBLE  AbsErr = 1D-6, _
                 RelErr = 1D-6               ' Abs. and relative errors

  DIM AS INTEGER Flag = 1                    ' Error flag

  DIM AS DOUBLE  StepSize = (Tstop - Tstart) / Nstep
  DIM AS DOUBLE  T = Tstart
  DIM AS DOUBLE  Tout
  DIM AS INTEGER I
  DIM AS DOUBLE  Yc1, Yc2

  PRINT
  PRINT " TEST02"
  PRINT " Solve a vector equation:"
  PRINT
  PRINT "  Y'(1) =   Y(2)"
  PRINT "  Y'(2) = - Y(1)"
  PRINT

  PRINT "  T          Y1         Y1 exact         Y2          Y2 exact"
  PRINT

  Yexact2 T, Yc1, Yc2

PRINT USING "##.##"; T;
  PRINT USING "    #.####^^^^"; Y(1); Yc1; Y(2); Yc2

  FOR I = 1 TO Nstep
    Tout = T + StepSize
    RKF45 @DiffEq2, Y(), Yp(), T, Tout, RelErr, AbsErr, Flag
    Yexact2 Tout, Yc1, Yc2
    PRINT USING "##.##"; T;
    PRINT USING "    #.####^^^^"; Y(1); Yc1; Y(2); Yc2
    T = Tout
  NEXT I
END SUB

SUB Test3

  CONST Neqn = 5                                 ' Number of equations

  DIM AS DOUBLE  Y(1 TO Neqn) = {1, 1, 1, 1, 1}  ' Functions and initial cond.
  DIM AS DOUBLE  Yp(1 TO Neqn)                   ' Derivatives

  DIM AS DOUBLE  Tstart = 0, Tstop = 1.5         ' Integration interval
  DIM AS INTEGER Nstep  = 11                     ' Number of steps

  DIM AS DOUBLE  AbsErr = 1D-6, _
                 RelErr = 1D-6                   ' Abs. and relative errors

  DIM AS INTEGER Flag = 1                        ' Error flag

  DIM AS DOUBLE  StepSize = (Tstop - Tstart) / Nstep
  DIM AS DOUBLE  T = Tstart
  DIM AS DOUBLE  Tout
  DIM AS INTEGER I

  PRINT
  PRINT " TEST03"
  PRINT " Solve a vector equation:"
  PRINT
  PRINT "  Y'(1) = Y(2)"
  PRINT "  Y'(2) = Y(3)"
  PRINT "  Y'(3) = Y(4)"
  PRINT "  Y'(4) = Y(5)"
  PRINT "  Y'(5) = (45 * Y(3) * Y(4) * Y(5) - 40 * Y(4)^3) / (9 * Y(3)^2)"

  PRINT
  PRINT "  T          Y1            Y2            Y3            Y4            Y5"
  PRINT

  PRINT USING "##.##"; T;
  PRINT USING "    #.####^^^^"; Y(1); Y(2); Y(3); Y(4); Y(5)

  FOR I = 1 TO Nstep
    Tout = T + StepSize
    RKF45 @DiffEq3, Y(), Yp(), T, Tout, RelErr, AbsErr, Flag
    PRINT USING "##.##"; T;
    PRINT USING "    #.####^^^^"; Y(1); Y(2); Y(3); Y(4); Y(5)
    T = Tout
  NEXT I
END SUB

CLS
PRINT
PRINT " PROGRAM TEST_RKF"
PRINT " Demonstrate the RKF45 ODE integrator."

Test1
PRINT : PRINT "Any key to continue... "
SLEEP

Test2
PRINT : PRINT "Any key to continue... "
SLEEP

Test3
PRINT : PRINT "Any key to continue... "
SLEEP

END 'of main program
