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()