Attribute VB_Name = "Module2"
'****************************************************************
'* View the solutions of a differential equations system having *
'* the form:   y'=f(x,y,z),  z'=g(x,y,z)                        *
'*                                                              *
'*                      Visual Basic 4.0 release by J-P Moreau  *
'*                               (www.jpmoreau.fr)              *
'* ------------------------------------------------------------ *
'* REFERENCE:                                                   *
'* "Graphisme dans le plan et dans l'espace avec Turbo Pascal   *
'*  4.0 By R. Dony - MASSON, Paris 1990, pages 254-260 "        *
'*  [BIBLI 12].                                                 *  
'* ------------------------------------------------------------ *
'* TUTORIAL:                                                    *
'* It can be proved that any differential equation of order > 1 *
'* may be replaced by a 1st order differential equations system.*
'* For example, the 2nd order equation y" = f(x,y,y') is equi-  *
'* valent to the system:                                        *
'*                         y' = z                               *
'*                         z' = f(x,y,z)                        *
'*                                                              *
'* This program draws the solutions (integral curves) of a sys- *
'* tem with the form:                                           *
'*                        y' = f(x,y,z)                         *
'*                        z' = g(x,y,z)                         *
'*                                                              *
'* The proposed example is:    y' = x*x*y*y - 1                 *
'*                             z' = x*x+y*y - 4                 *
'*                                                              *
'* The steps to follow are:                                     *
'*                                                              *
'*     1. choose a physical coordinates window x1,x2,y1,y2      *
'*     2. choose an increment step, dx                          *
'*     3. for each point of plane (x,y) inside previous window, *
'*        calculate by Runge-Kutta method and draw the solution *
'*        curve y=f(x) passing through this point.              *
'*                                                              *
'* Obviously, these curves will eventually go out of screen. It *
'* can be managed by the clipping capability of module Gr2D.    *
'* However to limit computation time, a curve is stopped when:  *
'*                                                              *
'*     - the following point will exit graph window,            *
'*     - the number of segments already drawn > 300,            *
'*     - the curve is approaching a critical point where        *
'*     - f(x,y) and g(x,y) ==> 0.                               *
'*                                                              *
'* Other possible examples:                                     *
'*                                                              *
'*     1) y' = x + y            z' = x*y                        *
'*     2) y' = y(y*y-1)         z' = x(x*x-1)                   *
'*     3) y' = -x               z' = x*x + y*y - 1              *
'*     4) y' = cos(x)-x*sin(y)  z' = x*sin(y) + y*sin(x)        *
'*                                                              *
'* ------------------------------------------------------------ *
'* NOTE: the variable switch (in subroutine Runge-Kutta) is used*
'*       to "switch" direction. For each start point in plane   *
'*       (x,y), the integral curve can be drawn in two opposite *
'*       directions. When switch=1, the curve is drawn in one   *
'*       direction until a stop condition is met; at this moment*
'*       switch is forced to -1 and the other part of the curve *
'*       is drawn until a new stop condition is met. Then switch*
'*       is again put to 1 and the process will continue with a *
'*       new starting point.                                    *
'****************************************************************
Public nbreseg As Integer
Public f1 As Double, f2 As Double
Public f3 As Double, f4 As Double
Public h As Double


' Example:  y' = x*x*y*y-1,  z' = x*x+y*y-4
' Other functions can be defined by user
Function F(X As Double, y As Double)
  F = X * X * y * y - 1#
End Function

' Example:  y' = x*x*y*y-1,  z' = x*x+y*y-4
' Other functions can be defined by user
Function G(X As Double, y As Double)
  G = X * X + y * y - 4#
End Function

'These values can be changed by user
Sub Data()
  Dim s As String * 12
  f1 = -2.5: f2 = 2.5: f3 = -2.5: f4 = 2.5
  s = InputBox("Input increment value:", "DATA", "0.05")
  h = Val(s)
  nbreseg = 300
  'client drawing zone in pixels
  MaxX = 9100: MaxY = 6700
End Sub

' Integration by Runge-Kutta method
Sub RungeKutta(ByVal y0#, ByVal z0#)
   
  Dim y As Double, z As Double
  Dim y1 As Double, z1 As Double
  Dim seg As Integer, switch As Integer
  Dim pointcrit As Double
  Dim k1 As Double, k2 As Double, k3 As Double, k4 As Double
  Dim l1 As Double, l2 As Double, l3 As Double, l4 As Double
  
  switch = 1
  Do
    y = y0: z = z0
    MoveXY y, z
    y1 = y: z1 = z: seg = 0
    Do
       LineXY y1, z1
       seg = seg + 1
       y = y1: z = z1
       k1 = switch * h * F(y, z)
       l1 = switch * h * G(y, z)
       y = y1 + k1 / 2#: z = z1 + l1 / 2#
       k2 = switch * h * F(y, z)
       l2 = switch * h * G(y, z)
       y = y1 + k2 / 2#: z = z1 + l2 / 2#
       k3 = switch * h * F(y, z)
       l3 = switch * h * G(y, z)
       y = y1 + k3 / 2#: z = z1 + l3 / 2#
       k4 = switch * h * F(y, z)
       l4 = switch * h * G(y, z)
       y1 = y + switch * h / 6# * (k1 + 2 * k2 + 2 * k3 + k4)
       z1 = z + switch * h / 6# * (l1 + 2 * l2 + 2 * l3 + l4)
       pointcrit = Sqr(F(y, z) * F(y, z) + G(y, z) * G(y, z))
       y1 = y + switch * h * F(y, z)
       z1 = z + switch * h * G(y, z)
    Loop Until y1 < f1 Or y1 > f2 Or z1 < f3 Or z1 > f4 Or seg >= nbreseg Or pointcrit <= 0.0001
    switch = -switch
  Loop Until switch = 1
End Sub 'RungeKutta()


Sub Initgraph()
  'extern MaxX,MaxY (see gr2d.bas)
  Form1.AutoRedraw = False
  Form1.Cls
  Fenetre f1, f2, f3, f4
  Cloture 750, MaxX - 100, 100, MaxY - 100
  Graduate (f2 - f1) / 10, (f4 - f3) / 10
  Grid (f2 - f1) / 10, (f4 - f3) / 10
  Axes
  Bordure
End Sub

'visualize Runge-Kutta solutions
'see example in functions F and G
Sub Draw_RungeKutta()
  'extern MaxX, MaxY (see gr2d.bas)
  Dim xi As Double, xj As Double
  Dim xpos As Double, ypos As Double

  Data
  Initgraph
  
  xi = 0#
  Do
    xi = xi + h
    xpos = f1 + xi * (f2 - f1)
    xj = 0#
    Do
      xj = xj + h
      ypos = f3 + xj * (f4 - f3)
      RungeKutta xpos, ypos
    Loop Until xj >= 1# - h
  Loop Until xi >= 1# - h
  
  s$ = "Sweep Increment : " & Str$(h)
  Display 250, MaxY + 300, s$
  
  Form1.Command2.Caption = "Continue"

End Sub 'Draw_RungeKutta()