原创力文档- 1、原创力文档（book118）网站文档一经付费（服务费），不意味着购买了该文档的版权，仅供个人/单位学习、研究之用，不得用于商业用途，未经授权，严禁复制、发行、汇编、翻译或者网络传播等，侵权必究。。
- 2、本站所有内容均由合作方或网友上传，本站不对文档的完整性、权威性及其观点立场正确性做任何保证或承诺！文档内容仅供研究参考，付费前请自行鉴别。如您付费，意味着您自己接受本站规则且自行承担风险，本站不退款、不进行额外附加服务；查看《如何避免下载的几个坑》。如果您已付费下载过本站文档，您可以点击 这里二次下载。
- 3、如文档侵犯商业秘密、侵犯著作权、侵犯人身权等，请点击“版权申诉”（推荐），也可以打举报电话：400-050-0827(电话支持时间：9:00-18:30)。
- 4、该文档为VIP文档，如果想要下载，成为VIP会员后，下载免费。
- 5、成为VIP后，下载本文档将扣除1次下载权益。下载后，不支持退款、换文档。如有疑问请联系我们。
- 6、成为VIP后，您将拥有八大权益，权益包括：VIP文档下载权益、阅读免打扰、文档格式转换、高级专利检索、专属身份标志、高级客服、多端互通、版权登记。
- 7、VIP文档为合作方或网友上传，每下载1次， 网站将根据用户上传文档的质量评分、类型等，对文档贡献者给予高额补贴、流量扶持。如果你也想贡献VIP文档。上传文档
查看更多
用Euler法改进的Euler法及四阶的龙格库塔法求解初值问题
微分方程数值解课程设计题目
1（30分）分别用Euler法、改进的Euler法及四阶的龙格库塔法求解初值问题：
h分别取0.1和0.2，要求计算并绘出图形，最后比较三种算法的精度。
（1）.先构建初值问题的函数
M文件：
function z=fun(x,y)
z=y-2*x/y;
（2）.Euler法
M文件：
function E=Euler(fun,x0,y0,h,N)
x=zeros(1,N+1);y=zeros(1,N+1);
x(1)=x0;y(1)=y0;
for n=1:N
x(n+1)=x(n)+h;
y(n+1)=y(n)+h*feval(fun,x(n),y(n));
end
T=[x;y]
.当h=0.1时
Euler(fun,0,1,0.1,10)
T = Columns 1 through 10
0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000
1.0000 1.1000 1.1918 1.2774 1.3582 1.4351 1.5090 1.5803 1.6498 1.7178
Column 11
1.0000
1.7848
②.当h=0.2时
Euler(fun,0,1,0.2,5)
T = 0 0.2000 0.4000 0.6000 0.8000 1.0000
1.0000 1.2000 1.3733 1.5315 1.6811 1.8269
（3）.改进的Euler法：
M文件：Euler_modify.m
function E=Euler_modify(fun,x0,y0,h,N)
x=zeros(1,N+1);y=zeros(1,N+1);
x(1)=x0;y(1)=y0;
for n=1:N
x(n+1)=x(n)+h;
z0=y(n)+h*feval(fun,x(n),y(n));
y(n+1)=y(n)+h/2*(feval(fun,x(n),y(n))+feval(fun,x(n+1),z0));
end
T=[x;y]
.当h=0.1时
Euler_modify(fun,0,1,0.1,10)
T = Columns 1 through 9
0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000
1.0000 1.0959 1.1841 1.2662 1.3434 1.4164 1.4860 1.5525 1.6165
Columns 10 through 11
0.9000 1.0000
1.6782 1.7379
②.当h=0.2时
Euler_modify(fun,0,1,0.2,5)
T = 0 0.2000 0.4000 0.6000 0.8000 1.0000
1.0000 1.1867 1.3483 1.4937 1.6279 1.7542
（4）.四阶R_K方法：
M文件：
function [x,y]=Rk_N4(f,x0,y0,h,N)
x=zeros(1,N+1);y=zeros(1,N+1);
x(1)=x0;y(1)=y0;
for n=1:N
x(n+1)=x(n)+h;
k1=h*feval(f,x(n),y(n));
k2=h*feval(f,x(n)+1/2*h,y(n)+1/2*k1);
k3=h*feval(f,x(n)+1/2*h,y(n)+1/2*k2);
k4=h*feval(f,x(n)+h,y(n)+k3);
y(n+1)=y(n)+1/6*(k1+2*k2+2*k3+k4);
end
.当h=0.1时
[x,y]=Rk_N4(fun,0,1,0.1,10)
x = Columns 1 through 9
0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000
Columns 10 through 1
文档评论（0）