Attribute VB_Name = "Module5"
'************************************************************************
'*                                                                      *
'* Test examples for methods to solve first order ordinary systems of   *
'* differential equations.                                              *
'*                                                                      *
'* We supply the right hand sides, their alpha-numeric description and  *
'* the exact analytic solution, if known.                               *
'*                                                                      *
'* When running the main test program, the user can select among the    *
'* examples below.                                                      *
'*                                                                      *
'*                          Visual Basic Release By J-P Moreau, Paris.  *
'*                                     (www.jpmoreau.fr)                *
'* -------------------------------------------------------------------- *
'* REF.: "Numerical Algorithms with C, By Gisela Engeln-Muellges        *
'*        and Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11].           *
'************************************************************************
' Note: here, only examples #0 and #3 are implemented in Basic
DefDbl A-H, O-Z
DefInt I-N

Option Base 0
  
  
  Sub dgl(bspn As Integer, n, x, y(), f())
    If bspn = 0 Then
        f(0) = y(0) * y(1) + Cos(x) - 0.5 * Sin(2# * x)
        f(1) = y(0) ^ 2 + y(1) ^ 2 - (1# + Sin(x))
    ElseIf bspn = 3 Then
        f(0) = y(1)
        f(1) = y(2)
        f(2) = y(3)
        f(3) = y(4)
        f(4) = (45# * y(2) * y(3) * y(4) - 40# * y(3) ^ 3) / (9# * y(2) ^ 2)
    End If
  End Sub

  Sub settxt(bspn As Integer, n, dgltxt() As String)
    If bspn = 0 Then
      n = 2
      dgltxt(0) = " y1' = y1 * y2 + cos(x) - 0.5 * sin(2.0*x)"
      dgltxt(1) = " y2' = y1 * y1 + y2 * y2 - (1 + sin(x))"
    ElseIf bspn = 3 Then
      n = 5
      dgltxt(0) = " y1' = y2"
      dgltxt(1) = " y2' = y3"
      dgltxt(2) = " y3' = y4"
      dgltxt(3) = " y4' = y5"
      dgltxt(4) = " y5' = (45 * y3 * y4 * y5 - 40 * y4 * y4 * y4) / (9 * y3 * y3)"
    End If
  End Sub

'Other examples given in C

' ----------------------- DE system number 0 -----------------------------

'*************************************************************************
'* Compute the value f of the right hand side of a DE system at (x,y)    *
'*************************************************************************
'static void dgl0(REAL x, REAL *y, REAL *f)
'{
'  f(0) = y(0) * y(1) + COS(x) - HALf * SIN(TWO * x);
'  f(1) = y(0) * y(0) + y(1) * y(1) - (ONE + SIN(x));
'}

'*************************************************************************
'* alpha-numeric description of above system                             *
'*************************************************************************
'static char *dgltxt0(void)
'{
'  return
'    "y1' = y1 * y2 + cos(x) - 0.5 * sin(2.0*x)\n"
'    "y2' = y1 * y1 + y2 * y2 - (1 + sin(x))\n";
'}


' ----------------------- DE system number 1 -----------------------------

'*************************************************************************
'* Compute the value f of the right hand side of a DE system at (x,y)    *
'*************************************************************************
'static void dgl1(REAL x, REAL *y, REAL *f)
'{
'  f(0) = y(1);
'  f(1) = (ONE + y(1) * y(1)) / y(0);
'}

'*************************************************************************
'* alpha-numeric description of above system                             *
'*************************************************************************
'static char *dgltxt1(void)
'{
'  return
'    "y1' = y2\n"
'    "y2' = (1 + y2 * y2) / y1\n";
'}


' ----------------------- DE system number 2 -----------------------------

'*************************************************************************
'* Compute the value f of the right hand side of a DE system at (x,y)    *
'*************************************************************************
'static void dgl2(REAL x, REAL *y, REAL *f)
'{
'  f(0) = y(1);
'  f(1) = -fOUR * y(0) + EXP(x);
'}

'*************************************************************************
'* alpha-numeric description of above system                             *
'*************************************************************************
'static char *dgltxt2(void)
'{
'  return
'    "y1' = y2\n"
'    "y2' = -4 * y1 + exp(x)\n";
'}


'----------------------- DE system number 3 ------------------------------

'*************************************************************************
'* Compute the value f of the right hand side of a DE system at (x,y)    *
'*************************************************************************
'static void dgl3(REAL x, REAL *y, REAL *f)
'{
'  f(0) = y(1);
'  f(1) = y(2);
'  f(2) = y(3);
'  f(3) = y(4);
'  f(4) = ((REAL)45.0 * y(2) * y(3) * y(4) -
'          (REAL)40.0 * y(3) * y(3) * y(3)) / (NINE * y(2) * y(2));
'}

'*************************************************************************
'* alpha-numeric description of above system                             *
'*************************************************************************
'static char *dgltxt3(void)
'{
'  return
'    "y1' = y2\n"
'    "y2' = y3\n"
'    "y3' = y4\n"
'    "y4' = y5\n"
'    "y5' = (45 * y3 * y4 * y5 - 40 * y4 * y4 * y4) / (9 * y3 * y3)\n";
'}



'----------------------- DE system number 4 ------------------------------

'*************************************************************************
'* Compute the value f of the right hand side of a DE system at (x,y)    *
'*************************************************************************
'static void dgl4(REAL x, REAL *y, REAL *f)
'{
'  f(0) = -y(1);
'  f(1) = y(0);
'  f(2) = y(4) - y(1) * y(2);
'  f(3) = y(0) * y(3) - y(5);
'  f(4) = -y(2) - y(1) * y(4);
'  f(5) = y(3) + y(0) * y(5);
'}

'*************************************************************************
'* alpha-numeric description of above system                             *
'*************************************************************************
'static char *dgltxt4(void)
'{
'  return
'    "y1' = -y2\n"
'    "y2' = y1\n"
'    "y3' = y5 - y2 * y3\n"
'    "y4' = y1 * y4 - y6\n"
'    "y5' = -y3 - y2 * y5\n"
'    "y6' = y4 + y1 * y6\n";
'}

'*************************************************************************
'* Compute the value of the analytic solution y(x) for the above DE      *
'* for the initial value problem  y(0) = (1,0,0,1,e,0)                   *
'*************************************************************************
'static void loesung4(REAL x, REAL *y)
'{
'  y(0) = COS(x);
'  y(1) = SIN(x);
'  y(2) = SIN(x) * EXP(COS(x));
'  y(3) = COS(x) * EXP(SIN(x));
'  y(4) = COS(x) * EXP(COS(x));
'  y(5) = SIN(x) * EXP(SIN(x));
'}


'--------------------------- DE system number 5 --------------------------

'*************************************************************************
'* Compute the value f of the right hand side of a DE system at (x,y)    *
'*************************************************************************
'static void dgl5(REAL x, REAL *y, REAL *f)
'{
'  f(0) = (REAL)-500.5 * y(0) + (REAL)499.5 * y(1);
'  f(1) =  (REAL)499.5 * y(0) - (REAL)500.5 * y(1);
'}

'*************************************************************************
'* alpha-numeric description of above system                             *
'*************************************************************************
'static char *dgltxt5(void)
'{
'  return
'    "y1' = -500.5 * y1 + 499.5 * y2\n"
'    "y2' =  499.5 * y1 - 500.5 * y2\n";
'}

'*************************************************************************
'* Compute the value of the analytic solution y(x) for the above DE      *
'* for the initial value problem  y(0) = (4,2)                           *
'*************************************************************************
'static void loesung5(REAL x, REAL *y)
'{
'  y(0) = THREE * EXP(-x) + EXP((REAL)-1000.0 * x);
'  y(1) = THREE * EXP(-x) - EXP((REAL)-1000.0 * x);
'}


'------------------------- DE system number 6 ----------------------------

'*************************************************************************
'* Compute the value f of the right hand side of a DE system at (x,y)    *
'*************************************************************************
'static void dgl6(REAL x, REAL *y, REAL *f)
'{
'  f(0) = (REAL)15.0 * y(1);
'  f(1) = (REAL)-0.6 * y(0);
'}

'*************************************************************************
'* alpha-numeric description of above system                             *
'*************************************************************************
'static char *dgltxt6(void)
'{
'  return
'    "y1' = 15.0 * y2\n"
'    "y2' = -0.6 * y1\n";
'}

'*************************************************************************
'* Compute the value of the analytic solution y(x) for the above DE      *
'* for the initial value problem  y(0) = (0,1)                           *
'*************************************************************************
'static void loesung6(REAL x, REAL *y)
'{
'  y(0) = fIVE * SIN(THREE * x);
'  y(1) = COS(THREE * x);
'}



'------------------------- DE system number 7 ----------------------------

'*************************************************************************
'* Compute the value f of the right hand side of a DE system at (x,y)    *
'*************************************************************************
'static void dgl7(REAL x, REAL *y, REAL *f)
'{
'  f(0) = y(1);
'  f(1) = y(2);
'  f(2) = y(3);
'  f(3) = fOUR * y(3) - SIX * y(2) + fOUR * y(1) - y(0);
'}

'*************************************************************************
'* alpha-numeric description of above system                             *
'*************************************************************************
'static char *dgltxt7(void)
'{
'  return
'    "y1' = y2\n"
'    "y2' = y3\n"
'    "y3' = y4\n"
'    "y4' = 4.0 * y4 - 6.0 * y3 + 4.0 * y2 - y1\n";
'}

'*************************************************************************
'* Compute the value of the analytic solution y(x) for the above DE      *
'* for the initial value problem  y(0) = (0,0,0,6)                       *
'*************************************************************************
'static void loesung7(REAL x, REAL *y)
'{
'  y(0) = x * x * x * EXP(x);
'  y(1) = x * x * (x + THREE) * EXP(x);
'  y(2) = x * (x * x + SIX * x + SIX) * EXP(x);
'  y(3) = (x * x * x + NINE * x * x + (REAL)18.0 * x + SIX) * EXP(x);
'}


'------------------------- DE system number 8 ----------------------------

'*************************************************************************
'* Compute the value f of the right hand side of a DE system at (x,y)    *
'*************************************************************************
'static void dgl8(REAL x, REAL *y, REAL *f)
'{
'  f(0) = (REAL)-52.0 * y(0) + (REAL)50.0 * y(1);
'  f(1) =  (REAL)50.0 * y(0) - (REAL)52.0 * y(1);
'}

'*************************************************************************
'* alpha-numeric description of above system                             *
'*************************************************************************
'static char *dgltxt8(void)
'{
'  return
'    "y1' = -52 * y1 + 50 * y2\n"
'    "y2' =  50 * y1 - 52 * y2\n";
'}

'*************************************************************************
'* Compute the value of the analytic solution y(x) for the above DE      *
'* for the initial value problem                                         *
'*      y(0) = (1265588.55375228,-1265586.55375225)                      *
'*************************************************************************
'static void loesung8(REAL x, REAL *y)
'{
'  y(0) = (REAL)1.000000015       * EXP(        -TWO * x) +
'         (REAL)1265587.553752265 * EXP((REAL)-102.0 * x);
'  y(1) = (REAL)1.000000015 *       EXP(        -TWO * x) -
'         (REAL)1265587.553752265 * EXP((REAL)-102.0 * x);
'}


'------------------------- DE system number 9 ----------------------------

'*************************************************************************
'* Compute the value f of the right hand side of a DE system at (x,y)    *
'*************************************************************************
'static void dgl9(REAL x, REAL *y, REAL *f)
'{
'  f(0) = y(1);
'#if defined(BC3)
'  _fpreset();
'#endif
'  f(1) = -y(0) * COSH(x);
'}


'*************************************************************************
'* alpha-numeric description of above system                             *
'*************************************************************************
'static char *dgltxt9(void)
'{
'  return
'    "y1' = y2\n"
'    "y2' = -y1 * cosh(x)\n";
'}


'------------------------- DE system number 10  --------------------------

'*************************************************************************
'* Compute the value f of the right hand side of a DE system at (x,y)    *
'*************************************************************************
'static void dgl10(REAL x, REAL *y, REAL *f)
'{
'  f(0) = y(1);
'  f(1) = -y(0) * y(0) * y(0);
'}

'*************************************************************************
'* alpha-numeric description of above system                             *
'*************************************************************************
'static char *dgltxt10(void)
'{
'  return
'    "y1' = y2\n"
'    "y2' = -y1^3\n";
'}


'------------------------- DE system number 11 ---------------------------

'*************************************************************************
'* Compute the value f of the right hand side of a DE system at (x,y)    *
'*************************************************************************
'static void dgl11(REAL x, REAL *y, REAL *f)
'{
'#define dk  ((REAL)-10000.0)
'  f(0) = y(1);
'  f(1) = dk * y(0) + (dk - ONE) * y(1);
'#undef  dk
'}

'*************************************************************************
'* alpha-numeric description of above system                             *
'*************************************************************************
'static char *dgltxt11(void)
'{
'  return
'    "y1' =  y1                                    (dk = -10000)\n"
'    "y2' =  dk * y1 + (dk - 1) * y2\n";
'}

'*************************************************************************
'* Compute the value of the analytic solution y(x) for the above DE      *
'* for the initial value problem  (y0) = (1,1)                           *
'*************************************************************************
'static void loesung11(REAL x, REAL *y)
'{
'  y(0) = EXP(-x);
'  y(1) = -EXP(-x);
'}

'end of file t_dgls.bas