计算机仿真的过程与方法

《通信系统仿真》实验报告
姓名杨利刚班级A0811 实验室203 组号28 学号28 实验日期
实验名称实验一计算机仿真的过程与方法实验成绩教师签字
一、实验目的
1、掌握计算机仿真的一般过程
2、掌握Matlab编程仿真的基本方法
3、掌握动态系统模型的状态方程求解方法
4、掌握基于概率模型的蒙特卡罗方法
二、实验原理
1、计算机仿真的一般过程
通信系统的计算机仿真就是根据物理系统的运行原理建立相应的数学描述，并进行计算机数值求解的过程。

系统的数学描述称为系统数学模型或仿真模型。

为了对系统数学模型进行计算机数值分析，还需要将数学模型以某种计算机语言表达出来，然后进行调试、运行，最后得出数值结果。

用计算机语言重新表达的数学模型称为系统的计算机仿真模型。

根据物理模型的不同特点、原理以及不同的系统仿真目标所得出的数学模型和相应求解算法也不尽相同。

通信系统的计算机仿真过程往往是多种形式数学模型和各种算法综合的数值计算过程。

对仿真模型和仿真结果的检验是仿真数据有效性的保证。

通常的验证方法是证伪，而不是证实。

通过模型的相互比较就能够查找出错误根源，进而改进和修正模型。

2、基于动态系统模型的状态方程求解方法
动态系统，就是有记忆系统的数学描述是状态方程。

对动态系统建模，就是根据研究对象的物理模型找出相应的状态方程的过程。

所谓对动态系统的仿真，就是利用计算机来对所得出的状态方程进行数值求解的过程。

3、基于概率模型的蒙特卡罗方法
对于不确定系统，使用基于概率模型的蒙特卡罗方法。

蒙特卡罗方法是一种基于随机试验和统计计算的数值方法，也称计算机随机模拟方法或统计模拟方法。

蒙特卡罗方法的数学基础是概率论中的大数定理和中心极限定理。

大数定理指出，随着独立随机试验次数增加，试验统计事件出现的概率将接近于该统计事件的概率。

蒙特卡罗方法的基本思想：当所求解问题是某种随机事件出现的概率，或某个随机变量的期望值时，通过某种实验的方法，以这种事件出现的概率来估计该随机事件的概率，或者得出这个随机变量的某些数字特征，并将其作为问题的解。

如所求解的问题不是一个随机事件问题，那么可以通过数学分析方法找出与之等价的随机事件模型，然后再利用蒙特卡罗方法去求解。

在建模和仿真中，应用蒙特卡罗方法主要有两部分工作：
1）用蒙特卡罗方法模拟某一过程时，产生所需要的各种概率分布的随机变量。

2）用统计方法把模型的数字特征估计出来，从而得到问题的数值解，即仿真结果。

三、实验内容
1、在例题1.1中，如果考虑落体受到空气阻力，且阻力与下落速度成正比，试修改数学模型和相应的仿真程序。

在考虑阻力的情况下，在相同高度同时下落的质量不同的物体仍然同时落地吗？请通过仿真验证并解释之。

实验代码：
% ch1example1prg1.m
g=9.8; % 重力加速度
k0=-1; %阻力系数
m=1; %物体质量
v=0; % 设定初始速度条件
s=0; % 设定初始位移条件
t=0; % 设定起始时间
dt=0.01; % 设置计算步长
N=200; % 设置仿真递推次数. 仿真时间等于N与dt的乘积
for k=1:N
a=g+k0*v/m; %加速的计算
v=v+a*dt; % 计算新时刻的速度
s(k+1)=s(k)+v*dt; % 新位移
t(k+1)=t(k)+dt; % 时间更新
end
% 理论计算, 以便与仿真结果对照
t_theory=0:0.01:N*dt; % 设置解析计算的时间点
v_theory=g*t_theory; % 解析计算的瞬时速度
s_theory=1/2*g*t_theory.^2; % 解析计算的瞬时位移
% 作图: 仿真结果与解析结果对比
t=0:dt:N*dt;
plot(t,s,'o', t_theory,s_theory, '+');
xlabel('时间 t'); ylabel('位移 s');
legend('仿真结果','理论结果');
运行结果：
0.2
0.4
0.6
0.8
1 1.2
1.4
1.6
1.8
2
024681012
14161820时间 t
位移 s
仿真结果理论结果
实验结论：
通过实验数据的验证，从相同高度下落的物体，在考虑阻力的情况下，质量大的物体先落地。

2、在例题1.1中，除了考虑落体受到空气阻力，如果再考虑空气对物体的浮力，那么如何进一步修改例题1.1中的数学模型和相应的仿真程序呢？请通过仿真验证你的模型并解释之。

实验代码：
% ch1example1prg1.m
clear all ,close all ,clc
g=9.8; % 重力加速度 k0=-1; %阻力系数 m=1; %物体质量 V=0.1 %物体体积 p=1.29%空气密度
v=0; % 设定初始速度条件 s=0; % 设定初始位移条件 t=0; % 设定起始时间 dt=0.01; % 设置计算步长
N=200; % 设置仿真递推次数. 仿真时间等于N 与dt 的乘积 for k=1:N
a=g+k0*v/m-p*V*g/m; %加速的计算 v=v+a*dt; % 计算新时刻的速度
s(k+1)=s(k)+v*dt; % 新位移 t(k+1)=t(k)+dt; % 时间更新 end
% 理论计算, 以便与仿真结果对照
t_theory=0:0.01:N*dt; % 设置解析计算的时间点 v_theory=g*t_theory; % 解析计算的瞬时速度 s_theory=1/2*g*t_theory.^2; % 解析计算的瞬时位移 % 作图: 仿真结果与解析结果对比 t=0:dt:N*dt; figure(1);
plot(t,s,'o', t_theory,s_theory, '+'); xlabel('时间 t'); ylabel('位移 s'); legend('仿真结果','理论结果'); 运行结果：
0.2
0.4
0.6
0.8
1 1.2
1.4
1.6
1.8
2
024681012
14161820时间 t
位移 s
仿真结果理论结果
实验结论：
通过实验数据的验证，从相同高度下落的物体，在考虑阻力和浮力的情况下，相同质量的物体，体积大的后落地。

3、请用解析法求出抛物线y=2-x2与y=0所围区域的面积表达式，然后采用蒙特卡罗方法仿真计算该区域的近似面积，并从仿真实验中观察随机试验次数与仿真结果精度之间的大致关系。

实验代码：
% ch1example3prg1.m
clear all,close all,clc
x=-sqrt(2):0.01:sqrt(2);
y=2-x.*x;% 计算半径为1的圆周上的点，以便作出圆周观察
m=0; % 在圆内在落点计数器
N=1000; % 设置试验次数
x1=2*sqrt(2)*rand(N,1)-sqrt(2);% 产生均匀分布于[-1, +1]直接的两个独立随机数x1,y1
y1=2*rand(N,1);
for n=1:N % 循环进行重复试验并统计
p1=x1(1:n);
q1=y1(1:n);
if(2-x1(n)*x1(n))<y1(n) % 计算落点到坐标原点的距离,判别落点是否在圆内
m=m+1; % 如果落入圆中,计数器加1 end
%plot(p1,q1,'.',x,y,'-k',[-2 -2 2 2 -2],[-2 2 2 -2 -2],'-k'); plot(p1,q1,'.',x,y,'-k',[-sqrt(2) sqrt(2) sqrt(2) -sqrt(2) -sqrt(2)],[0 0 2 2 0],'-k');
axis equal; % 坐标纵横比例相同
axis([-2 2 -0.7 2.2]); % 固定坐标范围
text(-1,-0.2,['试验总次数 n=',num2str(n)]);% 显示试验结果
text(-1,-0.4,['落入圆中数 m=',num2str(m)]);
text(-1,-0.6,['近似圆面积 S_c=',num2str(m/n*4*sqrt(2))]);
set(gcf,'DoubleBuffer','on'); % 双缓冲避免作图闪烁
drawnow; % 显示结果
end
运行结果：
-2
-1.5
-1
-0.5
0.5
1
1.5
2
-0.5
0.5
1
1.5
2
试验总次数 n=1000
落入圆中数 m=330近似圆面积 S c =1.8668
实验结论：
通过对实验数据的综合分析，可以得出实验次数越多，结果将越精确。

四、实验报告
通过分析实验结果可知：
从相同高度下落的物体，在考虑阻力的情况下，质量不同的物体不会同时落地。

从相同高度下落的物体，在考虑阻力和浮力的情况下，相同质量的物体，不同体积落地时间先后不一样，体积小的落地。

通过对实验数据的综合分析，可以得出实验次数越多，结果将越精确。