2022-10-31 10:27阅读: 81评论: 0推荐: 0

【计算机仿真技术】代码记录

计算机仿真技术代码记录

第一节课课程练习

报童问题：报童每卖一份报可赚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^{'''} + 0.75 x^{''} + x^{'} + 5 x = 5

$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}^{'} = - x_{1} + u \\ x_{2}^{'} = - 3 x_{2} + u \\ x_{3}^{'} = - 4 \times 3 + u \\ y = x 1 + x_{2} - 2 \times 3 \\ u = 2 t \end{array}

$\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}}^{'} = x_{2} + 1 \\ {x_{2}}^{'} = x_{1} - 2 x_{2} - 1 \\ x_{1} (0) = x_{2} (0) = 0 \end{array}

$\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^{'} = 1 - y \\ x^{'} = - 0.1 x + 1 - y + 0.5 w \\ z^{'} = - 2 z + 2 x \\ y^{'} = - 10 y + 10 z \end{array}

$\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{array}{l} D(z) = \frac{E_2(z)}{E_1(z)}=\frac{z+0.21}{z+0.632} \end{array}$

非线性部分：

\begin{matrix} (1) & f (x) = {\begin{cases} - 1 & x <= - 1 \\ x & - 1 < x <= 1 \\ 1 & x > 1 \end{cases} \end{matrix}

$\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 [G_{h} \frac{1}{s}] = \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 [G_{h} \frac{1}{s + 1}] = \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}

$\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