《常微分方程》课内实验任务书
学生姓名:张学阳1009300132
及学号：
学院:理学院
班级:数学101
课程名称：常微分方程
题目:Euler方法及其改进方法
指导教师
李鹏松教授
姓名及职称：
朱秀丽讲师
方向实验师
2012年10月17日
《常微分方程》课内实验
实验一 Euler 方法及其改进方法
一、实验目的
1．通过用Matlab 编程运用Euler 方法及其改进方法求解常微分方程初值问题，更进一步掌握常微分方程及其数值解法课程的理论内容，加深对数值解法的理解。
2．熟悉Matlab 编程环境。
二、实验学时和类型
本次实验为2学时、设计性实验。
三、实验内容
1.实验题目
运用Euler 方法及其改进方法求解常微分方程初值问题，具体问题在课内实验习题中学生任意抽取。
2.实验原理
Euler 方法：???+==+)
,()(100n n n n y x hf y y x y y
Euler 改进方法：[]??
?
??++==+++),(),(2)(11100n n n n n n y x f y x f h
y y x y y 3.设计思想
在能够获得精确解的题目抽取题目，分别应用Euler 方法、Euler 改进方法求解数值结果，将所得结果列表或者画图，参照精确结果，对Euler 方法、Euler 改进方法的计算误差进行分析。
4.参考程序源代码
%fun 为目标函数字符串 %x0为自变量初始值。 %y0为fun(x0);
%bou=[a,b]自变量区间 %h 为步长
fun='1/x^2-y^2';
bou=[1,6];
a=bou(1);
b=bou(2);
x0=1;
y0=0;
h=0.1;
n=ceil((b-a)/h);
xx=linspace(a,b,n+1)'; yy=zeros(1,n+1)';
lengthx=length(xx);
xx(1)=x0;yy(1)=y0;
for i=2:n+1
x=xx(i-1);y=yy(i-1);
k=eval(fun);
yy(i)=yy(i-1)+h*k;
end
Ys=dsolve('Dy=1/x^2-y^2','y(1)=0','x');
for i=1:lengthx
x=xx(i);
exacty(i)=eval(Ys);
end
YY=exacty';
yend=[xx,yy,YY]
p=plot(xx,yend(:,2),'k-o','LineWidth',1,...
'MarkerEdgeColor','k',...
'MarkerFaceColor','g',...
'MarkerSize',4);
hold on;
pp=plot(xx,yend(:,3),'r-.+','LineWidth',0.8,...
'MarkerEdgeColor','r',...
'MarkerFaceColor','m',...
'MarkerSize',6);
legend([p,pp],'Eula','jiequejie');
%fun为目标函数字符串
%x0为自变量初始值。
%y0为fun(x0);
%bou=[a,b]自变量区间
%h为步长
fun='1/x^2-y^2';
bou=[1,6];
a=bou(1);
b=bou(2);
x0=1;
y0=0;
h=0.1;
n=ceil((b-a)/h);
xx=linspace(a,b,n+1)'; yy=zeros(1,n+1)';
lengthx=length(xx);
xx(1)=x0;yy(1)=y0;
for i =2:n+1
x=xx(i-1);y=yy(i-1);
k1=eval(fun);
x=xx(i);y=yy(i-1)+h*k1;
k2=eval(fun);
yy(i)=yy(i-1)+h/2*(k1+k2);
end
Ys=dsolve('Dy=1/x^2-y^2','y(1)=0','x');
for i=1:lengthx
x=xx(i);
exacty(i)=eval(Ys);
end
YY=exacty';
yend=[xx,yy,YY]
p=plot(xx,yend(:,2),'k-o','LineWidth',1,...
'MarkerEdgeColor','k',...
'MarkerFaceColor','g',...
'MarkerSize',4);
hold on;
pp=plot(xx,yend(:,3),'r-.+','LineWidth',0.8,...
'MarkerEdgeColor','r',...
'MarkerFaceColor','m',...
'MarkerSize',6);
legend([p,pp],'gjeula','jingquejie');
n=ceil((b-a)/h);
xx=linspace(a,b,n+1)'; yy=zeros(1,n+1)';
xx(1)=x0;yy(1)=y0;
四、注意事项
课后请认真完成实验报告,可以和同学一起研究学习，但坚决杜绝互相抄袭实验报告，发现作弊者实验成绩按零分处理并上报理学院学生工作办公室。
每位同学最后需要上交实验报告和实验指导任务书纸制文件各1份，并提交电子版存档。
CENTRAL SOUTH UNIVERSITY 数值分析实验报告
Euler 方法与改进的Euler 方法的应用 一、问题背景 在工程和科学技术的实际问题中，常需求解微分方程，但常微分方程中往往 只有少数较简单和典型的常微分方程（例如线性常系数常微分方程等）可求出其 解析解,对于变系数常微分方程的解析求解就比较困难，而一般的非线性常微分方 程的求解困难就更不用说了。大多数情况下，常微分方程只能用近似方法求解。 这种近似解法可分为两大类：一类是近似解析法，如级数解法、逐次逼近法等;另 一类是数值解法，它给出方程在一些离散点上的近似值。 二、数学模型 在具体求解微分方程时，需具备某种定解条件，微分方程和定解条件合在一 起组成定解问题。定解条件有两种：一种是给出积分曲线在初始点的状态，称为 初始条件，相应的定解问题称为初值问题。另一类是给出积分曲线首尾两端的状 态，称为边界条件，相应的定解问题称为边值问题。在本文中主要讨论的是给定 初值条件的简单Euler 方法和改进的Euler 方法来求解常微分方程。 三、算法及流程 Euler 方法是最简单的一种显式单步法。对于方程 ()y x f dx dy ,= 考虑用差商代替导数进行计算，取离散化点列 nh x x n +=0，L n ,2,1,0= 则得到方程的近似式 ()()()()n n n n x y x f h x y x y ,1≈-+ 即 ()n n n n y x hf y y ,1+=+ 得到简单Euler 方法。具体计算时由0x 出发，根据初值，逐步递推二得到系列离 散数值。 简单Euler 方法计算量小，然而精度却不高，因而我们可以构造梯形公式 ()()[]η=++ =+++0111,,2 y y t f y t f h y y n n n n n n 其中()N a b h -=。这是一个二阶方法，比Euler 方法精度高。但是上述公式右边 有1+n y ，因而是隐式差分方程，可以用迭代方法计算1+n y 。初值可以由Euler 公式
姓名：樊元君学号：02 日期： 一、实验目的 掌握MATLAB语言、C/C++语言编写计算程序的方法、掌握改进欧拉法与四阶龙格-库塔求解一阶常微分方程的初值问题。掌握使用MATLAB程序求解常微分方程问题的方法。 ： 二、实验内容 1、分别写出改进欧拉法与四阶龙格-库塔求解的算法，编写程序上机调试出结果，要求所编程序适用于任何一阶常微分方程的数值解问题，即能解决这一类问题，而不是某一个问题。 实验中以下列数据验证程序的正确性。 求，步长h=。 * 2、实验注意事项 的精确解为，通过调整步长，观察结果的精度的变化 ^ )
三、程序流程图： ●改进欧拉格式流程图： ~ |
●四阶龙格库塔流程图： ] 四、源程序： ●改进后欧拉格式程序源代码： function [] = GJOL(h,x0,y0,X,Y) format long h=input('h='); … x0=input('x0='); y0=input('y0='); disp('输入的范围是：'); X=input('X=');Y=input('Y='); n=round((Y-X)/h); \
i=1;x1=0;yp=0;yc=0; for i=1:1:n x1=x0+h; yp=y0+h*(-x0*(y0)^2);%yp=y0+h*(y0-2*x0/y0);% · yc=y0+h*(-x1*(yp)^2);%yc=y0+h*(yp-2*x1/yp);% y1=(yp+yc)/2; x0=x1;y0=y1; y=2/(1+x0^2);%y=sqrt(1+2*x0);% fprintf('结果=%.3f,%.8f,%.8f\n',x1,y1,y); : end end ●四阶龙格库塔程序源代码： function [] = LGKT(h,x0,y0,X,Y) 。 format long h=input('h='); x0=input('x0='); y0=input('y0='); disp('输入的范围是：'); " X=input('X=');Y=input('Y='); n=round((Y-X)/h); i=1;x1=0;k1=0;k2=0;k3=0;k4=0; for i=1:1:n ~ x1=x0+h; k1=-x0*y0^2;%k1=y0-2*x0/y0;% k2=(-(x0+h/2)*(y0+h/2*k1)^2);%k2=(y0+h/2*k1)-2*(x0+h/2)/(y0+h/2*k1);% k3=(-(x0+h/2)*(y0+h/2*k2)^2);%k3=(y0+h/2*k2)-2*(x0+h/2)/(y0+h/2*k2);% k4=(-(x1)*(y0+h*k3)^2);%k4=(y0+h*k3)-2*(x1)/(y0+h*k3);% … y1=y0+h/6*(k1+2*k2+2*k3+k4);%y1=y0+h/6*(k1+2*k2+2*k3+k4);% x0=x1;y0=y1; y=2/(1+x0^2);%y=sqrt(1+2*x0);% fprintf('结果=%.3f,%.7f,%.7f\n',x1,y1,y); end · end
P77 31.利用改进欧拉方法计算下列初值问题，并画出近似解的草图：dy + =t = t y y ≤ ≤ ,2 ;5.0 0,3 )0( )1(= ,1 ? dt 代码： %改进欧拉法 function Euler(t0,y0,inv,h) n=round(inv(2)-inv(1))/h; t(1)=t0; y(1)=y0; for i=1:n y1(i+1)=y(i)+h*fun(t(i),y(i)); t(i+1)=t(i)+h; y(i+1)=y(i)+1/2*h*(fun(t(i),y(i))+ fun(t(i+1),y1(i+1))) end plot(t,y,'*r') function y=fun(t,y); y=y+1; 调用：Euler(0,3,[0,2],0.5) 得到解析解：hold on; y=dsolve('Dy=y+1','(y(0)=3)','t'); ezplot(y,[0,2]) 图像：
dy y =t - t y ;2.0 t = ≤ )0( 0,5.0 ,4 )2(2= ≤ ? ,2 dt 代码： function Euler1(t0,y0,inv,h) n=round(inv(2)-inv(1))/h; t(1)=t0; y(1)=y0; for i=1:n y1(i+1)=y(i)+h*fun(t(i),y(i)); t(i+1)=t(i)+h; y(i+1)=y(i)+1/2*h*(fun(t(i),y(i))+ fun(t(i+1),y1(i+1))) end plot(t,y,'*r') function y=fun(t,y); y=y^2-4*t; 调用： Euler1(0,0.5,[0,2],0.2) 图像：
法 function [x,y]=naeuler(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2); y(1)=y0; for n=1:length(x)-1 y(n+1)=y(n)+h*feval(dyfun,x(n),y(n)); end x=x';y=y'; x1=0::1;y1=(1+2*x1).^; plot(x,y,x1,y1) >> dyfun=inline('y-2*x/y'); [x,y]=naeuler(dyfun,[0,1],1,;[x,y] ans = 2.隐式Euler法 function [x,y]=naeulerb(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2); y(1)=y0; for n=1:length(x)-1 y(n+1)=iter(dyfun,x(n+1),y(n),h); end x=x';y=y'; x1=0::1;y1=(1+2*x1).^; plot(x,y,x1,y1) function y=iter(dyfun,x,y,h) y0=y;e=1e-4;K=1e+4; y=y+h*feval(dyfun,x,y); y1=y+2*e;k=1; while abs(y-y1)>e y1=y; y=y0+h*feval(dyfun,x,y); k=k+1; if k>K error('迭代发散'); end end >> dyfun=inline('y-2*x/y');
[x,y]=naeulerb(dyfun,[0,1],1,;[x,y] ans = 3.改进Euler法 function [x,y]=naeuler2(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2); y(1)=y0; for n=1:length(x)-1 k1=feval(dyfun,x(n),y(n)); y(n+1)=y(n)+h*k1; k2=feval(dyfun,x(n+1),y(n+1)); y(n+1)=y(n)+h*(k1+k2)/2; end x=x';y=y'; x1=0::1;y1=(1+2*x1).^; plot(x,y,x1,y1) >> dyfun=inline('y-2*x/y'); [x,y]=naeuler2(dyfun,[0,1],1,;[x,y] ans =
实验九欧拉方法 信息与计算科学金融崔振威201002034031 一、实验目的： 1、掌握欧拉算法设计及程序实现 二、实验内容： 1、p364-9.2.4、p386-9.5.6 三、实验要求： 主程序： 欧拉方法（前项）： function [x,y]=DEEuler(f,a,b,y0,n); %f:一阶常微分方程的一般表达式的右端函数 %a:自变量的取值下限 %b:自变量的取值上限 %y0：函数的初值 %n:积分的步数 if nargin<5,n=50; end h=(b-a)/n; x(1)=a;y(1)=y0; for i=1:n x(i+1)=x(i)+h; y(i+1)=y(i)+h*feval(f,x(i),y(i)); end format short 欧拉方法（后项）： function [x,y]=BAEuler(f,a,b,y0,n); %f:一阶常微分方程的一般表达式的右端函数 %a:自变量的取值下限 %b:自变量的取值上限 %y0：函数的初值 %n:积分的步数 if nargin<5,n=50; end h=(b-a)/n; x(1)=a;y(1)=y0; for i=1:n x(i+1)=x(i)+h; y1=y(i)+h*feval(f,x(i),y(i));
y(i+1)=y(i)+h*feval(f,x(i+1),y1); end format short 梯形算法： function [I,step,h2] = CombineTraprl(f,a,b,eps) %f 被积函数 %a,b 积分上下限 %eps 精度 %I 积分结果 %step 积分的子区间数 if(nargin ==3) eps=1.0e-4; end n=1; h=(b-a)/2; I1=0; I2=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h; while abs(I2-I1)>eps n=n+1; h=(b-a)/n; I1=I2; I2=0; for i=0:n-1 x=a+h*i; x1=x+h; I2=I2+(h/2)*(subs(sym(f),findsym(sym(f)),x)+subs(sym(f),findsym(sym(f)),x1)); end end I=I2; step=n; h2=(b-a)/n; 改进欧拉方法： function [x,y]=MoEuler(f,a,b,y0,n); %f:一阶常微分方程的一般表达式的右端函数 %a:自变量的取值下限 %b:自变量的取值上限 %y0：函数的初值 %n:积分的步数 if nargin<5,n=50; end
题目再现： 2. 当病人采取服用口服药或肌肉注射来治疗疾病时，药物虽然瞬间进入了体内，但它一般都集中与身体的某一部位，靠其表面与肌体接触而逐步被吸收。假定身体系统是一个单房室系统，设t 时刻体内药物的总量为x(t)，则x(t)满足： 问题分析：运用欧拉公式求微分方程的数值解。微分方程为： 11dx ,(0)0dt k t k De kx x -=-= 1，欧拉折线法： 设h=0.1,即n=201时，1k 0.6=，k 0.2=，D=200.MATLAB 程序如下所示： clear f=sym('0.6*200*exp(-1*0.6*t)-0.2*x '); a=0; b=20; h=0.1; n=(b-a)/h+1; t=0; x=0; szj=[t,x]; for i=1:n-1 11dx ,(0)0 dt k t k De kx x -=-=10011100 ,(,)()k t k k k k k k k t x x x h f t x x h k De kx t t h -++==??=+=+-??=+?
x=x+h*subs(f,{'t','x'},{t,x}); t=t+h; szj=[szj;t,x]; end szj ; x=dsolve(‘Dx=120*exp(-0.6*t)-0.2*x’,'x(0)=0','t') T=[0:0.1:20]; X=subs(x,T); plot(szj(:,1),szj(:,2),'or-',T,X,’b-’) 输出结果为： 解析解：x= 其中红线代表的是数值解，蓝线代表的是解析解。可见，数值解的误差还是比较小的。
生物信息技术0801 徐聪U200812594 #include
北京理工大学珠海学院实验报告 ZHUHAI CAMPAUS OF BEIJING INSTITUTE OF TECHNOLOGY 班级2012电气2班学号120109021010姓名陈冲指导教师张凯成绩 实验题目（实验五）欧拉法实验地点及时间JD501 2014/1/2（6-7节） 一、实验目的 1.掌握用程序语言来编辑函数。 2.学会用MATLAB编写Euler.m以及TranEuler.m函数。 二、实验环境 Matlab软件 三、实验内容 1、以书中第124页题目11为例编辑程序来实现计算结果。 2、使用MATLAB进行编写： 第一步：编写Euler.m函数，代码如下 编写TranEuler.m函数，代码如下 第二步：利用上述函数编辑命令：（可见实验结果中的截图）
在此之前先建立一个名为f.m 的M 文件，代码如下 function z=f(x); z=8-3y; 再编辑代码： 得到了欧拉法的结果：y （0.4）=2.47838030901267 编辑另一段命令： 得到改进欧拉法的结果：y （0.4）=2.46543714659780 在此基础上，我还编辑龙格库达的命令窗口代码，如下： 四、实验题目 用欧拉法和改进欧拉法求解初值问题'83,(0)2y y y =-=，试取步长0.2h =计算(0.4)y 的近似值。 五、实验结果
六、总结 通过这次实验我掌握了将得到的解进一步精确，而且要学会比较这几种方法的精确性，显然，四阶龙格库达比改进欧拉发精确，改进欧拉发比欧拉法精确。 实验难度不大，要比较n的取值不同，产生的影响不同。
西华数学与计算机学院上机实践报告 课程名称：计算方法A 年级：2010级 上机实践成绩： 指导教师：严常龙 姓名：李国强 上机实践名称：解常微分方程初值问题 学号：362011********* 上机实践日期：2013.12.25 上机实践编号：8 上机实践时间：14:00 一、目的 1．通过本实验加深对欧拉法、改进欧拉法、线性多步法的构造过程的理解； 2．能对上述四种方法提出正确的算法描述编程实现，观察计算结果的改善情况。 二、内容与设计思想 自选常微分方程的初值问题，分别用欧拉法、改进欧拉法求解。 分别用以上两种方法求解常微分方程初值问题： 2 '()1,([0.0,1.4],0.1)(0)0.0 y x y h y ?=+=?=?求解区间取步长 三、使用环境 操作系统：Win 8 软件平台：Visual C++ 6.0 四、核心代码及调试过程 #include
int i; float y[n+1]; y[0]=m;//赋初值 printf("改进欧拉法\n"); for(i=0;i 设系统方程为：y t y y /2)1(-=,1)0(=y ,用改进欧拉法求解各离散点y 的数值解，步长 10,1.0≤≤=t h ，解析解为t y 21+= 。 解：改进欧拉法 ),(1n n n p n y t hf y y +=+ )],(),([5.0111p n n n n n c n y t f y t f h y y +++++= 已知 n n n n n y t y y t f /2),(-= n n n n n n n p n y ht y h y t y h y y /2)1()/2(1-+=-+=+ 1 111111/5.0/)5.01()]/2()/2[(5.0+++++++-+-+=-+-+=n n n n n n n n n n n n n c n y ht hy y ht y h y t y y t y h y y 程序： h=0.1; t=0:h:1; N=length(t); y=ones(1,N); ey=ones(1,N); zy=ones(1,N); for k=1:N-1 y(1,k+1)=(1+h)*y(1,k)-(2*h*(k-1)/(N-1))./y(1,k);%预估公式 ey(1,k+1)=(1+h)*ey(1,k)-(2*h*(k-1)/(N-1))./ey(1,k);%欧拉公式 y(1,k+1)=(1+0.5*h)*y(1,k)-(h*(k-1)/(N-1))./y(1,k)+0.5*h*y(1,k+1)-(h*k/(N-1))./y(1,k+1);%改进欧拉 zy(1,k+1)=(1+2*k/(N-1)).^0.5;%解析解 end plot(t,zy,'-xk',t,y,':ob',t,ey,'-.*r','linewidth',1.0); xlabel('t'); ylabel('y'); 截图： 7. 取h =0.2，用改进欧拉方法求解下列初值问题。 '(0)5 y y ì?í?=?020x ＃) 第七题: #include "fstream.h" #include "math.h" int main() { double x0,x1,y0,yp,yc,h; ifstream infile("in.dat"); ofstream outfile("out.dat"); infile>>x0>>x1>>h>>y0;//依次输入x 的初值x0，x 的终值x1，步长h ，y （0）的值y0 outfile< 0 5 0.2 7.03239 0.4 9.89185 0.6 13.9148 0.8 19.5739 1 27.5337 1. 2 38.7287 1.4 54.4735 1.6 76.6169 1.8 107.759 2 151.558 2.2 213.156 2.4 299.787 2.6 421.626 2.8 592.981 3 833.977 3.2 1172.91 3. 4 1649.6 3.6 2320.01 3.8 3262.89 4 4588.97 4.2 6453.97 4.4 9076.93 4.6 12765.9 4.8 17954.1 5 25250.8 5.2 35513 5.4 49945.8 5. 6 70244.3 5.8 98792.3 差分法欧拉格式浅谈 科学计算中常常求解常微分方程的定解这类问题的最简形式是一阶方程的初值问题 ???=='00 )(),(y x y y x f y （1） 这里假定右函数),(y x f 适当光滑，譬如关于y 满足Lipschitz 条件，以保证上述初值问题的解y(x)存在且唯一。 虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特殊类型的方程。求解从实际问题当中归结出来的微分方程主要靠数值解法。 差分方法是一类重要的数值解法。这类方法回避解y(x)的函数表达式，而是寻求它在一系列离散节点 <<<< 《常微分方程》课内实验任务书 学生姓名:张学阳1009300132 及学号： 学院:理学院 班级:数学101 课程名称：常微分方程 题目:Euler方法及其改进方法 指导教师 李鹏松教授 姓名及职称： 朱秀丽讲师 方向实验师 2012年10月17日 《常微分方程》课内实验 实验一 Euler 方法及其改进方法 一、实验目的 1．通过用Matlab 编程运用Euler 方法及其改进方法求解常微分方程初值问题，更进一步掌握常微分方程及其数值解法课程的理论内容，加深对数值解法的理解。 2．熟悉Matlab 编程环境。 二、实验学时和类型 本次实验为2学时、设计性实验。 三、实验内容 1.实验题目 运用Euler 方法及其改进方法求解常微分方程初值问题，具体问题在课内实验习题中学生任意抽取。 2.实验原理 Euler 方法：???+==+) ,()(100n n n n y x hf y y x y y Euler 改进方法：[]?? ? ??++==+++),(),(2)(11100n n n n n n y x f y x f h y y x y y 3.设计思想 在能够获得精确解的题目抽取题目，分别应用Euler 方法、Euler 改进方法求解数值结果，将所得结果列表或者画图，参照精确结果，对Euler 方法、Euler 改进方法的计算误差进行分析。 4.参考程序源代码 %fun 为目标函数字符串 %x0为自变量初始值。 %y0为fun(x0); %bou=[a,b]自变量区间 %h 为步长 fun='1/x^2-y^2'; bou=[1,6]; a=bou(1); b=bou(2); x0=1; y0=0; h=0.1; n=ceil((b-a)/h); xx=linspace(a,b,n+1)'; yy=zeros(1,n+1)'; lengthx=length(xx); xx(1)=x0;yy(1)=y0; for i=2:n+1 x=xx(i-1);y=yy(i-1); k=eval(fun); yy(i)=yy(i-1)+h*k; end Ys=dsolve('Dy=1/x^2-y^2','y(1)=0','x'); for i=1:lengthx x=xx(i); exacty(i)=eval(Ys); end YY=exacty'; yend=[xx,yy,YY] p=plot(xx,yend(:,2),'k-o','LineWidth',1,... 'MarkerEdgeColor','k',... 'MarkerFaceColor','g',... 'MarkerSize',4); hold on; pp=plot(xx,yend(:,3),'r-.+','LineWidth',0.8,... 'MarkerEdgeColor','r',... 'MarkerFaceColor','m',... 'MarkerSize',6); legend([p,pp],'Eula','jiequejie'); %fun为目标函数字符串 %x0为自变量初始值。 %y0为fun(x0); %bou=[a,b]自变量区间 %h为步长 fun='1/x^2-y^2'; bou=[1,6]; a=bou(1); b=bou(2); x0=1; y0=0; 用改进的欧拉方法和四阶龙格-库塔方法解初值问题 一、题目： 取步长2.0=h ，分别用改进的欧拉方法和四阶龙格—库塔方法解初值问题 ? ??=≤≤+=.1)0(;10,'y x y x y 并比较结果。 二、基本思想： 1.改进的欧拉方法 改进的欧拉方法用梯形公式计算()()()()+1+1-=,.n n x n n x y x y x f x y x dx ?中的积分，并以 n y 和+1n y 分别表示()n y x 和()+1n y x 的近似值，就得到 ()()()+1+1+1=+,+,.2 n n n n n n h y y f x y f x y （梯形公式） 梯形公式也是一个一步法公式。由于公式右端也含有未知的+1n y ，故被称作是隐式的。隐式格式实际上是关于+1n y 的一个函数方程。为了避免解方程，可以采用欧拉方法计算初始值，再由梯形公式计算。这样建立起来的计算格式称为改进的欧拉格式： ()()()()+1+1+1=+,,=+,+,.2 p n n n n n n n n n y y hf x y h y y f x y f x y ????? 梯形公式也可以采用迭代法求解。如果仍然采用欧拉方法计算迭代初值，那么计算格式就是 []()[]()[]()() 0+1+1+1+1+1=+,,=+,+,.2n n n n k k n n n n n n y y hf x y h y y f x y f x y ????? 由于已假定(),f x y 满足里普希兹条件，所以有 [][][]()[]() [][]+1-1-1+1+1+1,+1+1+1+1+1-=-,-.22 k k k k k k n n n n n n n n h hL y y f x y f x y y y ≤ 从而，迭代的收敛条件是0<<1.2hL 2.四阶龙格-库塔方法 龙格—库塔方法不是用求导数的办法，而是用计算不同点上()y x f ,的函数值，然后对这些函数值作线性拟合，构造近似公式。组合的原则是使得近似公式与泰勒展开式有尽可能多的项吻合，以达到较高的精度。 改进的欧拉格式 第六章数值积分 6.1数值积分基本概念 6.1.1引言 在区间上求定积分 (6.1.1) 是一个具有广泛应用的古典问题，从理论上讲，计算定积分可用Newton-Leibniz公式 (6.1.2) 其中F(x)是被积函数f(x)的原函数.但实际上有很多被积函数找不到用解析式子表达的原函数，例如等等，表面看它们并不复杂， 但却无法求得F(x).此外，有的积分即使能找到F(x)表达式，但式子非常复杂，计算也很困难.还有的被积函数是列表函数，也无法用(6.1.2)的公式计算.而数值积分则只需计算f(x) 在节点xi(i=0，1,…,n)上的值，计算方便且适合于在计算机上机械地实现. 本章将介绍常用的数值积分公式及其误差估计、求积公式的代数精确度、收敛性和稳定性以及Romberg求积法与外推原理等. 6.1.2插值求积公式 根据定积分定义，对及都有 (极限存在)若不取极限，则积分I(f)可近似表示为 (6.1.3) 这里称为求积节点，与f无关，称为求积系数，(6.1.3)称为机械求积公式. 为了得到形如(6.1.3)的求积公式，可在上用Lagrange插值多项式 ，则得 其中 (6.1.4) 这里求积系数由插值基函数积分得到，它与f(x)无关.如果求积公式(6.1.3)中的系数由(6.1.4)给出，则称(6.1.3)为插值求积公式.此时可由插值余项得到 (6.1.5) 这里ξ∈，(6.1.5)称为插值求积公式余项. 当n=1时，，此时 由(6.1.4)可得 于是 (6.1.6) 称为梯形公式.从几何上看它是梯形A bB(见图6-1)的面积近似曲线y=f(x)下的曲边梯形面积，公式(6.1.6)的余项为 (6.1.7) clear clc X=[1 2 3 4 5]; A=[ 1 2 1 0 0; 4 0 0 1 0; 0 4 0 0 1]; C=[2 3 0 0 0 ]; b=[8;16;12]; t=[3 4 5]; B0=A(:,t); while 1 CB0=C(:,t); XN01=X; for i=1:length(t); for j=1:length(X); if XN01(j)==t(i) XN01(j)=0; end end end j=1; for i=1:length(X); if XN01(i)~=0 XN0(j)=XN01(i); j=j+1; end end for j=1:length(XN0); CN0(j)=C(XN0(j)); end N0=[]; for i=1:length(XN0); N0=[N0,A(:,XN0(i))]; end xiN0=CN0-CB0*B0*N0; j=1; z=[]; for i=1:length(xiN0) if xiN0(i)>0 z(j)=i; j=j+1; end end if length(z)+1==1; break; end n=1; for i=1:length(z) if z(i)>z(n) n=i; end end k=XN0(z(n));%换入变量 B=B0*b; P=B0*A(:,k); j=1; for i=1:length(P) if P(i)>0 x(j)=i; j=j+1; end end y=1; for i=1:length(x) if B(x(y))/P(x(y))>B(x(i))/P(x(i)) y=i; end end y1=x(y); y=t(y1);%换出变量 for i=1:length(t) if t(i)==y m=i; break; end end t(m)=k; P2=B0*A(:,k); q=P2(y1); P2(y1)=-1; P2=-P2./q; E=[1 0 0;0 1 0;0 0 1]; E(:,m)=P2; B0=E*B0; end CB0*B0*b Euler法解常微分方程 Euler法解常微分方程算法： Step 1 分别取积分上限、积分下限、步长 Step 2计算判断是否成立，成立转到Step 3，否则继续进行Step 4 Step 3 计算 Step 4 Euler法解常微分方程算程序： function euler2(fun,y0,A,h) %fun--y' %y0---初值 %A----x取值范围 %a----x左区间端点值 %b----x右区间端点值 %h----给定步长 x=min(A); b=max(A); y=y0; while x Step 3 （1）做显性Euler预测 （2）将带入 Step 4计算判断是否成立，成立返回Step 3，否则继续进行Step 5 Step 5 改进Euler法解常微分方程算程序： function gaijineuler2(fun,y0,A,h) %fun--y' %y0---初值 %A----x取值范围 %a----x左区间端点值 %b----x右区间端点值 %h----给定步长 a=min(A); b=max(A); x=a:h:b; y(1)=y0; for i=1:length(x)-1 w1=feval(fun,x(i),y(i)); y(i+1)=y(i)+h*w1; w2=feval(fun,x(i+1),y(i+1)); y(i+1)=y(i)+h*(w1+w2)/2; end x=x' y=y' 例：用改进Euler法计算下列初值问题（取步长h=0.25） 输入：fun=inline('-x*y^2') gaijineuler2(fun,2,[0 5],0.25) 得到： x = 0.2500 0.5000 0.7500 1.0000 1.2500 1.5000 1.7500 2.0000 2.2500 2.5000 2.7500 1)改进欧拉法求解常微分方程的初值问题 #include #include 1.Euler法 function[x,y]=naeuler(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2); y(1)=y0; for n=1:length(x)-1 y(n+1)=y(n)+h*feval(dyfun,x(n),y(n)); end x=x';y=y'; function[x,y]=naeulerb(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2); y(1)=y0; for n=1:length(x)-1 y(n+1)=iter(dyfun,x(n+1),y(n),h); end x=x';y=y'; x1=0:0.2:1;y1=(1+2*x1).^0.5; plot(x,y,x1,y1) >>dyfun=inline('y-2*x/y'); [x,y]=naeuler(dyfun,[0,1],1,0.2);[x,y] ans= 0 1.0000 0.2000 1.2000 0.4000 1.3733 0.6000 1.5315 0.8000 1.6811 1.0000 1.8269 2.隐式Euler法 function y=iter(dyfun,x,y,h) y0=y;e=1e-4;K=1e+4; y=y+h*feval(dyfun,x,y); y1=y+2*e;k=1; while abs(y-y1)>e y1=y; y=y0+h*feval(dyfun,x,y); k=k+1; if k>K error('迭代发散'); end end >>dyfun=inline('y-2*x/y'); [x,y]=naeulerb(dyfun,[0,1],1,0.2);[x,y] ans= 0 1.0000 0.2000 1.1641 0.4000 1.3014 0.6000 1.4146 0.8000 1.5019 1.0000 1.5561 3.改进Euler法 function[x,y]=naeuler2(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2); y(1)=y0; for n=1:length(x)-1 k1=feval(dyfun,x(n),y(n)); y(n+1)=y(n)+h*k1; k2=feval(dyfun,x(n+1),y(n+1)); y(n+1)=y(n)+h*(k1+k2)/2; end x=x';y=y'; >>dyfun=inline('y-2*x/y');matlab 欧拉算法 附截图
数值计算上机第七题关于改进欧拉方法
差分法欧拉格式浅谈
Euler方法及其改进方法
用改进的欧拉方法和四阶龙格-库塔方法解初值问题
证明隐式Euler方法稳定性
改进单纯形法matlab程序
MATLABEuler法解常微分方程
改进欧拉(c语言程序)
欧拉法matlab程序
相关主题
文本预览