实验一面向微分方程的数值积分法仿真 2.docx
《实验一面向微分方程的数值积分法仿真 2.docx》由会员分享，可在线阅读，更多相关《实验一面向微分方程的数值积分法仿真 2.docx（14页珍藏版）》请在冰豆网上搜索。
实验一面向微分方程的数值积分法仿真2
实验一面向微分方程的数值积分法仿真
一、实验目的
1．掌握数值积分法的基本概念、原理及应用；
2．用龙格-库塔法解算微分方程，增加编写仿真程序的能力；
3．分析数值积分算法的计算步长与计算精度、速度、稳定性的关系；
4.对数值算法中的“病态问题”进行研究。
二、实验内容
1、已知系统微分方程及初值条件
取步长
，试分别用欧拉方程法和RK4法求
时的
值，并将求得的值与解析解
比较（将三个解绘于同一坐标中，且用数值进行比较），说明造成差异的原因。
（①编程完成；②选用MATLABode函数完成。
）
程序代码如下：
t0=0;
tf=2;
h=0.1;
y1=1;
y2=1;
y3=1;
t1=0;
t2=0;
t3=0
n=round（tf-t0）/h;
fori=1:
n
y1（i+1）=y1（i）+h*（2*h+y1（i））;
t1=[t1,t1（i）+h];
end
fori=1:
n
k1=y2（i）+t2（i）;
k2=y2（i）+h*k1/2+t2（i）+h/2;
k3=y2（i）+h*k2/2+t2（i）+h/2;
k4=y2（i）+h*k3+t2（i）+h;
y2（i+1）=y2（i）+h*（k1+2*k2+2*k3+k4）/6;
t2=[t2,t2（i）+h];
end
fori=1:
n
y3（i+1）=2*exp（t3（i））-t3（i）-1;
t3=[t3,t3（i）+h];
end
plot（t1,y1,'r',t2,y2,'g',t3,y3,'k'）
实验结果如下；
分析：
红线为用欧拉法得到的结果，绿线为用四阶龙格—库塔法得到的结果，蓝线为根据解析方程得到的结果。
其差异原因主要有两个：
1、二者的方法不同，欧拉法是根据一阶微分方程计算得到的，龙格—库塔法是根据四阶微分方程得到的；2、由于步长取为0.1，所以得到的图像与解析解之间存在差异，若将步长取小，则得到的解将更靠近解析解。
2、已知系统的传递函数为
在单位阶跃输入下，系统响应的解析解为
试分别用欧拉方程法和RK4法对系统进行仿真（编程完成）：
1）比较两种数值积分解与解析解得逼近程度；（绘图）
程序代码如下：
a=[010;001;-22.06-27-10];
b=[0;0;1];
c=[40.600];
X1=[0;0;0];t=0;Y1=0;
X=0;
u=1;
Y2=0;Y3=0;
X2=[0;0;0];
x=0;
h=0.1;
t0=0;
tf=2;
t1=0;t2=0;t3=0;
N=（tf-t0）/h;
fori=1:
N
k1=a*X1+b;
k2=b+a*（h*k1/2+X1）;
k3=b+a*（h*k2/2+X1）;
k4=b+a*（h*k3+X1）;
X1=X1+h*（k1+2*k2+2*k3+k4）/6;
Y1=[Y1,c*X1];
t1=[t1,t1（i）+h];
end
fori=1:
N
x=X2（:
i）+h*（a*X2（:
i）+b*u）;
y=c*x;
X2=[X2,x];
Y2=[Y2,y];
t2=[t2,t2（i）+h];
end
fori=1:
N
y=1.84-4.95*i*exp（-1.88*i）-1.5*exp（-1.88*i）-0.34*exp（-6.24*i）;
Y3=[Y3,y];
t3=[t3,t3（i）+h];
end
plot（t1,Y1,'r',t2,Y2,'g',t3,Y3,'b'）
当h=0.01时的结果
当h=0.01时的结果
分析：
这是我得到的结果，发现两个方法得到的结果与实际结果都有较大差距，当是龙格—库塔法更接近实际的结果。
2）改变步长，分析步长对数值解精度的影响；
改变步长后，发现只是两根仿真得到的曲线靠近了，但是与实际曲线仍然是差距很大，这是经过仔细的检查和讨论我觉得程序还是对的，不知道错在哪里了。
3）不断加大步长，分析计算稳定性的变化。
当取h=0.5时，得到的结果：
加大步长后结果得到的结果不稳定，不能够很好的对系统进行仿真，另外，由于系统步长选择偏大，根据解析解得到的结果也与实际值有了一定的差距，但是如果步长取得不一样又无法比较。
3、求下图所示系统的阶跃响应
的数值解。
（，
，
，
，
）分析
、
对系统响应的影响。
（①编程用RK4求解；②Simulink）
程序代码如下：
k=1;
a=conv（[100],conv（[0.251],[0.251]））
b=[2*kk]
X0=[0000]
v=1;n0=4;tf=10;h0=0.25;r=1;
V=v;
n=n0;
T0=0;
Tf=tf;
h=h0;
R=r;
b=b/a
（1）;a=a/a
（1）;A=a（2:
n+1）;
A=[rot90（rot90（eye（n-1,n）））;-fliplr（A）];
B=[zeros（1,n-1）,1]';
m1=length（b）;
C=[fliplr（b）,zeros（1,n-m1）];
Ab=A-B*C*V;
X=X0';y=0;t=T0;
N=round（Tf-T0）/h;
fori=1:
N
k1=Ab*X+B*R;
k2=Ab*（X+h*k1/2）+B*R;
k3=Ab*（X+h*k2/2）+B*R;
k4=Ab*（X+h*k3）+B*R;
X=X+h*（k1+2*k2+2*k3+k4）/6;
y=[y,C*X];
t=[t,t（i）+h];
end
[t',y']
plot（t,y）
得到结果如下：
K=0.5,v=1时的结果如下：
K=1，v=1时，
K=1，V=0.5时
K=1,v=5时，
K=2，v=1时
分析：
当k取值增大，v值不变时，系统输出的波头增多，而且也变陡，稳态精度降低，当k增加到一定程度时系统便发散了（即不稳定了）。
当v值增大，k值不变时，波头也是变多变陡，当v值增大到一定程度时系统便不稳定了。
4、已知系统状态状态方程为
采用RK4法，步长分别取
，求系统的零输入响应，并绘图分析各状态变量的响应状态及产生的原因。
（提示：
病态系统）
程序代码如下：
a=[-2119-20;19-2120;40-40-40];
x=[1;0;-1];X=x;t=0;
t0=0;tf=2;h=0.01;
n=round（tf-t0）/h;
fori=1:
n
x=X（:
i）+h*a*X（:
i）;
X=[X,x];
t=[t,t（i）+h];
end
b=X（1,:
）;c=X（2,:
）;d=X（3,:
）;
plot（t,b,'r',t,c,'g',t,d,'b'）
当h=0.01时得到的结果
当h=0.02时得到的结果
当h=0.04时得到的结果
分析：
如图，当h=0.01时，在t=0.2s以后系统输出便趋于平稳，当取h=0.02时，系统输出振荡剧烈，趋于稳定的时间也变长，当取h=0.04后，系统输出呈发散振荡形式。
当h=0.06后系统仍然是发散的，即当h的取值改变时，原先稳定的系统变得不稳定了，这便是病态系统。
但是这个结果与书上的不同。
三、实验报告要求
记录完成实验内容所采取的步骤、方法和结果。
回答思考题。
四、预习要求及思考题
要求实验前，必须预习实验知识点，按实验内容的要求的确定仿真方案，完成程序设计，以便在实验时进行调试、分析。
思考题：
1、在进行仿真计算时，是否选用的数值积分法的阶次越高越好？
答：
阶次不是越高越好。
（1）阶次越高计算公式也越复杂，每一步需要计算的次数也更多，需要的时间也更长。
（2）每一阶的龙格库塔法都有对应的稳定区域，当hλ（λ为系统的特征方程根）超过稳定区域时，即使阶次高得到的结果也未必是稳定的。
当hλ得值接近稳定边界时，误差也会增大。
2、选用数值积分法进行仿真的原则。
答：
（1）精度仿真结果的精度主要受三项误差的影响：
1）截断误差：
由算法本身的精度阶次所决定。
2）舍入误差：
由计算机字长决定。
3）累积误差：
由以上两项误差随计算时间长短累积情况决定。
（2）计算速度计算速度取决于所用的数值方法和计算步长。
（3）稳定性数值稳定性主要与计算步长h有关，不同的数值方法对h都有不同的稳定性限制范围，且与被仿真对象的时间常数也有关系。
升级会员 