微分方程.docx
《微分方程.docx》由会员分享，可在线阅读，更多相关《微分方程.docx（13页珍藏版）》请在冰豆网上搜索。
微分方程
微分方程数值方法上机报告
BasicTheoryofOrdinaryDifferentialEquationsExperimentReport
教务处
2014年6月
上机作业
题目一:
对于初值问题u’=u,u（0）=1,在单区间【0,1】上，用Euler法，梯形法及RK方法进行计算，取步长h=0.1,0.2,0.5.
试比较，
（1）用同样的步长，三种方法中哪一个精度最好；
（2）对同一种方法以不同步长进行计算，哪一个结果最好.
所用公式：
Euler法公式：
梯形法公式：
RK法公式：
编写matlab程序：
1.Euler.m
functionUn=Euler（h）
%Un=Euler（h）
symsu
f=u;
U=zeros（1,1/h）;
U
（1）=1;
fork=1:
1/h
U（k+1）=U（k）+h*subs（f,u,U（k））;
end
Un=U;
End
2.tixingfa.m
functionUn=tixingfa（h）
%Un=tixingfa（h）
symsu
f=u;
U=zeros（1,1/h）;
U
（1）=1;
fork=1:
1/h
U（k+1）=U（k）+h*subs（f,u,U（k））;
U（k+1）=U（k）+h/2*（subs（f,u,U（k））+subs（f,u,U（k+1）））;
end
Un=U;
end
3.RK.m
functionUn=RK（h）
%Un=RK（h）
symsu
f=u;
U=zeros（1,1/h）;
U
（1）=1;
fork=1:
1/h
K1=subs（f,u,U（k））;
K2=subs（f,u,U（k）+1/2*h*K1）;
K3=subs（f,u,U（k）+1/2*h*K2）;
K4=subs（f,u,U（k）+h*K3）;
U（k+1）=U（k）+1/6*h*（K1+2*K2+2*K3+K4）;
end
Un=U;
end
图像结果显示：
编写a1.m用于画图显示Euler法不同步长计算结果比较，
t0=0:
0.001:
1;
t1=0:
0.1:
1;
t2=0:
0.2:
1;
t3=0:
0.5:
1;
U0=exp（t0）;
U1=Euler（0.1）;
U2=Euler（0.2）;
U3=Euler（0.5）;
holdon
plot（t0,U0,'k'）
plot（t1,U1,'b','Marker','x'）
plot（t2,U2,'--r','Marker','+'）
plot（t3,U3,':
r','Marker','*'）
legend（'Un=e^t','buchang0.1','buchang0.2','buchang0.5','Location','NorthWest'）
title（''Euler法不同步长计算结果比较'）
gridon
xlable（'t'）
ylable（'Un'）
holdoff
运行结果如下；图
由图可知对于Euler法步长取0.1计算结果最好。
编写a2.m用于画图显示梯形法不同步长计算结果比较，
t0=0:
0.001:
1;
t1=0:
0.1:
1;
t2=0:
0.2:
1;
t3=0:
0.5:
1;
U0=exp（t0）;
U1=tixingfa（0.1）;
U2=tixingfa（0.2）;
U3=tixingfa（0.5）;
holdon
plot（t0,U0,'k'）
plot（t1,U1,'b','Marker','x'）
plot（t2,U2,'--r','Marker','+'）
plot（t3,U3,':
r','Marker','*'）
legend（'Un=e^t','buchang0.1','buchang0.2','buchang0.5','Location','NorthWest'）
title（'梯形法不同步长计算结果比较'）
gridon
xlable（'t'）
ylable（'Un'）
holdoff
运行结果如下：
图2
图上四条线比较接近，但还是可以看到红色虚线有所偏离，为更好的比较将其早末端放大得下图3
由图可知对于梯形法步长取0.1计算结果最好。
编写a3.m用于画图显示RK法不同步长计算结果比较，
t0=0:
0.001:
1;
t1=0:
0.1:
1;
t2=0:
0.2:
1;
t3=0:
0.5:
1;
U0=exp（t0）;
U1=RK（0.1）;
U2=RK（0.2）;
U3=RK（0.5）;
holdon
plot（t0,U0,'k'）
plot（t1,U1,'b','Marker','x'）
plot（t2,U2,'--r','Marker','+'）
plot（t3,U3,':
r','Marker','*'）
legend（'Un=e^t','buchang0.1','buchang0.2','buchang0.5','Location','NorthWest'）
title（'RK法不同步长计算结果比较'）
gridon
xlable（'t'）
ylable（'Un'）
holdoff
运行结果如下：
图4
图上四条线比较接近，但还是可以看到红色虚线有所偏离，为更好的比较将其早末端放大得下图图5
由图可知对于RK法步长取0.1计算结果最好。
编写a4.m用于画图显示三种方法同一步长的精度比较，不妨取步长0.1
t0=0:
0.001:
1;
t1=0:
0.1:
1;
U0=exp（t0）;
U1=Euler（0.1）;
U2=tixingfa（0.1）;
U3=RK（0.1）;
holdon
plot（t0,U0,'k'）
plot（t1,U1,'b','Marker','x'）
plot（t1,U2,'--r','Marker','+'）
plot（t1,U3,':
r','Marker','*'）
legend（'Un=e^t','Euler','tixingfa','RK','Location','NorthWest'）
title（'取步长0.1比较三种方法的精度'）
gridon
xlable（'t'）
ylable（'Un'）
holdoff
运行结果如下：
图6
很明显，Euler法精度不及另外两种方法，但梯形法与RK法都与真实结果比较接近，先做放大处理
选择一个节点如t=0.8放大图像如下图7
由图可知还是RK法的精度较好。
结论：
综上图形结果显示，同样的步长Eluer法精度不及RK法和梯形法，RK法与梯形法比较接近但还是RK法精度好一些。
同一种方法取不同步长计算，步长越小计算结果越好。
题目二:
对于初值问题u’=u-2t/u,u（0）=1,在单区间【0,1】上，用Adams四阶预估-校正算法的PECE模式求其数值解，取步长h=0.1利用计算结果估计数值解的局部误差主项。
（真解u（t）=sqrt（1+2t））
所用公式：
Adams四阶预估-校正（PECE）公式：
编写matlab程序：
function[Un,e]=di13ti（）
%[Un,e]=di13ti（）
f0=u-2*t/u;
v=[t,u];
U=zeros（1,11）;
T=0:
0.1:
1;
f=zeros（1,11）;
h=0.1;
U
（1）=1;
f
（1）=1;
fork=1:
3
K1=subs（f0,v,[（k-1）*h,U（k）]）;
K2=subs（f0,v,[（k-1）*h+1/2*h,U（k）+1/2*h*K1]）;
K3=subs（f0,v,[（k-1）*h+1/2*h,U（k）+1/2*h*K2]）;
K4=subs（f0,v,[（k-1）*h+h,U（k）+h*K3]）;
U（k+1）=U（k）+1/6*h*（K1+2*K2+2*K3+K4）;
f（k+1）=U（k+1）-2*T（k+1）/U（k+1）;
end
fork=5:
11
U（k）=U（k-1）+h/24*（55*f（k-1）-59*f（k-2）+37*f（k-3）-9*f（k-4））;
f（k）=U（k）-2*T（k）/U（k）;
U（k）=U（k-1）+h/24*（9*f（k）+19*f（k-1）-5*f（k-2）+f（k-3））;
f（k）=U（k）-2*T（k）/U（k）;
end
Un=U;
fork=1:
11
zz（k）=sqrt（1+2*T（k））;
end
e=U-zz;
end
运行结果及图形显示：
运行得：
>>[Un,e]=di13ti（）
Un=
Columns1through8
1.00001.09541.18321.26491.34161.41421.48321.5492
Columns9through11
1.61251.67331.7321
e=
1.0e-05*
Columns1through8
00.04170.07890.11640.05710.02710.01270.0042
Columns9through11
-0.0013-0.0054-0.0088
编写a5.m用于画图显示计算结果，
t0=0:
0.001:
1;
t1=0:
0.1:
1;
U0=sqrt（1+2*t0）;
U1=di13ti（）;
holdon
plot（t0,U0,'k'）
plot（t1,U1,'--r','Marker','x'）
legend（'U=sqrt（1+2*t）','U1','Location','NorthWest'）
gridon
xlable（'t'）
ylable（'Un'）
holdoff
运行结果如下：
图8
放大如下图9
由图可知Adams四阶预估-校正算法的PECE模式求其数值解精度相当好
结果：
表格一Adams四阶预估-校正算法的PECE模式数值解及其局部误差项
数值解
T0
T1
T2
T3
T4
T5
T6
T7
T8
T9
T10
数值解
1.0000
1.1094
1.1832
1.2649
1.3416
1.4142
1.4832
1.5492
1.6125
1.6733
1.7321
误差项
0
0.000000417
0.000000789
0.000001161
0.000000571
0.000000271
0.000000127
0.000000042
-0.000000013
-0.000000054
-0.000000088
综上图形结果及表格显示，Adams四阶预估-校正算法的PECE模式数值解与精确解相当接近。
升级会员 