实验四常微分方程初值问题数值解法Word文档格式.docx
《实验四常微分方程初值问题数值解法Word文档格式.docx》由会员分享，可在线阅读，更多相关《实验四常微分方程初值问题数值解法Word文档格式.docx（17页珍藏版）》请在冰点文库上搜索。
处的近似值
，并与问题的解析解
相比较。
3、用Matlab绘图函数plot（x,y）绘制P285实验任务
（2）的精确解和近似解的图形。
4、使用matlab中的ode45求解P285实验任务
（2），并绘图。
四、实验结果
1、Euler算法程序、改进Euler算法程序；
2、用Euler算法程序、改进Euler算法求解P167例题1的运行结果；
3、Runge－Kutta算法程序；
4、用Runge－Kutta算法求解P178例题、P285实验任务
（2），计算结果如下（其中
表示数值解，
表示解析解，结果保留八位有效数字）：
0.05
0.1
0.15
0.2
0.25
-0.02219009
-0.03877057
-0.04975651
-0.05516258
-0.05500309
-0.03877058
-0.05500310
0.3
0.35
0.4
0.45
0.5
-0.04929202
-0.03804297
-0.02126924
0.00101623
0.02880079
-0.04929202
-0.03804298
-0.02126925
0.00101622
0.02880078
5、P285实验任务
（2）精确解与近似解的图形比较
（图片贴到此处）
6、用matlab中的ode45求解P285实验任务
（2）
MATLABdesktopkeyboardshortcuts,suchasCtrl+S,arenowcustomizable.
Inaddition,manykeyboardshortcutshavechangedforimprovedconsistency
acrossthedesktop.
Tocustomizekeyboardshortcuts,usePreferences.Fromthere,youcanalso
restorepreviousdefaultsettingsbyselecting"
R2009aWindowsDefaultSet"
fromtheactivesettingsdrop-downlist.Formoreinformation,seeHelp.
Clickhereifyoudonotwanttoseethismessageagain.
>
ydot_fun=inline（'
y-2*x./y'
'
x'
y'
）;
[x,y]=euler_f（ydot_fun,0,1,0.1,10）
x=
Columns1through4
00.00000.00000.0000
Columns5through8
0.00000.00000.00000.0000
Columns9through11
0.00000.00001.0000
y=
1.00001.00001.81821.4722
1.02891.77961.63321.5217
1.77111.00871.7982
[x,y]=euler_r（ydot_fun,0,1,0.1,10）
1.00001.90911.29971.5776
1.39991.69091.56691.6146
1.20581.51861.5414
-y+x+1'
[x,y]=euler_f（ydot_fun,0,1,0.1,5）
Columns5through6
0.00000.0000
1.00001.00001.00001.0000
1.00001.0000
[x,y]=euler_r（ydot_fun,0,1,0.1,5）
1.50001.5625
y^2'
[x,y]=Runge_Kutta4（ydot_fun,0,1,0.1,5）
00.00000.00000.0000
1.00001.21941.70151.1445
1.03231.0669
（-y+x^2+4*x-1）./2'
[x,y]=Runge_Kutta4（ydot_fun,0,0,0.05,10）
0.00000.00000.0000
0-0.6823-0.3692-0.8554
-0.6671-0.5772-0.4665-0.0910
-0.30080.75860.9418
ydot_fun=inline（'
[x,y]=Runge_Kutta4（ydot_fun,0,0,0.05,10）;
sprintf（'
%.8f,%.8f\n'
x,y）
ans=
0.00000000,0.05000000
0.10000000,0.15000000
0.20000000,0.25000000
0.30000000,0.35000000
0.40000000,0.45000000
0.50000000,0.00000000
-0.02219009,-0.03877057
-0.04975651,-0.05516258
-0.05500309,-0.04929202
-0.03804297,-0.02126924
0.00101623,0.02880079
xx=[0:
0.05:
0.5];
z=exp（-xx./2）+xx.^2-1;
sprintf（'
xx,z）
-0.02219009,-0.03877058
-0.05500310,-0.04929202
-0.03804298,-0.02126925
0.00101622,0.02880078
xx=[0:
holdon
plot（x,y,'
o'
plot（xx,z,'
*'
holdoff
[x,y]=ode45（ydot_fun,[0,0.5],0）;
[x'
;
]
00.72600.45210.1781
0-0.1419-0.8501-0.1371
0.90420.53430.16450.7947
-0.0158-0.7267-0.1287-0.8058
Columns9through12
0.42490.57590.72680.8778
-0.3417-0.2993-0.0323-0.0400
Columns13through16
0.02880.02880.02880.0288
-0.5763-0.9828-0.1403-0.8426
Columns17through20
-0.1789-0.4532-0.3373-0.1758
Columns21through24
-0.2876-0.8809-0.2155-0.8392
Columns25through28
-0.8617-0.8665-0.0809-0.5467
Columns29through32
-0.3663-0.6119-0.5031-0.5134
Columns33through36
0.02880.02880.02880.0287
-0.5891-0.0564-0.8055-0.3352
Columns37through40
0.02880.02880.02870.0287
-0.9452-0.6404-0.3229-0.7744
Columns41through44
0.02870.02880.02880.0287
-0.8238-0.2486-0.9740-0.9953
Columns45through48
0.02870.02870.02870.0287
-0.5211-0.8719-0.68590.2165
Columns49through52
0.02870.52160.01440.5072
0.71370.33560.59700.5270
Column53
0.0000
0.6419
）
五、讨论分析
（通过最后绘制的图形直观分析近似解与准确解之间的误差，自己补充）
从图中可以看到4阶的Runge-Kutta画出来的点基本和原函数的对应点相重合，精度之高，从数据来看，起码有9位有效数字。
六、改进实验建议
（自己补充）
可以通过输出数据的精度控制函数
实验四常微分方程初值问题数值解法
（这里只提供解答格式请同学自己按照上机实验报告格式填写）
1、饿uler法求解初值问题
function[x,y]=euler_f（ydot_fun,x0,y0,h,N）
%Euler（向前）公式，其中
%ydot_fun---一阶微分方程的函数
%x0,y0---初始条件
%h---区间步长
%N---区间的个数
%x---Xn构成的向量
%y---Yn构成的向量
x=zeros（1,N+1）;
y=zeros（1,N+1）;
x
（1）=x0;
y
（1）=y0;
forn=1:
N
x（n+1）=x（n）+h;
y（n+1）=y（n）+h*feval（ydot_fun,x（n）,y（n））;
end
用欧拉法计算p167例1
x=Columns1through11
00.00000.00000.00000.00000.0000
0.00000.00000.00000.00001.0000
y=.00001.00001.81821.47221.02891.7796
1.63321.52171.77111.00871.7982
2、改进Euler公式
function[x,y]=euler_r（ydot_fun,x0,y0,h,N）
%改进Euler公式，其中
ybar=y（n）+h*feval（ydot_fun,x（n）,y（n））;
y（n+1）=y（n）+h/2*（feval（ydot_fun,x（n）,y（n））+feval（ydot_fun,x（n+1）,ybar））;
用改进欧拉公式计算P167例题1
x=Columns1through6
Columns7through11
y=Columns1through6
1.00001.90911.29971.57761.39991.6909
1.56691.61461.20581.51861.5414
用欧拉法、改进欧拉法计算P171例题3
x=00.00000.00000.00000.00000.0000
y=1.00001.00001.00001.00001.00001.0000
用改进欧拉法计算P171例题3
x=00.00000.00000.00000.00000.0000
y=1.00001.00001.00001.00001.50001.5625
3、标准四阶Runge_Kutta公式
function[x,y]=Runge_Kutta4（ydot_fun,x0,y0,h,N）
%标准四阶Runge_Kutta公式，其中
k1=h*feval（ydot_fun,x（n）,y（n））;
k2=h*feval（ydot_fun,x（n）+1/2*h,y（n）+1/2*k1）;
k3=h*feval（ydot_fun,x（n）+1/2*h,y（n）+1/2*k2）;
k4=h*feval（ydot_fun,x（n）+h,y（n）+k3）;
y（n+1）=y（n）+1/6*（k1+2*k2+2*k3+k4）;
用四阶龙格库塔计算P178例题
y=1.00001.21941.70151.14451.03231.0669
用龙格库塔方法计算P285实验任务
（2）
0.00000.00000.00000.00000.0000
y=Columns1through6
0-0.6823-0.3692-0.8554-0.6671-0.5772
-0.4665-0.0910-0.30080.75860.9418
z=exp（-xx./2）+xx.^2-1;
holdon
holdoff
其图形如图一所示
图一精确解与龙格库塔计算结果比较
4、用Matlab提供的ode45求解P285实验任务
（2）
[x,y]=ode45（ydot_fun,[0,0.5],0）;
plot（x,y,’*’）
ans=Columns1through6
00.72600.45210.17810.90420.53430-0.1419-0.8501-0.1371-0.0158-0.7267
0.16450.79470.42490.57590.72680.8778
-0.1287-0.8058-0.3417-0.2993-0.0323-0.0400
0.02880.02880.02880.02880.02880.0288
-0.5763-0.9828-0.1403-0.8426-0.1789-0.4532
-0.3373-0.1758-0.2876-0.8809-0.2155-0.8392
-0.8617-0.8665-0.0809-0.5467-0.3663-0.6119
0.02880.02880.02880.02880.02880.0287
-0.5031-0.5134-0.5891-0.0564-0.8055-0.3352
0.02880.02880.02870.02870.02870.0288
-0.9452-0.6404-0.3229-0.7744-0.8238-0.2486
0.02880.02870.02870.02870.02870.0287
-0.9740-0.9953-0.5211-0.8719-0.68590.2165
0.02870.52160.01440.50720.0000
0.71370.33560.59700.52700.6419
图二用ode45求解得到的近似解散点图
升级会员 