龙格库塔法.frm

完整的数值模拟软件

VERSION 5.00
Begin VB.Form Form17 
   BackColor       =   &H00FFFFFF&
   Caption         =   "龙格库塔法"
   ClientHeight    =   8295
   ClientLeft      =   60
   ClientTop       =   345
   ClientWidth     =   11235
   LinkTopic       =   "Form21"
   ScaleHeight     =   8295
   ScaleWidth      =   11235
   StartUpPosition =   3  '窗口缺省
   Begin VB.CommandButton Command1 
      Caption         =   "龙格-库塔算法"
      BeginProperty Font 
         Name            =   "宋体"
         Size            =   10.5
         Charset         =   134
         Weight          =   700
         Underline       =   0   'False
         Italic          =   0   'False
         Strikethrough   =   0   'False
      EndProperty
      Height          =   615
      Left            =   360
      TabIndex        =   0
      Top             =   240
      Width           =   1815
   End
End
Attribute VB_Name = "Form17"
Attribute VB_GlobalNameSpace = False
Attribute VB_Creatable = False
Attribute VB_PredeclaredId = True
Attribute VB_Exposed = False
Private Sub Command1_Click()
Cls
Dim y1(100000), y2(100000), y3(100000), y4(100000)
X = 0
y1(1) = 1.2
y2(1) = 0
y3(1) = 0
y4(1) = -1.04936
u = 0.012129
v = 1 - u
h = InputBox("请输入h:", "")
h = Val(h)
n = 6.2 / h
i = 2
Do While i <= n
   
   X = X + h
   r1 = Sqr((y1(i - 1) + u) * (y1(i - 1) + u) + y2(i - 1) * y2(i - 1))
   r1 = r1 ^ 3
   r2 = Sqr((y1(i - 1) - v) * (y1(i - 1) - v) + y2(i - 1) * y2(i - 1))
   r2 = r2 ^ 3
   k1 = y3(i - 1)
   l1 = y4(i - 1)
   m1 = y1(i - 1) + 2 * y4(i - 1) - v * (y1(i - 1) + u) / r1 - u * (y1(i - 1) - v) / r2
   n1 = y2(i - 1) - 2 * y3(i - 1) - v * y2(i - 1) / r1 - u * y2(i - 1) / r2
   
   r1 = Sqr((y1(i - 1) + h * k1 / 2 + u) * (y1(i - 1) + h * k1 / 2 + u) + (y2(i - 1) + h * l1 / 2) * (y2(i - 1) + h * l1 / 2))
   r1 = r1 ^ 3
   r2 = Sqr((y1(i - 1) + h * k1 / 2 - v) * (y1(i - 1) + h * k1 / 2 - v) + (y2(i - 1) + h * l1 / 2) * (y2(i - 1) + h * l1 / 2))
   r2 = r2 ^ 3
   k2 = y3(i - 1) + h * m1 / 2
   l2 = y4(i - 1) + h * n1 / 2
   m2 = y1(i - 1) + h * k1 / 2 + 2 * (y4(i - 1) + h * n1 / 2) - v * (y1(i - 1) + h * k1 / 2 + u) / r1 - u * (y1(i - 1) + h * k1 / 2 - v) / r2
   n2 = y2(i - 1) + h * l1 / 2 - 2 * (y3(i - 1) + h * m1 / 2) - v * (y2(i - 1) + h * l1 / 2) / r1 - u * (y2(i - 1) + h * l1 / 2) / r2
      
   r1 = Sqr((y1(i - 1) + h * k2 / 2 + u) * (y1(i - 1) + h * k2 / 2 + u) + (y2(i - 1) + h * l2 / 2) * (y2(i - 1) + h * l2 / 2))
   r1 = r1 ^ 3
   r2 = Sqr((y1(i - 1) + h * k2 / 2 - v) * (y1(i - 1) + h * k2 / 2 - v) + (y2(i - 1) + h * l2 / 2) * (y2(i - 1) + h * l2 / 2))
   r2 = r2 ^ 3
   k3 = y3(i - 1) + h * m2 / 2
   l3 = y4(i - 1) + h * n2 / 2
   m3 = y1(i - 1) + h * k2 / 2 + 2 * (y4(i - 1) + h * n2 / 2) - v * (y1(i - 1) + h * k2 / 2 + u) / r1 - u * (y1(i - 1) + h * k2 / 2 - v) / r2
   n3 = y2(i - 1) + h * l2 / 2 - 2 * (y3(i - 1) + h * m2 / 2) - v * (y2(i - 1) + h * l2 / 2) / r1 - u * (y2(i - 1) + h * l2 / 2) / r2
   
   r1 = Sqr((y1(i - 1) + h * k3 + u) * (y1(i - 1) + h * k3 + u) + (y2(i - 1) + h * l3) * (y2(i - 1) + h * l3))
   r1 = r1 ^ 3
   r2 = Sqr((y1(i - 1) + h * k3 - v) * (y1(i - 1) + h * k3 - v) + (y2(i - 1) + h * l3) * (y2(i - 1) + h * l3))
   r2 = r2 ^ 3
   k4 = y3(i - 1) + h * m3
   l4 = y4(i - 1) + h * n3
   m4 = y1(i - 1) + h * k3 + 2 * (y4(i - 1) + h * n3) - v * (y1(i - 1) + h * k3 + u) / r1 - u * (y1(i - 1) + h * 3 - v) / r2
   n4 = y2(i - 1) + h * l3 - 2 * (y3(i - 1) + h * m3) - v * (y2(i - 1) + h * l3) / r1 - u * (y2(i - 1) + h * l3) / r2
   
   y1(i) = y1(i - 1) + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6
   y2(i) = y2(i - 1) + h * (l1 + 2 * l2 + 2 * l3 + l4) / 6
   y3(i) = y3(i - 1) + h * (m1 + 2 * m2 + 2 * m3 + m4) / 6
   y4(i) = y4(i - 1) + h * (n1 + 2 * n2 + 2 * n3 + n4) / 6
   i = i + 1
Loop

Rem 输出坐标轴
Cls
x0 = 7000
y0 = 5000
Line (200, y0)-(15000, y0)
Line (x0, 500)-(x0, 9000)

t0 = 2000

i = 1
Do While i <= n - 1
  Line (x0 + y1(i) * t0, y0 - y2(i) * t0)-(x0 + y1(i + 1) * t0, y0 - y2(i + 1) * t0)
  i = i + 1
Loop
End Sub