计算机仿真技术实验报告-实验三.doc
《计算机仿真技术实验报告-实验三.doc》由会员分享，可在线阅读，更多相关《计算机仿真技术实验报告-实验三.doc（13页珍藏版）》请在冰豆网上搜索。
计算机仿真技术实验报告
实验三利用数值积分算法的仿真实验
实验三利用数值积分算法的仿真实验
一．实验目的
1）熟悉MATLAB的工作环境；
2）掌握MATLAB的.M文件编写规则，并在命令窗口调试和运行程序；
3）掌握利用欧拉法、梯形法、二阶显式Adams法及四阶龙格库塔法构建系统仿真模型的方法,并对仿真结果进行分析。
二．实验内容
系统电路如图2.1所示。
电路元件参数：
直流电压源，电阻，电感，电容。
电路元件初始值：
电感电流，电容电压。
系统输出量为电容电压。
连续系统输出响应的解析解为：
（2-1）
其中，，。
三、要求
1） 利用欧拉法、梯形法、二阶显式Adams法及显式四阶Runge-Kutta法构建系统仿真模型，并求出离散系统的输出量响应曲线；
2） 对比分析利用欧拉法、梯形法、二阶显式Adams法及显式四阶Runge-Kutta法构建系统仿真模型的仿真精度与模型运行的稳定性问题；
3） 分别编写欧拉法、梯形法、二阶显式Adams法及显式四阶Runge-Kutta法的.m函数文件，并存入磁盘中。
.m函数文件要求输入参数为系统状态方程的系数矩阵、仿真时间及仿真步长。
编写.m命令文件，在该命令文件中调用已经编写完成的上述.m函数文件，完成仿真实验；
4）subplot和plot函数将输出结果画在同一个窗口中，每个子图加上对应的标题。
四.实验原理
（1）连续系统解析解
连续系统输出响应的解析解为：
其中，，
（2）原系统的传递函数
根据所示电路图，我们利用电路原理建立系统的传递函数模型，根据系统的传递函数是在零初始条件下输出量的拉普拉斯变换与输入量的拉普拉斯变换之比，可得该系统的传递函数：
（3）系统的仿真模型
在连续系统的数字仿真算法中,较常用的有欧拉法、梯形法、二阶显式Adams法及显式四阶Runge-Kutta法等。
欧拉法、梯形法和二阶显式Adams法是利用离散相似原理构造的仿真算法，而显式四阶Runge-Kutta法是利用Taylor级数匹配原理构造的仿真算法。
对于线性系统，其状态方程表达式为:
其中：
是系统的n维状态向量
是系统的m维输入向量
是系统的r维输出向量
A为阶参数矩阵，又称动态矩阵，B为阶输入矩阵，C为阶输出矩阵，D为阶交联矩阵。
根据图所示电路，系统状态方程模型:
式中，状态变量，输出变量，系数矩阵为：
，，。
（1）欧拉法
利用前向欧拉法构建线性系统的仿真模型为：
式中，为积分步长，为单位矩阵。
利用后向欧拉法构建线性系统的仿真模型为：
对于前向欧拉法，系数矩阵为：
，，，D=0。
对于后向欧拉法，系数矩阵为：
，，。
（2）梯形法
利用梯形法构建线性系统的仿真模型为：
对图所示的系统，利用梯形法构造的系统差分方程具有形式：
其系数矩阵为：
，，，，D=0。
（3）二阶显式Adams法
利用二阶显式Adams法构建线性系统的仿真模型为：
式中：
二阶显式Adams法为多步计算方法，利用多步计算方法对系统进行仿真时，需要与之具有相同计算精度的单步计算方法辅助计算。
二阶显式Adams法的计算精度为二阶，可以采用梯形法或改进的Euler法等辅助计算。
利用改进的Euler法构建线性系统的仿真模型为：
其中，。
由式计算出和后，便可以转入由二阶显式Adams法构造的离散系统模型计算，即系统差分方程。
其计算方程为：
（）
（4）显式四阶Runge-Kutta法
利用显式四阶Runge-Kutta法构建线性系统的仿真模型为：
五．实验过程
1.实验程序
（1）前向欧拉法
function[]=RLC（R,L,C,U,t,h）
R=10;
L=0.01;
C=1.0e-6;
U=1;
t=0.01;
h=2.0e-4;
m=fix（t/h）;
n=2;
A=[-R/L-1/L;1/C0];
B=[1/L;0];
D=[01];
E=[10;01];
%前向欧拉法%
fori=1:
1:
n
x1（1:
n,1）=0;
end
fork=1:
m
x1（1:
n,k+1）=x1（1:
n,k）+（A*x1（1:
n,k）+B）*h;
end
fork=1:
1:
m
y1（k）=D*x1（1:
n,k）;
end
%解析解%
p=R/（2*L）;
w=sqrt（1/（L*C）-（R/（2*L））^2）;
fork=1:
1:
m
y（k）=U*（1-exp（-p*（k-1）*h）*（cos（w*（k-1）*h）+sin（w*（k-1）*h）*p/w））;
end
%输出曲线%
fork=1:
1:
m
t（k）=（k-1）*h;
end
subplot（2,3,1）,plot（t,y,'g',t,y1,'r'）
legend（'y解析解,','y1前向欧拉'）
title（'前向欧拉法'）
（2）后向欧拉法
function[]=RLC（R,L,C,U,t,h）
R=10;
L=0.01;
C=1.0e-6;
U=1;
t=0.01;
h=2.0e-4;
m=fix（t/h）;
n=2;
A=[-R/L-1/L;1/C0];
B=[1/L;0];
D=[01];
E=[10;01];
%后向欧拉法%
fori=1:
1:
n
x2（1:
n,1）=0;
end
A1=inv（E-A*h）;
fork=1:
m
x2（1:
n,k+1）=A1*（x2（1:
n,k）+B*h）;
end
fork=1:
1:
m
y2（k）=D*x2（1:
n,k）;
end
%解析解%
p=R/（2*L）;
w=sqrt（1/（L*C）-（R/（2*L））^2）;
fork=1:
1:
m
y（k）=U*（1-exp（-p*（k-1）*h）*（cos（w*（k-1）*h）+sin（w*（k-1）*h）*p/w））;
end
%输出曲线%
fork=1:
1:
m
t（k）=（k-1）*h;
end
subplot（2,3,2）,plot（t,y,'g',t,y2,'r'）
legend（'y解析解,','y2后向欧拉'）
title（'后向欧拉法'）
（3）梯形法
function[]=RLC（R,L,C,U,t,h）
R=10;
L=0.01;
C=1.0e-6;
U=1;
t=0.01;
h=2.0e-4;
m=fix（t/h）;
n=2;
A=[-R/L-1/L;1/C0];
B=[1/L;0];
D=[01];
E=[10;01];
%梯形法%
fori=1:
1:
n
x3（1:
n,1）=0;
end
A2=inv（E-A*h/2）;
fork=1:
m
x3（1:
n,k+1）=A2*（x3（1:
n,k）+B*h+A*x3（1:
n,k）*h/2）;
end
fork=1:
1:
m
y3（k）=D*x3（1:
n,k）;
end
%解析解%
p=R/（2*L）;
w=sqrt（1/（L*C）-（R/（2*L））^2）;
fork=1:
1:
m
y（k）=U*（1-exp（-p*（k-1）*h）*（cos（w*（k-1）*h）+sin（w*（k-1）*h）*p/w））;
end
%输出曲线%
fork=1:
1:
m
t（k）=（k-1）*h;
end
subplot（2,3,3）,plot（t,y,'g',t,y3,'r'）
legend（'y解析解,','y3梯形法'）
title（'梯形法'）
（4）二阶显式Adams法
function[]=RLC（R,L,C,U,t,h）
R=10;
L=0.01;
C=1.0e-6;
U=1;
t=0.01;
h=2.0e-4;
m=fix（t/h）;
n=2;
A=[-R/L-1/L;1/C0];
B=[1/L;0];
D=[01];
E=[10;01];
%二阶显示Adams法%
fori=1:
1:
n
x4（1:
n,1）=0;
end
fork=1:
m
x4（1:
n,k+1）=A2*（x4（1:
n,k）+B*h+A*x4（1:
n,k）*h/2）;
end
fork=3:
m
fm1=23*（A*x4（1:
n,k）+B）;
fm2=-16*（A*x4（1:
n,k-1）+B）;
fm3=5*（A*x4（1:
n,k-2）+B）;
x4（1:
n,k+1）=x4（1:
n,k）+（fm1+fm2+fm3）*h/12;
end
fork=1:
1:
m
升级会员 