1  建立数学模型：微分方程、状态方程、传递函数

2  建立仿真模型：选择仿真算法，将数学模型转化为仿真模型

3  编写运行程序：选择某种语言，写出计算机仿真程序，进行运行，获得仿真结果

一、数学模型

根据前一篇文章的数学模型

不懂 级并串 的参考上一篇文章： 级联、串联、并联求传递函数的方框图和状态方程

二、算法介绍-四阶龙格库塔法（RUNGE-KUTTA）

省略推导公式，直接给出龙格库塔公式：

其中，h为求解步长。

三、VB编程

1.编译环境：VB6.0

2.新建工程

新建CommandButton和Textbox

双击“Command1”，写程序

3.主程序

Private Sub Command1_Click()
Dim h As Single
Dim t As Single
Dim I As Integer
h = 0.1
'h是步长0.1
t = 0
For I = 1 To 200
'I 仿真200步
t = t + h
k11 = h * f1(x1, x2, x3, x4, t)
k21 = h * f2(x1, x2, x3, x4, t)
k31 = h * f3(x1, x2, x3, x4, t)
k41 = h * f4(x1, x2, x3, x4, t)
k12 = h * f1(x1 + k11 / 2, x2 + k21 / 2, x3 + k31 / 2, x4 + k41 / 2, t + h / 2)
k22 = h * f2(x1 + k11 / 2, x2 + k21 / 2, x3 + k31 / 2, x4 + k41 / 2, t + h / 2)
k32 = h * f3(x1 + k11 / 2, x2 + k21 / 2, x3 + k31 / 2, x4 + k41 / 2, t + h / 2)
k42 = h * f4(x1 + k11 / 2, x2 + k21 / 2, x3 + k31 / 2, x4 + k41 / 2, t + h / 2)
k13 = h * f1(x1 + k12 / 2, x2 + k22 / 2, x3 + k32 / 2, x4 + k42 / 2, t + h / 2)
k23 = h * f2(x1 + k12 / 2, x2 + k22 / 2, x3 + k32 / 2, x4 + k42 / 2, t + h / 2)
k33 = h * f3(x1 + k12 / 2, x2 + k22 / 2, x3 + k32 / 2, x4 + k42 / 2, t + h / 2)
k43 = h * f4(x1 + k12 / 2, x2 + k22 / 2, x3 + k32 / 2, x4 + k42 / 2, t + h / 2)
k14 = h * f1(x1 + k13, x2 + k23, x3 + k33, x4 + k43, t + h)
k24 = h * f2(x1 + k13, x2 + k23, x3 + k33, x4 + k43, t + h)
k34 = h * f3(x1 + k13, x2 + k23, x3 + k33, x4 + k43, t + h)
k44 = h * f4(x1 + k13, x2 + k23, x3 + k33, x4 + k43, t + h)
x1 = x1 + (k11 + 2 * k12 + 2 * k13 + k14) / 6
x2 = x2 + (k21 + 2 * k22 + 2 * k23 + k24) / 6
x3 = x3 + (k31 + 2 * k32 + 2 * k33 + k34) / 6
x4 = x4 + (k41 + 2 * k42 + 2 * k43 + k44) / 6
Next I
Text1.Text = x4
End Sub

4.子函数

要计算传递函数的结果，需要将其转化为状态方程，子函数是根据状态方程编写的。

级联法函数

Function f1(x1, x2, x3, x4, t As Single)
f1 = 1 - 7 * x1 - 14 * x2 - 8 * x3
End Function
Function f2(x1, x2, x3, x4, t As Single)
f2 = x1
End Function
Function f3(x1, x2, x3, x4, t As Single)
f3 = x2
End Function
Function f4(x1, x2, x3, x4, t As Single)
f4 = 2 * x2 + 6 * x3
End Function

串联法函数

Function f1(x1, x2, x3, x4, t As Single)
f1 = u - x1
End Function
Function f2(x1, x2, x3, x4, t As Single)
f2 = 1 + 2 * x1 - 2 * x2
End Function
Function f3(x1, x2, x3, x4, t As Single)
f3 = 2 * x2 - 4 * x3
End Function
Function f4(x1, x2, x3, x4, t As Single)
f4 = x3
End Function

并联法函数

Function f1(x1, x2, x3, x4, t As Single)
f1 = 1 - x1
End Function
Function f2(x1, x2, x3, x4, t As Single)
f2 = 1 - 2 * x2
End Function
Function f3(x1, x2, x3, x4, t As Single)
f3 = 1 - 4 * x3
End Function
Function f4(x1, x2, x3, x4, t As Single)
f4 = 4 / 3 * x1 - x2 - 1 / 3 * x3
End Function

5.运行验证

令u=1，即上述子函数中的“1”。

三种方法所得结果相似，侧面证明程序编写正确。

读者可以验证u=t时的不同结果

当遇到写代码无从下手时，可以考虑查询相关算法的例程，相互对照进行借鉴。例如我在其他帖子上看到了用四阶龙格库塔法求一阶微分方程的例子，只有一个未知数，所以用一次5个公式就好了。而本例有4个未知数，所以用5*4=20个公式，以此类推。