Attribute VB_Name = "Module1"
'************************************************************************
'*                                                                      *
'* Solve an ordinary system of first order differential equations using *
'* -------------------------------------------------------------------- *
'* automatic step size control                                          *
'* ----------------------------                                         *
'*                                                                      *
'* Programming language: ANSI C                                         *
'* Author:               Klaus Niederdrenk (FORTRAN)                    *
'* Adaptation:           Juergen Dietel, Computer Center, RWTH Aachen   *
'* Source:               existing C, Pascal, QuickBASIC and FORTRAN     *
'*                       codes                                          *
'* Date:                 6.2.1992, 10.2.1995                            *
'*                                                                      *
'*                       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].           *
'************************************************************************

DefDbl A-H, O-Z
DefInt I-N

Option Base 0

Dim yhilf()
Dim k1() As Double, k2() As Double
Dim k3() As Double, k4() As Double
Dim k5() As Double, k6() As Double
Dim k7() As Double
Dim g6() As Double, g7() As Double

Dim y_bad(), y_good()


'debug only
Sub PrintVec(Name As String, n, V())
  Form1.Print Name
  For i = 0 To n - 1
    Form1.Print " "; V(i);
  Next i
  Form1.Print
End Sub


'Maximum norm of a difference vector ........
Function dist_max(vector1(), vector2(), n)
'************************************************************************
'* Compute the maximum norm of the difference of two [0..n-1] vectors   *
'*                                                                      *
'* global name used:                                                    *
'* ================                                                     *
'*   None                                                               *
'************************************************************************
' double  abstand     reference value for computation of distance
' double  hilf        distance of two vector elements
  abstand = 0#
  For i = n - 1 To 0 Step -1
    hilf = Abs(vector1(i) - vector2(i))
    If hilf > abstand Then abstand = hilf
  Next i
  dist_max = abstand
End Function

' Runge-Kutta embedding formulas of 2nd, 3rd degree ....................
Sub ruku23(x, y(), bspn As Integer, n, h, y2(), y3())
'************************************************************************
'* Compute 2nd and 3rd order approximates y2, y3 at x + h starting with *
'* a solution y at x by using Runge-Kutta embedding formulas on the     *
'* first order system of n differential equations   y' = f(x,y) , as    *
'* supplied by  dgl().                                                  *
'*                                                                      *
'* Input parameters:                                                    *
'* =================                                                    *
'* x    x-value of left end point                                       *
'* y    y-values at x                                                   *
'* bspn # example                                                       *
'* n    number of differential equations                                *
'* dgl  function that evaluates the right hand side of the system       *
'*      y' = f(x,y)                                                     *
'* h    step size                                                       *
'*                                                                      *
'* yhilf,k1,k2,k3: auxiliary vectors defined in module awp.             *
'*                                                                      *
'* Output parameters:                                                   *
'* ==================                                                   *
'* y2   2nd order approximation for y at x + h                          *
'* y3   3rd order approximation for y at x + h                          *
'*                                                                      *
'************************************************************************
  ReDim yhilf(n)
  ReDim k1(n) As Double
  ReDim k2(n) As Double
  ReDim k3(n) As Double
  Call dgl(bspn, n, x, y, k1)
  For i = 0 To n - 1
    yhilf(i) = y(i) + h * k1(i)
  Next i
  Call dgl(bspn, n, x + h, yhilf, k2)
  For i = 0 To n - 1
    yhilf(i) = y(i) + 0.25 * h * (k1(i) + k2(i))
  Next i
  Call dgl(bspn, n, x + 0.5 * h, yhilf, k3)
  For i = 0 To n - 1
    y2(i) = y(i) + 0.5 * h * (k1(i) + k2(i))
    y3(i) = y(i) + h / 6# * (k1(i) + k2(i) + 4# * k3(i))
  Next i
End Sub

' England's Einbettungs formulas of 4th and 5th degree ................
Sub engl45(x, y(), bspn As Integer, n, h, y4(), y5())
' double    x                   starting point of integration ........
' double    y(0:n-1)            initial value at x ...................
' integer   bspn,               # example
'           n                   number of differential equations .....
' double    h                   step size ............................
' double   y4(0:n-1),           4th order approximation for y at x + h
'          y5(0:n-1)            5th order approximation for y at x + h
' auxiliary vectors
ReDim yhilf(n)
ReDim k1(n) As Double
ReDim k2(n) As Double
ReDim k3(n) As Double
ReDim k4(n) As Double
ReDim k5(n) As Double
ReDim k6(n) As Double
'************************************************************************
'* Compute 4th and 5th order approximates y4, y5 at x + h starting with *
'* a solution y at x by using the England embedding formulas on the     *
'* first order system of n differential equations   y' = f(x,y) , as    *
'* supplied by  dgl().                                                  *
'*                                                                      *
'* Input parameters:                                                    *
'* =================                                                    *
'* x    initial x-value                                                 *
'* y    y-values at x, type pVEC                                        *
'* n    number of differential equations                                *
'* dgl  function that evaluates the right hand side of the system       *
'*      y' = f(x,y)                                                     *
'* h    step size                                                       *
'*                                                                      *
'* yhilf, k1..K6: auxiliary vectors                                     *
'*                                                                      *
'* Output parameters:                                                   *
'* ==================                                                   *
'* y4   4th order approximation for y at x + h (pVEC)                   *
'* y5   5th order approximation for y at x + h (pVEC)                   *
'*                                                                      *
'************************************************************************
  Call dgl(bspn, n, x, y, k1)
  For i = 0 To n - 1
    yhilf(i) = y(i) + 0.5 * h * k1(i)
  Next i
  Call dgl(bspn, n, x + 0.5 * h, yhilf, k2)
  For i = 0 To n - 1
    yhilf(i) = y(i) + (0.25 * h * (k1(i) + k2(i)))
  Next i
  Call dgl(bspn, n, x + 0.5 * h, yhilf, k3)
  For i = 0 To n - 1
    yhilf(i) = y(i) + h * (-k2(i) + 2# * k3(i))
  Next i
  Call dgl(bspn, n, x + h, yhilf, k4)
  For i = 0 To n - 1
    yhilf(i) = y(i) + h / 27# * (7# * k1(i) + 10# * k2(i) + k4(i))
  Next i
  Call dgl(bspn, n, x + (2# / 3#) * h, yhilf, k5)
  For i = 0 To n - 1
    yhilf(i) = y(i) + h / 625# * (28# * k1(i) - 125# * k2(i) + 546# * k3(i) + 54# * k4(i) - 378# * k5(i))
  Next i
  Call dgl(bspn, n, x + h / 5#, yhilf, k6)
  For i = 0 To n - 1
    y4(i) = y(i) + h / 6# * (k1(i) + 4# * k3(i) + k4(i))
    y5(i) = y(i) + h / 336# * (14# * k1(i) + 35# * k4(i) + 162# * k5(i) + 125# * k6(i))
  Next i
End Sub

' embedding formulas of Prince-Dormand of 4./5. order .........................
Sub prdo45(x, y(), bspn As Integer, n, h, y4(), y5(), steif1 As Integer, steifanz As Integer, steif2 As Integer)
' double    x,                      starting point of integration .....
'           y(0:n-1)                initial value at x ................
' integer   bspn,                   # example .........................
'           n                       number of DEs .....................
' double    h,                      step size .........................
'           y4(0:n-1),              solution of 4th order at x+h ......
'           y5(0:n-1)               solution of 5th order at x+h ......
' auxiliary flags
' integer   steif1, steifanz,steif2
' auxiliary vectors
ReDim yhilf(n)
ReDim k1(n) As Double, k2(n) As Double
ReDim k3(n) As Double, k4(n) As Double
ReDim k5(n) As Double, k6(n) As Double
ReDim k7(n) As Double
ReDim g6(n) As Double, g7(n) As Double
'************************************************************************
'* Compute 4th and 5th order approximates y4, y5 at x + h starting with *
'* a solution y at x by using the Prince-Dormand embedding formulas on  *
'* the first order system of n differential equations y' = f(x,y) , as  *
'* supplied by  dgl().                                                  *
'* Simultaneously we perform two tests for stiffness whose results are  *
'* stored in steif1 and steif2.                                         *
'*                                                                      *
'* Input parameters:                                                    *
'* =================                                                    *
'* x    initial x-value                                                 *
'* y    y-values at x (pVEC)                                            *
'* n    number of differential equations                                *
'* dgl  function that evaluates the right hand side of the system       *
'*      y' = f(x,y)                                                     *
'* h    step size                                                       *
'*                                                                      *
'* yhilf, k1..k7,g6,g7: auxiliary vectors.                              *
'*                                                                      *
'* Output parameters:                                                   *
'* ==================                                                   *
'* y4   4th order approximation for y at x + h                          *
'* y5   5th order approximation for y at x + h                          *
'*                                                                      *
'***********************************************************************}
  Dim steifa As Integer  'Flag which is set if the second test for stiffness
                         'Shampine und Hiebert) is positive; otherwise the
                         'flag is erased.

  Call dgl(bspn, n, x, y, k1)
  For i = 0 To n - 1
    yhilf(i) = y(i) + 0.2 * h * k1(i)
  Next i
  Call dgl(bspn, n, x + 0.2 * h, yhilf, k2)
  For i = 0 To n - 1
    yhilf(i) = y(i) + 0.075 * h * (k1(i) + 3# * k2(i))
  Next i
  Call dgl(bspn, n, x + 0.3 * h, yhilf, k3)
  For i = 0 To n - 1
    yhilf(i) = y(i) + h / 45# * (44# * k1(i) - 168# * k2(i) + 160# * k3(i))
  Next i
  Call dgl(bspn, n, x + 0.8 * h, yhilf, k4)
  For i = 0 To n - 1
    yhilf(i) = y(i) + h / 6561# * (19372# * k1(i) - 76080# * k2(i) + 64448# * k3(i) - 1908# * k4(i))
  Next i
  Call dgl(bspn, n, x + (8# / 9#) * h, yhilf, k5)
  For i = 0 To n - 1
    g6(i) = y(i) + h / 167904# * (477901# * k1(i) - 1806240# * k2(i) + 1495424# * k3(i) + 46746# * k4(i) - 45927# * k5(i))
  Next i
  Call dgl(bspn, n, x + h, g6, k6)
  For i = 0 To n - 1
    g7(i) = y(i) + h / 142464# * (12985# * k1(i) + 64000# * k3(i) + 92750# * k4(i) - 45927# * k5(i) + 18656# * k6(i))
  Next i
  Call dgl(bspn, n, x + h, g7, k7)
  For i = 0 To n - 1
    y5(i) = g7(i)
    y4(i) = y(i) + h / 21369600# * (1921409# * k1(i) + 969088# * k3(i) + 13122270# * k4(i) - 5802111# * k5(i) + 1902912# * k6(i) + 534240# * k7(i))
  Next i

' Test for stiffness via dominant eigenvalue

  If dist_max(k7, k6, n) > 3.3 * dist_max(g7, g6, n) Then steif1 = 1

' one step in steffness test of Shampine & Hiebert

  For i = 0 To n - 1
    g6(i) = h * (2.2 * k2(i) + 0.13 * k4(i) + 0.144 * k5(i))
    g7(i) = h * (2.134 * k1(i) + 0.24 * k3(i) + 0.1 * k6(i))
  Next i

  If dist_max(g6, g7, n) < dist_max(y4, y5, n) Then
    steifa = 1
  Else
    steifa = 0
  End If

  If (steifa > 0) Then
    steifanz = steifanz + 1
    If steifanz >= 3 Then steif2 = 1
  Else
    steifanz = 0
  End If

End Sub

' Find the maximum norm of a REAL vector ............................
Function norm_max(vektor(), n) As Double
'      double  vektor(0:n-1)                 vector .................
'      integer n                             length of vector .......
' ************************************************************************
' *          Return the maximum norm of a [0..n-1] vector v              *
' *                                                                      *
' ************************************************************************
Dim norm As Double                            ' local max.
'   double betrag                             ' magnitude of a component
    norm = 0#
    For i = 0 To n - 1
      betrag = Abs(vektor(i))
      If betrag > norm Then norm = betrag
    Next i
    norm_max = norm
End Function

' 1st order DESs with automatic step size control ...........................
Sub awp(x, xend, bspn As Integer, n, y(), epsabs, epsrel, h, methode, fmax As Integer, aufrufe As Integer, fehler As Integer)
'  double x,                      initial/final x value ..............
'         xend                    desired end point ..................
' integer bspn,                   # example
'         n                       number of DEs ......................
'  double y(0:n-1),               initial/final y value ..............
'         epsabs,                 absolute error bound ...............
'         epsrel,                 relative error bound ...............
'         h                       initial/final step size ............
' integer methode,                desired method (3, 6, 7) ...........
'         fmax,                   maximal # of calls of  dgl() .......
'         aufrufe,                actual # of calls of  dgl() ........
'         fehler                  error code .........................
'************************************************************************
'* Compute the solution y of a system of first order ordinary           *
'* differential equations       y' = f(x,y)   at xend from the given    *
'* initial data (x0, y0).                                               *
'* We use automatic step size control internally so that the error of   *
'* y (absolutely or relatively) lies within the given error bounds      *
'* epsabs and epsrel.                                                   *
'*                                                                      *
'* Input parameters:                                                    *
'* =================                                                    *
'* x        initial value for x                                         *
'* y        initial values for y(0:n-1)                                 *
'* bspn     # example                                                   *
'* n        number of differential equations                            *
'* dgl      function that evaluates the right hand side of the system   *
'*          y' = f(x,y)  (see t_dgls) (Removed from list of parameters) *
'* xend     end of integration; xend > x0                               *
'* h        initial step size                                           *
'* epsabs   absolute error bound; >= 0; if = 0 we only check the        *
'*          relative error.                                             *
'* epsrel   relative error bound; >= 0; if = 0 we check only the        *
'*          absolute eror.                                              *
'* fmax     max number of evaluations of right hand side in dgl()       *
'* methode  chooses the method                                          *
'*          = 3: Runge-Kutta method of 2nd/3rd order                    *
'*          = 6: England formula of 4th/5th order                       *
'*          = 7: Formula of Prince-Dormand of 4th/5th order             *
'*                                                                      *
'* Output parameters:                                                   *
'* ==================                                                   *
'* x        final x-value for iteration. If fehler = 0  we usually have *
'*            x = xend.                                                 *
'* y        final y-values for the solution at x                        *
'* h        final step size used; leave for subsequent calls            *
'* aufrufe  actual number of calls of dgl()                             *
'*                                                                      *
'* Return value (fehler):                                               *
'* =====================                                                *
'* = 0: all ok                                                          *
'* = 1: both error bounds chosen too small for the given mach. constant *
'* = 2: xend <= x0                                                      *
'* = 3: h <= 0                                                          *
'* = 4: n <= 0                                                          *
'* = 5: more right hand side calls than allowed: aufrufe > fmax,        *
'*      x and h contain the current values when stop occured.           *
'* = 6: improper input for embedding formula                            *
'* = 7: lack of available memory (not used here)                        *
'* = 8: Computations completed, but the Prince Dormand formula stiff-   *
'*      ness test indicates possible stiffness.                         *
'* = 9: Computations completed, but both Prince Dormand formula stiff-  *
'*      ness tests indicate possible stiffness. Use method for stiff    *
'*      systems instead '                                               *
'* =10: aufrufe > fmax, see error code 5; AND the Prince Dormand formula*
'*      indicates stiffness; retry using a stiff DE solver '            *
'*                                                                      *
'************************************************************************
Dim MACH_EPS As Double
Dim MACH_2 As Double
Dim mach_1 As Double  'machine constant dependent variable which
                      'avoids using too little steps near xend.
Dim xend_h As Double  '|xend| - MACH_2, carrying same sign as xend
Dim ymax As Double    'Maximum norm of newest approximation of max
                      'order.
Dim hhilf As Double   'aux storage for the latest value of h
                      'produced by step size control. It is saved
                      'here in order to avoid to return a `h' that
                      'resulted from an arbitrary reduction at the
                      'end of the interval.
Dim diff As Double    'distance of the two approximations from the
                      'embedding formula.
Dim s As Double       'indicates acceptance level for results from
                      'embeding formula.

'approximate solution of low order
ReDim y_bad(n)
'ditto of high order
ReDim y_good(n)
   

Dim amEnde As Integer 'flag that shows if the end of the interval
                      'can be reached with the actual step size.
Dim fertig As Integer 'flag indicating end of iterations.

Dim steif1 As Integer 'Flag, that is set in prdo45() if its
                      'stiffness test (dominant eigenvalue)
                      'indicates so. Otherwise no changes.
Dim steifanz As Integer 'counter for number of successive successes
                        'of stiffness test of Shampine and Hiebert in
                        'prdo45().
Dim steif2 As Integer 'Flag, set in prdo45(), when the stiffness
                      'test of  Shampine and Hiebert wa successful
                      'three times in a row; otherwise no changes.
   
'initialize some variables
  fehler = 0
  MACH_EPS = 1.2E-16
  MACH_2 = 100 * MACH_EPS
  mach_1 = MACH_EPS ^ 0.75
  amEnde = 0
  fertig = 0
  steif1 = 0
  steif2 = 0
  steifanz = 0
  aufrufe = 1
  ymax = norm_max(y, n)

  If xend >= 0# Then
    xend_h = xend * (1# - MACH_2)
  Else
    xend_h = xend * (1# + MACH_2)
  End If

' ----------------------- check inputs ----------------------
  If epsabs <= MACH_2 * ymax And epsrel <= MACH_2 Then
    fehler = 1
    Exit Sub
  End If
  If xend_h < x Then
    fehler = 2
    Exit Sub
  End If
  If h < MACH_2 * Abs(x) Then
    fehler = 3
    Exit Sub
  End If
  If n <= 0 Then
    fehler = 4
    Exit Sub
  End If
  If methode <> 3 And methode <> 6 And methode <> 7 Then
    fehler = 6
    Exit Sub
  End If
 
' **********************************************************************
' *                                                                    *
' *                       I t e r a t i o n s                          *
' *                                                                    *
' **********************************************************************
  If x + h > xend_h Then
                                    'almost at end point ?
    hhilf = h                       'A shortened step might be
    h = xend - x                    'enough.
    amEnde = 1
  End If

  While fertig = 0                  'solve DE system by integrating from
                                    'x0 to xend by suitable steps.
    'choose method
    If methode = 3 Then
      Call ruku23(x, y, bspn, n, h, y_bad, y_good)
    ElseIf methode = 6 Then
      Call engl45(x, y, bspn, n, h, y_bad, y_good)
    ElseIf methode = 7 Then
      Call prdo45(x, y, bspn, n, h, y_bad, y_good, steif1, steifanz, steif2)
    End If

    aufrufe = aufrufe + methode

    diff = dist_max(y_bad, y_good, n)

    If (diff < MACH_2) Then         'compute s
      s = 2#
    Else
      ymax = norm_max(y_good, n)
      s = Sqr(h * (epsabs + epsrel * ymax) / diff)
      If methode <> 3 Then s = Sqr(s)
    End If

    If s > 1# Then                  'integration acceptable ?
      For i = 0 To n - 1            'accept highest order solution
        y(i) = y_good(i)            'move x
      Next i
      
      x = x + h

      If amEnde <> 0 Then           'at end of interval ?
        fertig = 1                  'stop iteration
        If methode = 7 Then
          If steif1 > 0 Or steif2 > 0 Then fehler = 8
          If steif1 > 0 And steif2 > 0 Then fehler = 9
        End If
      ElseIf aufrufe > fmax Then    'too many calls of   dgl() ?
        hhilf = h                   'save actual step size
        fehler = 5                  'report error and stop
        fertig = 1
        If methode = 7 And (steif1 > 0 Or steif2 > 0) Then fehler = 10
      Else                          'Integration was successful
                                    'not at the interval end ?
        h = h * XMin(2#, 0.98 * s)  'increase step size for next
                                    'step properly, at most by
                                    'factor two. Value `0.98*s' is
                                    'used in order to avoid that
                                    'the theoretical value s is
                                    'exceeded by accidental
                                    'rounding errors.
        If x + h > xend_h Then      'nearly reached xend ?
          hhilf = h                 '=> One further step with
          h = xend - x              'reduced step size might be
          amEnde = 1                'enough.
          If h < mach_1 * Abs(xend) Then  'very close to xend ?
            fertig = 1                    'finish iteration.
          End If
        End If
      End If
    Else                            'step unsuccessful ?
                                    'before repeating this step
      h = h * XMax(0.5, 0.98 * s)   'reduce step size properly, at
                                    'most by factor 1/2 (for factor
      amEnde = 0                    '0.98: see above).
    End If

  Wend
 
  h = hhilf           'return the latest step size computed by step
                      'size control and  error code to the caller.

End Sub


Function XMax(a, b)
  If a >= b Then
    XMax = a
  Else
    XMax = b
  End If
End Function

Function XMin(a, b)
  If a <= b Then
    XMin = a
  Else
    XMin = b
  End If
End Function

' -------------------------- END  uawp.bas ----------------------------