【计算机仿真技术】代码记录
计算机仿真技术代码记录
第一节课课程练习
报童问题：报童每卖一份报可赚a元，每天卖K份报的概率为 $ P_k $ ，若订的报卖不出去，每退一份报赔钱b元,试求报童每天期望收益达到最大的订报量Z。
代码：
Private Sub Command1_Click()
Dim s As Single
Dim z As Integer
s = 40
For z = 1 To 9
s1=0
For t = 1 To 100
m = Rnd()
If m < 0.025 Then K = 1
If m < 0.075 And m>=0.025 Then K =2
If m < 0.175 And m>=0.075 Then K =3
If m < 0.35 And m>=0.175 Then K =4
If m < 0.65 And m>=0.35 Then K =5
If m < 0.825 And m>=0.65 Then K =6
If m < 0.925 And m>=0.825 Then K =7
If m < 0.975And m>=0.925 Then K =8
If m <1 And m>=0.675 Then K =9
If z > K Then
S1 = S1 + (z - K) * 40
Else
S1 = S1 + (K - z) * 20
End If
Next
If s >= S1 / 100 Then s = S1 / 100
z0 = z
End If
Next
Text1.Text = s
Text2.Text = z0
End Sub
运行结果：
s = 34.4
z0 = 4
第二节课程练习
采用4阶龙格库塔法仿真如下微分方程,步长0.1s,仿真100步。
\[0.075 x^{\prime \prime \prime}+0.75 x^{\prime \prime}+x^{\prime}+5 x=5
\]
程序如下：
Function f1(x1, x2, x3, x4, t As Single)
f1 = x2
End Function
Function f2(x1, x2, x3, x4, t As Single)
f2 = x3
End Function
Function f3(x1, x2, x3, x4, t As Single)
f3 = (5 - 5 * x1 - x2 - 0.75 * x3) / 0.075
End Function
Function f4(x1, x2, x3, x4, t As Single)
f4 = x1
End Function
Private Sub Command1_Click()
Dim h As Single
Dim t As Single
Dim I As Integer
'Dim u As Single
h = 0.1
'h是步长0.1
t = 0
'u = 1
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阶龙格库塔法仿真如下微分方程,步长0.1s,仿真100步。
\[\begin{array}{l}
x_1^{\prime} =-x_1+u \\
x_2^{\prime} =-3 x_2 +u \\
x_3^{\prime} =-4 \times 3+u \\
y=x 1+x_2-2 \times 3 \\
u = 2t
\end{array}
\]
代码：
Function f1(x1, x2, x3, x4, t As Single)
f1 = -x1 + 2 * t
End Function
Function f2(x1, x2, x3, x4, t As Single)
f2 = -3 * x2 + 2 * t
End Function
Function f3(x1, x2, x3, x4, t As Single)
f3 = -4 * x3 + 2 * t
End Function
Function f4(x1, x2, x3, x4, t As Single)
f4 = x1 + x2 - 2 * x3
End Function
Private Sub Command1_Click()
Dim h As Single
Dim t As Single
Dim I As Integer
'Dim u As Single
h = 0.01
'h是步长0.1
t = 0
'u = 1
For I = 1 To 100
'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 = x1
End Sub
运行结果：
x1 = 0.7484009
第四节课堂练习
仿真如下微分方程,步长0.1s,仿真100步。
\[\begin{array}{l}
{x_1}^{\prime}=x_2 +1 \\
{x_2}^{\prime} = x_1 -2x_2 -1 \\
x_1(0)=x_2(0)=0
\end{array}
\]
代码如下：
Public Sub Command1_Click()
Dim x1(1000) As Single
Dim x2(1000) As Single
For n = 1 To 100
x1(n + 1) = x1(n) + 0.01 * x2(n) + 0.01 + 0.00005 * (-x1(n) - 2 * x2(n) - 1)
x2(n + 1) = 0.9802 * x2(n) - 0.0099 * (x1(n) + 1) - 0.0000497 * (x2(n) + 1)
Next
Text1.Text = x1(100)
Text2.Text = x2(100)
End Sub
运行结果：
x1(100) = 0.6284227
x2(100) = -0.6284116
第5节课堂练习
分别使用双线性变换法、RK4法、Z域法、时域法对如下控制系统进行仿真，步长0.1s,仿真200步。
简单变换下有：
\[\begin{array}{l}
w^{\prime}=1-y \\
x^{\prime}=-0.1 x+1-y+0.5 \mathrm{w} \\
z^{\prime}=-2 z+2 x \\
y^{\prime}=-10 y+10 z
\end{array}
\]
仿真代码：
Function f1(x1, x2, x3, x4, t As Single)
f1 = 1 - x4
End Function
Function f2(x1, x2, x3, x4, t As Single)
f2 = 0.5 * x1 - 0.1 * x2 + 1 - x4
End Function
Function f3(x1, x2, x3, x4, t As Single)
f3 = 2 * x2 - 2 * x3
End Function
Function f4(x1, x2, x3, x4, t As Single)
f4 = 10 * x3 - 10 * x4
End Function
Private Sub Command1_Click()
'RK4法
Dim h As Single
Dim t As Single
Dim I As Integer
h = 0.01
'h是步长0.1
t = 0
'u = 1
For I = 1 To 200
'I 仿真100步
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
Private Sub Command2_Click()
'双线性方法
Dim u(1000) As Single
Dim w(1000) As Single
Dim x(1000) As Single
Dim z(1000) As Single
Dim y(1000) As Single
t = 0.01
For n = 1 To 200
u(n + 1) = 1 - y(n)
w(n + 1) = w(n) + t * u(n + 1) / 2 + t * u(n) / 2
x(n + 1) = ((0.5 * t + 2) * w(n + 1) + (0.5 * t - 2) * w(n) - (0.1 * t - 2) * x(n)) / (0.1 * t + 2)
z(n + 1) = (2 * t * x(n + 1) + 2 * t * x(n) - (2 * t - 2) * z(n)) / (2 * t + 2)
y(n + 1) = (10 * t * z(n + 1) + 10 * t * z(n) + (2 - 10 * t) * y(n)) / (2 + 10 * t)
Next
Text2.Text = y(201)
End Sub
Private Sub Command3_Click()
'Z域法
Dim u(1000) As Single
Dim w(1000) As Single
Dim x(1000) As Single
Dim z(1000) As Single
Dim y(1000) As Single
t = 0.01
For n = 1 To 200
u(n + 1) = 1 - y(n)
w(n + 1) = w(n) + t * u(n)
x(n + 1) = Exp(-0.1 * t) * x(n) + w(n + 1) + (4 - 5 * Exp(-0.1 * t)) * w(n)
z(n + 1) = Exp(-2 * t) * z(n) + (1 - Exp(-2 * t)) * x(n)
y(n + 1) = Exp(-10 * t) * y(n) + (1 - Exp(-10 * t)) * z(n)
Next
Text3.Text = y(201)
End Sub
Private Sub Command4_Click()
'时域法
Dim u(1000) As Single
Dim w(1000) As Single
Dim x(1000) As Single
Dim z(1000) As Single
Dim y(1000) As Single
t = 0.01
For n = 1 To 200
u(n + 1) = 1 - y(n)
w(n + 1) = w(n) + t * (1 - y(n))
'x(n + 1) = Exp(-0.1 * t) * x(n) - (10 * Exp(-0.1 * t) - 10) * (0.5 * w(n) + 1 - y(n))
x(n + 1) = Exp(-0.1 * t) * x(n) - (10 * Exp(-0.1 * t) - 10) * w(n) + (10 * t - (100 - 100 * Exp(0.1 * t))) * (1 - y(n))
z(n + 1) = Exp(-2 * t) * z(n) + (1 - Exp(-2 * t)) * x(n)
y(n + 1) = Exp(-10 * t) * y(n) + (1 - Exp(-10 * t)) * z(n)
Next
Text4.Text = y(201)
End Sub
第六节课堂练习
对下图所示系统进行仿真。
其中：
数字控制器部分传递函数
\[\begin{array}{l}
D(z) = \frac{E_2(z)}{E_1(z)}=\frac{z+0.21}{z+0.632}
\end{array}
\]
非线性部分：
\[\begin{equation}
f(x)=
\begin{cases}
-1 & x<=-1 \\
x & -1<x<=1 \\
1 & x>1
\end{cases}
\end{equation}
\]
其他部分：
\[\begin{array}{l}
e_1(n+1)=1-y(n) \\
\frac{k(z)}{g(z)}=z\left[G_{h} \frac{1}{s}\right]=\frac{T_p}{z-1} \\
k(n+1, m+1)=k(n+1, m)+T p g(n+1, m) \\
\frac{h(z)}{k(z)}=z\left[G_{h} \frac{1}{s+1}\right]=\frac{1-e^{-T_p}}{z-e^{-T_p}} \\
h(n+1, m+1)=e^{-T_p} h(n+1, m)+(1-e^{-T_p}) k(n+1, m)
\end{array}
\]
对应代码：
Private Sub Command1_Click()
Dim e1(102) As Single
Dim e2(102) As Single
Dim u(102) As Single
Dim x(102) As Single
Dim y(102) As Single
Dim k(102, 11) As Single
Dim g(102, 11) As Single
Dim h(102, 11) As Single
tp = 0.01
For n = 1 To 100
e1(n + 1) = 1 - y(n)
e2(n + 1) = e1(n + 1) + 0.21 * e1(n) - 0.632 * e2(n)
u(n + 1) = e2(n)
If u(n + 1) <= 1 And u(n + 1) >= -1 Then
x(n + 1) = u(n + 1)
End If
If u(n + 1) > 1 Then
x(n + 1) = 1
End If
If u(n + 1) < -1 Then
x(n + 1) = -1
End If
For m = 1 To 10
g(n + 1, m) = x(n + 1)
k(n + 1, m + 1) = k(n + 1, m) + tp * g(n + 1, m)
h(n + 1, m + 1) = Exp(-tp) * h(n + 1, m) + (1 - Exp(-tp)) * k(n + 1, m)
Next
y(n + 1) = h(n + 1, 10)
Next
Text1.Text = y(101)
End Sub
倒数第二节课堂练习
- 产生100个[0,1]区间的随机数,并统计大于0.6的数的个数(RND\SQR\LOG)
- 产生1000个[5,8]区间的随机数,并统计 [6,7]区间数的个数
- L1:20(0,0.1) 正态(6△原则)(x=u+σy) L2:30(0,0.2) 正态 求L1+L2在[50.1,50.2]之间的概率(仿1000次)
\[\begin{array}{l} y_{1}=\sqrt{-2 \ln r_{1}} \cos 2 \pi r_{2} \\ y_{2}=\sqrt{-2 \ln r_{1}} \sin 2 \pi r_{2} \end{array} \]
- 产生500个λ=5的指数分布值,显示500个的均值
- 产生100个λ =3.5的泊松分布值,显示100个的均值
- 产生15000个P=0.2 m=25的二项分布值,显示值为1的个数
代码：
Private Sub Command1_Click()
n = 0
'第1题
For i = 1 To 100
x = Rnd()
If x > 0.6 Then
n = n + 1
End If
Next
Text1.Text = n
End Sub
Private Sub Command2_Click()
'第2题
n = 0
For i = 1 To 1000
x = 5 + 3 * Rnd()
If x >= 6 And x <= 7 Then
n = n + 1
End If
Next
Text2.Text = n
End Sub
Private Sub Command3_Click()
'第3题
n = 0
For i = 1 To 1000
Y1 = Sqr(-2 * Log(Rnd())) * Cos(2 * pi * Rnd())
Y2 = Sqr(-2 * Log(Rnd())) * Sin(2 * pi * Rnd())
X1 = 20.05 + (0.1 / 6) * Y1
X2 = 30.1 + (0.2 / 6) * Y2
If (X1 + X2) > 50.1 And (X1 + X2) < 50.2 Then
n = n + 1
End If
Next
Text3.Text = n / 1000
End Sub
Private Sub Command4_Click()
'第4题
Sum = 0
For i = 1 To 500
x = -1 * (1 / 5) * Log(1 - Rnd())
Sum = Sum + x
Next
Text4.Text = Sum / 500
End Sub
Private Sub Command5_Click()
'第5题
Sum = 0
For i = 1 To 100
r = 1
n = 0
Do While (r >= Exp(-3.5))
n = n + 1
r1 = Rnd()
r = r * r1
Loop
y = n - 1
Sum = Sum + y
Next
Text5.Text = Sum / 100
End Sub
Private Sub Command6_Click()
'第6题
Sum = 0
p = 0.2
m = 25
For cnt = 1 To 15000
i = 0
n = 0
For cnt1 = 1 To 25
r = Rnd()
If r > (1 - p) Then
n = n + 1
End If
i = i + 1
Next
y = n
If (1 <> y) Then
Sum = Sum + 1
End If
Next
Text6.Text = 15000 - Sum
End Sub
运行结果：
