原创力文档- 1、本文档共10页，可阅读全部内容。
- 2、原创力文档（book118）网站文档一经付费（服务费），不意味着购买了该文档的版权，仅供个人/单位学习、研究之用，不得用于商业用途，未经授权，严禁复制、发行、汇编、翻译或者网络传播等，侵权必究。
- 3、本站所有内容均由合作方或网友上传，本站不对文档的完整性、权威性及其观点立场正确性做任何保证或承诺！文档内容仅供研究参考，付费前请自行鉴别。如您付费，意味着您自己接受本站规则且自行承担风险，本站不退款、不进行额外附加服务；查看《如何避免下载的几个坑》。如果您已付费下载过本站文档，您可以点击 这里二次下载。
- 4、如文档侵犯商业秘密、侵犯著作权、侵犯人身权等，请点击“版权申诉”（推荐），也可以打举报电话：400-050-0827(电话支持时间：9:00-18:30)。
...
常微分方程数值解实验报告
学院：数学与信息科学
专业：信息与计算科学
：思义
学号：201216524
课程：常微分方程数值解
....
...
实验一：常微分方程的数值解法
1、分别用Euler法、改进的Euler法（预报校正格式）和S—K法求
解初值问题。（h=0.1）并与真解作比较。
yyx1
y(0）1
1.1实验代码：
%欧拉法
function[x,y]=naeuler(dyfun,xspan,y0,h)
%dyfun是常微分方程，xspan是x的取值围，y0是初值，是步长
x=xspan(1):h:xspan(2);
y(1)=y0;
forn=1:length(x)-1
y(n+1)=y(n)+h*feval(dyfun,x(n),y(n));
end
%改进的欧拉法
function[x,m,y]=naeuler2(dyfun,xspan,y0,h)
%dyfun是常微分方程，xspan是x的取值围，y0是初值，是步长。
%返回值x为x取值，m为预报解，y为校正解
x=xspan(1):h:xspan(2);
y(1)=y0;
m=zeros(length(x)-1,1);
forn=1:length(x)-1
k1=feval(dyfun,x(n),y(n));
y(n+1)=y(n)+h*k1;
m(n)=y(n+1);
k2=feval(dyfun,x(n+1),y(n+1));
y(n+1)=y(n)+h*(k1+k2)/2;
end
%四阶S—K法
function[x,y]=rk(dyfun,xspan,y0,h)
%dyfun是常微分方程，xspan是x的取值围，y0是初值，是步长。
x=xspan(1):h:xspan(2);
y(1)=y0;
forn=1:length(x)-1
k1=feval(dyfun,x(n),y(n));
k2=feval(dyfun,x(n)+h/2,y(n)+(h*k1)/2);
k3=feval(dyfun,x(n)+h/2,y(n)+(h*k2)/2);
k4=feval(dyfun,x(n)+h,y(n)+h*k3);
....
...
y(n+1)=y(n)+(h/6)*(k1+2*k2+2*k3+k4);
end
%主程序
x=[0:0.1:1];y=exp(-x)+x;
dyfun=inline(-y+x+1);
[x1,y1]=naeuler(dyfun,[0,1],1,0.1);
[x2,m,y2]=naeuler2(dyfun,[0,1],1,0.1);
[x3,y3]=rk(dyfun,[0,1],1,0.1);
plot(x,y,r,x1,y1,+,x2,y2,*,x3,y3,o);
xlabel(x);ylabel(y);
legend(y为真解,y1为欧拉解,y2为改进欧拉解,y3为S—K解
,Location,NorthWest);
1.2实验结果：
x真解y欧拉解y1预报值m校正值y2S—K解y3
0.01.00001.00001.00001.0000
0.11.00481.00001.00001.00501.0048
0.21.01871.01001.01451.01901.0187
0.31.04081.02901.03711.04121.0408
0.4
文档评论（0）