计算机仿真试题及答案

南京林业大学研究生试卷课程号 33351 课程名称 计算机仿真 20 16 ~20 17 学年第 1 学期 一、（20分）某系统的传递函数为)4)(3)(2)(14(143)(2++++++=S S S S S S S G 试用级联法、串联法和并联法中的两种画出系统模拟结构图，并列出相应的状态方程和输出方程。

自拟仿真，比较两者结果。

答：选用级联法和并联法写出传递函数，公式如下：)4(1)3(43)2()31()41()1(/6/4122/4113/4371/41/1/43)4)(3)(2)(14(143)(4324322+∙+∙++∙++=++++++=++++++=s s s s s s s s s s s s s s s s s s s s G分别画出级联法、并联法的系统模拟结构图及其状态方程和输出方程，图式如下：432342312432114143'''641224113437'x x x y x x x x x x x x x x u x ++====----=4434321321211443'33543'243'41'x y x x x x x x u x x x u x x u x =-=--+=-+=-=使用龙格—库塔法程序进行仿真，令u=1，系统在t=0-2s 之间（仿真步长0.1s ）则输出y截图如下：比较后发现两者仿真结果接近。

二、（20分）已知系统的状态方程及输出方程 11231021232031233012332022022022x x x x x xx x x x x x x x x y x x x =+-=⎧⎪=-+=⎪⎨=+-+=⎪⎪=+-⎩拟用四阶龙格-库塔法求系统在t=0-2s 之间的输出值（仿真步长取0.01s 或0.1s ），试编写仿真程序，对仿真结果截图。

（单号C++，双号VB ）答：用c++进行四阶龙格-库塔法，系统在t=0-2s间的（仿真步长0.01）编程：#include<iostream>#include<cstring>using namespace std;float x1;float x2;float x3;float x4;float t;float f1(float x1,float x2,float x3,float x4,float t){float f1;return f1=3*x1+2*x2-x3;}float f2(float x1,float x2,float x3,float x4,float t){float f2;return f2=2*x1-x2+2*x3;}float f3(float x1,float x2,float x3,float x4,float t){float f3;return f3=x1+2*x2-x3+2;}float f4(float x1,float x2,float x3,float x4,float t){float f4;return f4=0;}int main(){float k11,k12,k13,k14,k21,k22,k23,k24,k31,k32,k33,k34,k41,k42,k43,k44;float t=0,h=0.01;double y;float x1=0,x2=0,x3=0,x4=0;for(int i=1;i<200;i++){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;}y=x1+2*x2-2*x3;cout<<"y="<<y<<endl;return 0;}仿真结果截图如下：三、（20分）已知单位反馈系统的开环传递函数为()()()()34150+++=sssss G试用Z 域离散相似法和时域离散相似法进行仿真，编写仿真程序，对仿真结果截图。

（输出初态为0，输入为3，仿真步长为0.02S,100-200个周期）（单号VB ，双号C++）答：先将函数()()()()34150+++=s s s s s G 化成为：50(1)1+150()(4)(3)s+4s+3s s G s s s s s +==∙∙++ 画出系统结构图：使用z 域离散相似法计算得：3/))*3exp(1(*50*)()*3exp(*)()1(4/))*4exp(3(*)()1()*4exp(*)()1(*3)()1(t n x t n y n y t n w n w t n x n x tn w n w --+-=+-+-++-=++=+ 仿真程序如下：Private Sub Command1_Click() Dim n As Single Dim t As SingleDim w(200) As Single Dim x(200) As Single Dim z(200) As Single Dim y(200) As Single Dim u(200) As Single t = 0.02 y(0) = 0For n = 0 To 199w(n + 1) = w(n) + t * 3x(n + 1) = x(n) * Exp(-4 * t) + w(n + 1) - w(n) * (3 + Exp(-4 * t)) / 4 y(n + 1) = y(n) * Exp(-3 * t) + x(n) * 50 * (1 - Exp(-3 * t)) / 3 NextText1.Text = y(200) End Sub仿真结果截图如下：使用时域相似法计算得：3/))*3exp(1(*)(*50)(*)*3exp()1())*4exp(*25.025.0(*)3)(()*4exp(*)()1(3*)()1(T n x n y T n y T n w T n x n x T n w n w --+-=+--++-=++=+ 仿真程序如下：Private Sub Command1_Click() Dim w(2000) As Single Dim x(2000) As Single Dim y(2000) As Single Dim T As Single Dim n As SingleT = 0.02y(0) = 0For n = 0 To 199w(n + 1) = w(n) + T * 3x(n + 1) = x(n) * Exp(-4 * T) + (w(n) + 3) * (0.25 - 0.25 * Exp(-4 * T))y(n + 1) = Exp(-3 * T) * y(n) + 50 * x(n) * (1 - Exp(-3 * T)) / 3NextText1.Text = y(200)End Sub仿真结果截图如下：四、（20分）投掷硬币实验，至少投掷多少次，才能使正面朝上的频率在0.4-0.6之间的概率不小于0.9 ？试编写仿真程序，对仿真结果截图。

（仿真8000-12000次）答：（仿真12000次）程序如下：Private Sub Command1_Click()Dim a As SingleDim i As IntegerDim j As IntegerDim t As IntegerDim s As SingleFor j = 1 To 12000n = 0For i = 1 To 25r = Rnd()If r >= 0.5 Thenn = n + 1End IfNextIf 0.4 <= n / 25 <= 0.6 Thenm = m + 1End IfIf m / 25 >= 0.9 Thens = s + 1End IfNextText1.Text = sEnd Sub仿真结果截图如下：五、（20分）某自选市场，平均每分钟有6位顾客到来，每位顾客购物时间服从（4.2，7.2）分钟均布（不包括交费时间）。