欧拉法仿真

仿真作业
1.普通欧拉公式法
clear;
x1=0;x2=0;
M=1.0;K=1.0;B=0.5;
U=1.0;
ST=20; T=0.1
NP=ST/T;
for i=1:1:NP
F1=x2;
F2=(-K/M)*x1-(B/M)*x2+(1/M)*U;
x1=x1+T*F1;
x2=x2+T*F2;
y(i)=x1;
end
plot((1:NP)*T,y,'m-')
hold on
2.梯形公式法
clear;
x10=0;x20=0;
M=1.0;K=1.0;B=0.5;
U=1.0;
ST=50; T=0.1
NP=ST/T;
for i=1:NP
F10=x20;
F20=(-K/M)*x10-(B/M)*x20+(1/M)*U;
x11=x10+T*F10;
x21=x20+T*F20;
F11=x21;
F21=-K*x11/M-B*x21/M+U/M;
x12=x10+T*(F10+F11)/2;
x22=x20+T*(F20+F21)/2;
x10=x12;
x20=x22;
y(i)=x12;
end
plot((1:NP)*T,y,'r-')
hold on
3.四阶龙格－库塔法（R-K法）
clear;
X10=0;X20=0;
M=1;K=1;B=0.5;U=1;
ST=30;T=0.1;
NP=ST/T;
for i=1:NP
F10=X20;
F20=-K/M*X10-B/M*X20+U/M;
X11=X10+T/2*F10;
X21=X20+T/2*F20;
F11=X21;
F21=-K/M*X11-B/M*X21+U/M;
X12=X10+T/2*F11;
X22=X20+T/2*F21;
F12=X22
F22=-K/M*X12-B/M*X22+U/M;
X13=X10+T*F12;
X23=X20+T*F22;
F13=X22;
F23=-K/M*X13-B/M*X23+U/M;
F1=(F10+2*F11+2*F12+F13)/6;
F2=(F20+2*F21+2*F22+F23)/6;
X14=X10+T*F1;
X24=X20+T*F2;
X10=X14;
X20=X24;
Y(i)=X14;
end
plot((1:NP)*T,Y,'r-')
hold on
4.1隐式欧拉公式法
clear;
X10=0;X20=0;
M=1;K=1;B=0.5;U=1;
ST=30;T=10;
NP=ST/T;
for i=1:NP
X11=X10+T*(X20-T*K*X10/M+T*U/M)/(1+T*T*K/M+T*B/M);
X21=(X20-T*K*X10/M+T*U/M)/(1+T*T*K/M+T*B/M);
X10=X11;
X20=X21;
Y(i)=X11;
end
plot((1:NP)*T,Y,'m-')
hold on
当我们继续增大仿真时间与仿真步距时：例如ST=30000,T=100,1000,10000时，观测到的仿真曲线都不发散。

4.2局部隐式欧拉公式
clear;
X10=0;X20=0;
M=1;K=1;B=1;U=1;
ST=30;T=0.1;
NP=ST/T;
for i=1:NP
X11=X10+T*X20;
X21=(X20-T*K*X11/M+T*U/M)/(1+T*B/M);
X10=X11;
X20=X21;
Y(i)=X11;
end
plot((1:NP)*T,Y,'M-')
hold on
分析：
1.显然，由仿真曲线图可以看到：
(1)使用普通欧拉公式法时，仿真曲线随仿真步距T的变化，当T=0.1时曲线很平
滑，当T=0.3，0.4时曲线失真，当T=0.5时，曲线发散。

(2)使用梯形公式法时，仿真曲线随仿真步距T的变化，当T=0.1时曲线很平滑，
T=0.9，1.3时曲线失真，当T=1.5时曲线就发散了。

(3)使用四阶龙格－库塔法时，仿真曲线随仿真步距T的变化，当T=0.1时曲线很
平滑，T=0.8，1.6，2.1时曲线失真，当T=2.7时曲线就发散了。

(4)使用隐式欧拉公式法时，仿真曲线随仿真步距T的变化，当T=0.1时曲线很平
滑，T=1，2，5，10时曲线失真但不发散，我进一步设置ST＝30000，然后分
别让T=100,300,1000,10000从得到的曲线看出失真但不发散。

(5)使用局部隐式欧拉公式法时，仿真曲线随仿真步距T的变化，当T=0.1时曲线
很平滑，当T=0.9,1.2,1.8时曲线失真，当T=2.6时，曲线发散。

综上可以得出结论：从稳定性考虑，使用隐式欧拉公式时，尽管T取得非常大，曲线不会出现发散的现象，其次稳定性较强的是四阶龙格库塔法和梯形公式法，最弱的是局部隐式欧拉公式法和普通欧拉公式法。

2.我在相应的程序最后加上一句“Y(i)=****”使输出数值显示在命令窗口。

（1）对于普通欧拉公式法，首先把T设定很小的值比如：T=0.01，输出为：1.0005，以此值作为真值；增大T使T=0.3，输出为：0.9848。

（2）对于梯形公式法，首先设定T=0.01，输出为：1.0005；增大T使T=0.3，输出为：1.0007；继续增大T使T=1.2，输出为：0.9956。

（3）对于四阶龙格－库塔法，首先设定T=0.01，输出为：1.0005；增大T使T=0.3，输出为：1.0003；继续增大T使T=1.2，输出为：1.0013。

（4）对于隐式欧拉公式法，首先设定T=0.01，输出为：1.0005；增大T使T=0.3，输出为：1.0000。

继续增大T使T=1.2，输出为：1.0000。

（5）对于局部隐式欧拉公式法，首先设定T=0.01，输出为：1.0005；增大T使T=0.3，输出为：1.0009。

综上，可以得出结论，从精确性考虑,使用四阶龙格－库塔法时精确度最高，即随T增大时误差最小，其次精确度较高的是梯形公式法、隐式欧拉公式法和局部隐式欧拉公式法，最差的是普通欧拉公式法。