計算機仿真龍格庫塔方法的軟件VB設計與實現

計算機仿真（龍格庫塔方法）的軟件 VB 設計與實現 Dim a(0 To 10) As Single, b(0 To 10) As Single, c(0 To 10) As Single, d(0 To 10) As Single Dim e(0 To 10) As Single, h(0 To 10) As Single, p(0 To 10) As Single, q(0 To 10) As Single Dim n1 As Byte, n2 As Byte n1 表示的是 y 的階數 ,n2 表示的是輸入函數的階數 Dim i As Integer Dim f(0 To 2000) As Single f(m)=y, 也就是各個時刻 的y值 Dim X0(0 To 12) As Single, X1(0 To 12) As Single Dim dt As Single, u As Single dt 為采樣周期 ,u 為輸入 Dim q1 As Single, p1 As Single, h1 As Single, e1 As Single Dim b0(0 To 10) As Single, bn(-4 To 10) As Single Dim m As Integer Dim max As Single Private Sub Combo1_Click() n1 表示的是 y 的階數 n1 = Combo1.ListIndex + 1 For i = n1 + 1 To 8 Text1(i).Visible = False Label1(i).Visible = False Label9(i).Visible = False Next For i = 0 To n1 Text1(i).Visible = True Label1(i).Visible = True Label9(i).Visible = True Next Combo2.Clear Combo2.Text = " 請選擇輸出 u 的最大階數 " For i = 0 To n1 Combo2.AddItem ((i) & " 階 ") Next End Sub n2 表示的是輸入函數的階數 Private Sub Combo2_Click() n2 = Combo2.ListIndex For i = n2 + 1 To 8 Text2(i).Visible = False Label2(i).Visible = False Label10(i).Visible = False Next For i = 0 To n2 Text2(i).Visible = True Label2(i).Visible = True Label10(i).Visible = True Next End Sub Private Sub Calculate_Click() dt = Val(Text3.Text) dt 是采樣周期 u = Val(Text4.Text) n1 = Combo1.ListIndex + 1 If n1 = 0 Then 對付忘記選擇最高階數時的情況 Call MsgBox(" 請選擇階數！ ", 48, " 未選擇階數 ") Exit Sub End If For i = 0 To 12 設X0的初值都是0,且i大于8，這一點很重要! X0(i) = 0 Next For i = 0 To n1 - 1 a(i)b(i)c(i)d(i) 等等都是計算的中間 a(i) = -Val(Text1(i).Text) / Val(Text1(n1).Text) 變量 Next For i = 0 To n2 b0(i) = Val(Text2(i).Text) / Val(Text1(n1).Text) Next If Combo3.ListIndex = 1 Then n1 = n1 + 1 For i = n1 - 1 To 1 Step -1 a(i) = a(i - 1) Next a(0) = 0 End If For i = -4 To -1 bn(i) = 0 Next b(0) = a(0) * a(n1 - 1) For i = 1 To n1 - 1 b(i) = a(i - 1) + a(i) * a(n1 - 1) Next c(0) = a(0) * b(n1 - 1) For i = 1 To n1 - 1 c(i) = b(i - 1) + a(i) * b(n1 - 1) Next d(0) = a(0) * c(n1 - 1) For i = 1 To n1 - 1 d(i) = c(i - 1) + a(i) * c(n1 - 1) Next For i = 0 To n1 - 1 e(i) = dt * a(i) + dt * dt / 2 * b(i) + dt * dt * dt / 6 * c(i) + dt * dt * dt * dt / 24 * d(i) h(i) = dt * dt / 2 * a(i) + dt * dt * dt / 6 * b(i) + dt * dt * dt * dt / 24 * c(i) p(i) = dt * dt * dt / 6 * a(i) + dt * dt * dt * dt / 24 * b(i) q(i) = dt * dt * dt * dt / 24 * a(i) Next If n2 = 0 Then n2 為 0 和不為 0 的公式是不一樣的，表示輸入函數 u 為 n2 階 b0(0) = Val(Text2(0).Text) / Val(Text1(n1).Text) m = 1 Do q1 = 0: p1 = 0: h1 = 0: e1 = 0 For i = 0 To n1 - 1 這些個也是中間變量 q1 = q1 + q(i) * X0(i + 1) p1 = p1 + p(i) * X0(i + 1) h1 = h1 + h(i) * X0(i + 1) e1 = e1 + e(i) * X0(i + 1) Next For i = 1 To n1 X1(i) = X0(i) + dt * X0(i + 1) + dt * dt / 2 * X0(i + 2) + dt * dt * dt / 6 * X0(i + 3) + dt * dt * dt * dt / 24 * X0(i + 4) Next X1(n1) = X1(n1) + e1 + (dt + dt * dt / 2 * a(n1 - 1) + dt * dt * dt / 6 * b(n1 - 1) + dt * dt * dt * dt / 24 * c(n1 - 1)) * b0(0) * u If n1 >= 4 Then X1(n1 - 3) = X1(n1 - 3) + q1 + dt * dt * dt * dt / 24 * b0(0) * u X1(n1 - 2) = X1(n1 - 2) + p1 + (dt * dt * dt / 6 + dt * dt * dt * dt / 24 * a(n1 - 1)) * b0(0) * u X1(n1 - 1) = X1(n1 - 1) + h1 + (dt * dt / 2 + dt * dt * dt / 6 * a(n1 - 1) + dt * dt * dt * dt / 24 * b(n1 - 1)) * b0(0) * u ElseIf n1 >= 3 Then X1(n1 - 2) = X1(n1 - 2) + p1 + (dt * dt * dt / 6 + dt * dt * dt * dt / 24 * a(n1 - 1)) * b0(0) * u X1(n1 - 1) = X1(n1 - 1) + h1 + (dt * dt / 2 + dt * dt * dt / 6 * a(n1 - 1) + dt * dt * dt * dt / 24 * b(n1 - 1)) * b0(0) * u ElseIf n1 >= 2 Then X1(n1 - 1) = X1(n1 - 1) + h1 + (dt * dt / 2 + dt * dt * dt / 6 * a(n1 - 1) + dt * dt * dt * dt / 24 * b(n1 - 1)) * b0(0) * u End If f(m) = X1(1) For i = 1 To n1 X0(i) = X1(i) Next m = m + 1 Loop Until m = 2000 循環(huán)計算 2000 次 Else For i = (n2 + 1) To 10 b0(i) = 0 Next For i = 0 To n1 - 1 bn(i) = b0(i) + b0(n1) * a(i) Next Do For i = 1 To n1 X1(i) = X0(i) + dt * X0(i + 1) + dt * dt / 2 * X0(i + 2) + dt * dt * dt / 6 * X0(i + 3) + dt * dt * dt * dt / 24 * X0(i + 4) + (h(n1 - i) * bn(n1 - 1) + p(n1 - i) * bn(n1 - 2) + q(n1 - i) * bn(n1 - 3)) * u Next For i = 1 To n1 X1(i) = X1(i) + e(n1 - i) * X0(1) + h(n1 - i) * X0(2) + p(n1 - i) * X0(3) + q(n1 - i) * X0(4) + (dt * bn(n1 - i) + dt * dt / 2 * bn(n1 - 1 - i) + dt * dt * dt / 6 * bn(n1 - 2 - i) + dt * dt * dt * dt / 24 * bn(n1 - 3 - i)) * u Next f(m) = X1(1) For i = 1 To n1 X0(i) = X1(i) Next m = m + 1 Loop Until m = 2000 End If max = f(1) For m = 2 To 1999 開始繪圖 If Abs(f(m)) > Abs(max) Then max = f(m) End If Next If max > 0 Then Picture1.Cls Picture1.Scale (-20 * dt, 1.2 * max)-(2020 * dt, -0.05 * max) Picture1.Line (0, 0)-(2002 * dt, 0) Picture1.Line (0, 0)-(0, 1.18 * max) Picture1.Line (2002 * dt, 0)-(1960 * dt, -0.025 * max) Picture1.Line (2002 * dt, 0)-(1960 * dt, 0.025 * max) Picture1.Line (0, 1.18 * max)-(10 * dt, 1.1 * max) Picture1.Line (0, 1.18 * max)-(-10 * dt, 1.1 * max) Else Picture1.Cls Picture1.Scale (-5 * dt, 0.2 * Abs(max))-(2002 * dt, 1.2 * max) Picture1.Line (0, 0)-(2002 * dt, 0) Picture1.Line (0, 0)-(0, 1.2 * max) Picture1.Line (2002 * dt, 0)-(1900 * dt, 0.1 * max) Picture1.Line (0, 1.2 * max)-(50 * dt, 1.1 * max) End If For m = 2 To 1999 Picture1.Line ((m - 1) * dt, f(m - 1))-(m * dt, f(m)), vbRed Next If (f(1999) - f(1800)) < 0.00001 Then Picture1.Line (0, f(1999))-(2000 * dt, f(1999)) Picture1.CurrentX = 0 Picture1.CurrentY = f(1999) + 0.1 * max Picturel.Print "y g =" & f(1999) End If If Combo3.ListIndex = 0 Then Picture1.Line (0, f(1999))-(2000 * dt, f(1999)) Picture1.CurrentX = 0 Picture1.CurrentY = f(1999) + 0.1 * max Picture1.Print "y g =" & f(1999) Else Picture1.Line (0, 0)-(2000 * dt, u * 2000 * dt) End If End Sub Private Sub Combo3_Click() Text4.Visible = True Label12.Visible = True If Combo3.ListIndex = 0 Then Label5.Caption ="請輸入階躍信號 u的系數：” Else Label5.Caption = " 請輸入斜坡信號 u 的系數： End If End Sub Private Sub dtxnzb_Click() Dim y00 As Single Dim thigema As Single Dim ts As Single Dim ess As Single If Combo3.List In dex = 0 Then yOO = f(1999) thigema = max - y00 m = 2000 Do Loop Un til Abs(f(m) - y00) >= (0.05 * yOO) ts = m * dt ess = u - y00 Framel.Visible = =True Label6.Capti on = ="超調量d :"& thigema Label7.Capti on = ="調節(jié)時間 ts:" & ts Label8.Capti on = ="穩(wěn)態(tài)誤差 ess:" & ess Else Framel.Visible = True Label6.Caption ="穩(wěn)態(tài)誤差 ess:" & (u * 1999 * dt - f(1999)) / (u * 1999 * dt) & "%" End If End Sub