常微分方程初值问题数值解法的比较

常微分方程初值问题数值解法的比较数值计算实践—课程 设计 领导形象设计圆作业设计ao工艺污水处理厂设计附属工程施工组织设计清扫机器人结构设计 报告课题名称常微分方程初值问题数值解法的比较完成时间2013-1-17姓名班级学号成绩一.实验目的及 内容 财务内部控制制度的内容财务内部控制制度的内容人员招聘与配置的内容项目成本控制的内容消防安全演练内容 1实验目的：（1）了解常微分方程初值问题的理论背景以及初值问题稳定性、收敛性的研究；（2）熟练掌握欧拉法、改进欧拉法、龙格-库塔法以及截断误差分析；（3）比较欧拉法、改进欧拉法及龙格-库塔法，能够选择合适的方法进行问题的研究计算；y'y2实验内容：求微分方程(0x1)（欧拉法求解）y(0)1'2xyy求微分方程y(0x1)（改进欧拉法求解）y(0)1y'1log(x1)求微分方程(0x1)（龙格-库塔求解）y(0)1根据实验的结果进行分析，了解一般方法的的优缺点，稳定性，收敛性以及截断误差的分析，针对相应问题拿出有效方法得出最优的结果。二．相关背景知识介绍以及初值问题稳定性的研究：在科学与 工程 路基工程安全技术交底工程项目施工成本控制工程量增项单年度零星工程技术标正投影法基本原理 问题中，常微分方程表述物理量的变化规律，应用非常广泛，比如，天体运动的轨迹，机器人控制，化学反应过程的描述和控制以及电路瞬态过程分析等。这些问题中要求解随时间变化的物理量，即位置函数y(t),t表示时间，而微分方程描述了未知函数与它的一阶或高阶导数之间的关系。'考虑一阶常微分方程的初值问题yf(x,y),x[x0,b],y(x0)y0，如果存在实数L0,使1得|f(x,y1)|f(x,y2)L|y1y2|,y1,y2R,则称f关于y满足利普希茨条件，L称为利普希茨常数。'yf(t,y),tt0对于常微分方程初值问题，考虑初值y0的扰动是问题的解y(t)发生偏差的情y(t0)y0形。若t时y(t)的偏差被控制在有界范围内，则称该初值问题是稳定的，否则该初值问题不稳定的。特别地，若t时y(t)的偏差收敛于零，则称该初值问题是渐进稳定的。y'y,tt0(tt0)对于初值问题稳定性的研究，易知其准确解为y(t)y0e，假定初值经过y(t0)y0^(tt0)扰动后变为y0y0，对于扰动后的解为y(t)(y0y0)e因此带来的扰动误差为^(tt0)y(t)y0e，因此考虑t时y(t)的值，它取决于。易知，若0，则原问题是稳定的；若0，则原问题是不稳定的；若0，则原问题是渐进稳定的。实际遇到的大多数常微分方程初值问题都是稳定的，因此在后面的讨论数值解法时这常常是默认条件。1．欧拉法：依据：积分曲线上一点(x,y)的切线斜率等于函数值。方法：推进法，初始点p0(x0,y0)出发，依照方向场在改点的方向推进到x1,x2,...,xn向前欧拉法的得到:(1)将y(x)在xn处泰勒展开''2'''3'y(xn)hy(xn)hy(x1)y(xh)y(x)y(x)h...取h的线性部分，得nnnn2!3!yn1ynhf(xn,yn)（2）将初值问题中得导数用向前差商来代替有'y(xn1)y(xn)y(xn1)y(xn)yf(xn,yn)，因此yn1ynhf(xn,yn)xn1xnhdy（3）将f(x,y)两边同时对x的区间[x,x]上积分dxnn1xn1xn1xn1xn1'ydxf(xn,yn)dx,即y(xn1)y(xn)f(xn,y(x))dx对右端f(xn,y(x))dx用左xnxnxnxn矩形公式得yn1ynhf(xn,yn)，此方法称向前欧拉法，也叫显示欧拉法。2xn1（4）对右端f(xn,y(x))dx用右矩形公式得yn1ynhf(xn1,yn1)，也xn叫隐式欧拉法。误差分析：1.称y(xn1)yn1为计算yn1时的局部截断误差；p12.如果数值方法的局部截断误差为O(h)，那么称这种数值方法的阶数是p，其实p为非负整数。通常情况下，步长h越小，p越高，则局部截断误差越小，计算精''2'''3'y(xn)hy(xn)hy(xn1)y(xnh)y(xn)y(xn)h...y(x)在xn初泰勒展开有2!3!2ynhf(xn,yn)O(h)2则有y(xn1)yn1O(h)可见欧拉方法是一阶方法，精度不是很高。2．改进欧拉方法：xn1梯形公式：对右端f(xn,y(x))dx用梯形公式得+显然梯形公式是隐式公式。xn改进欧拉公式：先用欧拉公式求的一个初步的近似值yn1，成为预测值，预测值yn1的精度可能达不到要求，在用梯形公式将他校正一次，记为yn1，这个结果成为校正值。预测：yn1ynhf(xn,yn)h校正：yn1yn[f(xn,yn)f(xn1,yn1)]2误差分析：记Tn1为改进欧拉公式在xn1处的截断误差，Tn1y(xn1)y(xn)h记Tn1y(xn1)y(xn)[f(xn,y(xn))f(xn1,y(xn1)]因此有2hTn1y(xn1)y(xn)Tn1[f(xn1,y(xn1))f(xn1,y(xn1))]，T1表示在x1出的局2nn3部截断误差。由此得，梯形公式的局部截断误差为Tn1O(h)，因此改进欧拉的截断误差为3Tn1y(xn1)y(xn)O(h)，可见改进欧拉的方法是二阶方法，改进欧拉方法优于欧拉方法。3．龙格—库塔法:'根据拉格朗日微分中值定理，y(xn1)y(xn)y()(xn1xn)y(xn)hf(,y())，***记khf(,y()),得到y(xn1)y(xn)k，这样，只要给出一种计算k的算法，就可以得3yn1ynK1到相应的计算公式。欧拉公式可以写为K1f(xn,yn)hyn1yn(K1K2)2改进欧拉公式可以写成K1f(xn,yn)因此推出一般的推出广式K2f(xn1,ynhK1)yn1ynh(C1K1C2K2...CPKP)K1f(xn,yn)K2f(xna2h,ynhb21K1)称为p阶龙格-库塔方法，简称p阶R-K方法。...p1KPf(xnaph,ynhbpiKi)i1r(xn,yn,h)ciKi,i1因为K1f(xn,yn),这里的ci,i,ij均为常数。i1Kif(xnih,ynhijKj)j1因为给定的系数不唯一，因此这里的常数有无穷多个解，下面是特殊情况下和一般情况下得结果（一阶龙格-库塔）当r=1时，这就是欧拉法。1（二阶龙格-库塔）当r=2时，ci,iij1，这就是改进欧拉法。2三阶和四阶龙格-库塔也只是在一般情况下得结果。hyn1yn(K1K2K3)6K1f(xn,yn)三阶hhK2f(x,yK1)n2n2K3f(xnh,ynhK12hK3)hyn1yn(K12K22K3K4)6K1f(xn,yn)hh四阶K2f(xn,ynK1)22hhK2f(xn,ynK2)22K3f(xnh,ynhK3)4Tn1y(xn1)y(xn)h(1k12k2)其局部截断误差是：三．程序代码欧拉法代码如下：（1）向前欧拉法：functiony=Euler1(fun,x0,y0,xN,N)%fun为一阶微分方程，x0，y0为初始条件，xN为取值范围的一个端点，N为区间个数x=zeros(1,N+1);x=zeros(1,N+1);x(1)=x0;y(1)=y0;h=(xN-x0)/N;%h为区间步长for(n=1:N)x(n+1)=x(n)+h;y(n+1)=y(n)+h*feval(fun,x(n),y(n));%根据向前欧拉公式计算y值endT=[x',y']（2）向后欧拉法：functiony=Euler2(fun,x0,y0,xN,N)%fun为一阶微分方程，x0，y0为初始条件，xN为取值范围的一个端点，N为区间个数x=zeros(1,N+1);x=zeros(1,N+1);x(1)=x0;y(1)=y0;h=(xN-x0)/N;%h为区间步长for(n=1:N)x(n+1)=x(n)+h;z0=y(n)+h*feval(fun,x(n),y(n));for(k=1:3)z1=y(n)+h*feval(fun,x(n+1),z0);if(abs(z1-z0)<1e-3)break;endz0=z1;endy(n+1)=z1;%%根据向后欧拉公式计算y值endT=[x',y']改进欧拉代码如下：function[x,y]=Gaijineuler(f,x0,y0,xZ,h)%f为一阶微分方程的函数；x0，y0为初始条件；xZ为取值范围的一个端点，h为区间步长n=fix((xZ-x0)/h);%计算分点数y(1)=y0;5x(1)=x0;fori=1:nx(i+1)=x0+i*h;yp=y(i)+h*feval(f,x(i),y(i));yc=y(i)+h*feval(f,x(i+1),yp);y(i+1)=(yp+yc)/2;%根据改进欧拉公式计算结果endx=x';y=y';1．3.龙格-库塔代码如下：（三阶龙格-库塔）functionR=Longgekuta3(f,a,b,aZ,h)%a，b为端点，h为步长，aZ为初值n=(b-a)/h;T=zeros(1,n+1);%定义向量Y=zeros(1,n+1);T=a:h:b;%计算各个分点Y(1)=aZ;%初值赋予forj=1:nk1=feval(f,T(j),Y(j));k2=feval(f,T(j)+h/2,Y(j)+k1*h/2);k3=feval(f,T(j)+h,Y(j)-h*k1+k2*2*h);Y(j+1)=Y(j)+(k1+4*k2+k3)*h/6;%根据公式计算endR=[T'Y'];（四阶龙格-库塔）functionR=Longgekuta4(f,a,b,aZ,h)%a，b为端点，h为步长，aZ为初值n=(b-a)/h;T=zeros(1,n+1);%定义向量Y=zeros(1,n+1);T=a:h:b;%计算各个分点Y(1)=aZ;%初值赋予forj=1:nk1=feval(f,T(j),Y(j));k2=feval(f,T(j)+h/2,Y(j)+k1*h/2);k3=feval(f,T(j)+h/2,Y(j)+k2*h/2);k4=feval(f,T(j)+h,Y(j)+k3*h);Y(j+1)=Y(j)+(k1+2*k2+2*k3+k4)*h/6;%根据公式计算endR=[T'Y'];四．数值结果：6输入：定义M文件并保存ffx.mEuler1('ffx',0,1,1,10)结果：T=01.00000.10001.10000.20001.21000.30001.33100.40001.46410.50001.61050.60001.77160.70001.94870.80002.14360.90002.35791.00002.5937ans=Columns1through71.00001.10001.21001.33101.46411.61051.7716Columns8through111.94872.14362.35792.5937向前欧拉公式结果：输入：Euler2('ffx',0,1,1,10)结果T=01.00000.10001.11100.20001.23440.30001.37160.40001.52400.50001.69330.60001.88140.70002.09040.80002.32270.90002.58071.00002.8674ans=Columns1through771.00001.11101.23441.37161.52401.69331.8814Columns8through112.09042.32272.58072.8674改进欧拉公式结果：输入：[x,y]=Gaijineuler('f',0,1,1,0.1)结果：x=00.10000.20000.30000.40000.50000.60000.70000.80000.90001.0000y=1.00001.09591.18411.26621.34341.41641.48601.55251.61651.67821.7379龙格-库塔计算结果：（三阶）输入：Longgekuta3('ff',0,1,1,0.1)结果：ans=Longgekuta3('ff',0,1,1,0.1)8ans=01.00000.10001.10480.20001.21880.30001.34110.40001.47110.50001.60820.60001.75200.70001.90210.80002.05800.90002.21951.00002.3863（四阶）输入：Longgekuta4('ff',0,1,1,0.1)结果：ans=01.00000.10001.10480.20001.21880.30001.34110.40001.47110.50001.60820.60001.752090.70001.90210.80002.05800.90002.21951.00002.3863五．计算结果分析:方法显示欧拉简单精度低隐式欧拉稳定性最好精度低，计算量大梯形公式精度提高计算量大中点公式精度提高，显式多一个初值，可能影响精度1.收敛性：若某算法对于任意固定的x=xi=x0+ih，当h0(同时i)时有yiy(xi)，则称该算法是收敛的。以下讨论的都是单步法（指在计算yn1时只用到它前一步的信息yn）欧拉法，龙格-库塔法都是单步法的例子；'yy,tt0例：就初值问题考察欧拉显式 格式 pdf格式笔记格式下载页码格式下载公文格式下载简报格式下载 的收敛性。y(t0)y0(tt0)解：该问题的精确解为y(t)y0e欧拉公式为yi1yihyi(1h)yixi1hhxi对于固定的x=xi=x0+ih，有yiy0(1h)y0[(1h)]xiy0ey(xi)p设单步法具有p阶精度，满足利普希茨条件，其整体的截断误差：y(xn)ynO(h)判断单步法的收敛性，归结为验证增量函数能否满足利普希茨条件。对于欧拉法，由于其增量函数就是f（x，y），因此当f（x，y）关于y满足利普希茨条件时候它就是收敛的。对于改进欧拉法和龙格-库塔方法，找到增量函数，满足利普希茨条件的时候，找到利普希茨常数，这些方法都是收敛的。2.稳定性：y'30y,例：考察初值问题在区间[0,0.5]上的解，分别用欧拉显，隐式格式，改进欧拉格式计y(0)1算数值解：10节点欧拉显式欧拉隐式改进欧拉法精确解ye30x01.00001.00001.00001.00000.1-2.00002.5*10^(-1)2.50004.9787*10^(-2)0.24.00006.2500*10^(-2)6.25002.4788*10^(-3)0.3-8.00001.5625*10^(-2)1.5625*10^(1)1.2341*10^(-4)0.41.6000*10^(1)3.9063*10^(-3)3.9063*10^(1)6.1442*10^(-6)0.5-3.2000*10^(1)9.7656*10^(-1)9.7656*10^(1)3.0590*10^(-7)若某算法在计算过程中任一步产生的误差在以后的计算中都逐步衰减，则称该算法是绝对稳定的.考察显示欧拉法：i1i1yi1yihyi(1h)y0,0y0y0,yi1yi1(1h)0，由此可见，要保证初始误差0以后逐步衰减hh必须满足：|1h|1i111考察后退欧拉法：yi1yihyi1,yi1yi,i10由此可见，1h1h据对绝对稳定绝对稳定绝对稳定的区域为|1h|1'h梯形公式：模型yy,yn1yn(ynyn1)，绝对稳定区域为2|2h||2h|相对稳定区间为h0yn1ynhf(xn,yn)(1h)yn,改进欧拉公式：h12yn1yn[f(xn,yn)f(xn1,yn1)[1h(h)]yn2212相对稳定区域为|1h(h)|1稳定区间为2h02121314四阶龙格-库塔：yn1[1h(h)(h)(h)]yn2624121314相对稳定区域：[1h(h)(h)(h)]12624相对稳定区间：2.78h0y'12x,考察初值问题在区间[0,1]上的解y(0)1节点改进欧拉法四节龙格-库塔法准确解1101110.21.1866671.1832291.1832160.41.3483121.3416671.3416410.61.4937041.4832811.4832400.81.6278611.6125141.61245211.7542051.7321421.732051变步长的龙格-库塔法：实际计算中，步长h的选择是一个十分重要的问题，若步长h取的太小，将增加许多不必要的的计算量，若步长取的太大，其计算竞速很难保证，因此，与数值积分一样，微分方程的(h)5数值解也有一个选择步长的问题。由四阶龙格-库塔法的性质，其局部截断误差为y(xn1)yn1ch这5h里的c与y(x)在[x,x1]内的值有关系。然后将步长折半，即取为步长，从x跨两步到x1，在求nn2nnh5h()()2h2h5一个近似值y，每跨一步的误差为c，因此y(x1)y12c()，比较之后n12nn2h()2y(xn1)yn111(h)误差的大约减少。y(xn1)yn11616六．计算时出现的问题，解决方法及体会：对于初值问题，在用不同方法计算离散点的函数值的时候，有精确值给定的情况下，不同方法跟精确值的差别很大，也就是不是一种好方法，这与实际问题的稳定性有很大的关系，只要所求的初值问题，能够落在使用方法的稳定区间内，所得到的结果和精确值的离散点的函数值相差不会很大，并且在稳定区间内，改进欧拉法优于欧拉法，三阶之后的龙格-库塔方法优于改进欧拉法，离散点的函数值更加接近精确值。所谓的欧拉法就是一阶的龙格-库塔方法，改进欧拉法也是在显示欧拉法预测和梯形公式校正的基础上作出的改进，也是二阶龙格-库塔的特殊情况。所谓的三阶，四阶甚至n阶龙格-库塔方法，也是在待定系数有无数解得情况下的一般情况。一阶常微分初值问题的收敛性，只要其增量函数关于y能满足利普希茨条件，这就收敛。而起在计算过程过有避免不了的误差，欧拉法精度相对来说比较低，龙格-库塔的阶数越高，精度就越高，所得的结果和精确值越接近，截断误差越小。同一个问题，想要减少截断误差，我们可以采用变步长的的方法，将步长变为h，这样得出的结果，根据不同的方法，截断误差会有相应的改变。2而在每种方法的程序设计方面，主要思想还是根据每种方法的公式，找到步长，分点数，用循环语句，得出相应的结果，只是在龙格-库塔方法里面的三阶，四阶的公式里面的待定系数有无数个解，我们给出的也是一般情况的三阶，四阶龙格-库塔方法，精度，误差相对于欧拉法，改进欧拉法，梯形公式都优于他们。数值计算实践—课程设计原始数据记录实验名称：常微分方程初值问题数值解法研究实验时间：2013年1月14日12输入：定义M文件并保存ffx.mEuler1('ffx',0,1,1,10)结果：T=01.00000.10001.10000.20001.21000.30001.33100.40001.46410.50001.61050.60001.77160.70001.94870.80002.14360.90002.35791.00002.5937ans=Columns1through71.00001.10001.21001.33101.46411.61051.7716Columns8through111.94872.14362.35792.5937向前欧拉公式结果：输入：Euler2('ffx',0,1,1,10)结果T=01.00000.10001.11100.20001.23440.30001.37160.40001.52400.50001.69330.60001.88140.70002.09040.80002.32270.90002.58071.00002.8674ans=13Columns1through71.00001.11101.23441.37161.52401.69331.8814Columns8through112.09042.32272.58072.8674改进欧拉公式结果：输入：[x,y]=Gaijineuler('f',0,1,1,0.1)结果：x=00.10000.20000.30000.40000.50000.60000.70000.80000.90001.0000y=1.00001.09591.18411.26621.34341.41641.48601.55251.61651.67821.7379龙格-库塔计算结果：（三阶）输入：Longgekuta3('ff',0,1,1,0.1)结果：ans=Longgekuta3('ff',0,1,1,0.1)14ans=01.00000.10001.10480.20001.21880.30001.34110.40001.47110.50001.60820.60001.75200.70001.90210.80002.05800.90002.21951.00002.3863（四阶）输入：Longgekuta4('ff',0,1,1,0.1)结果：ans=01.00000.10001.10480.20001.21880.30001.34110.40001.47110.50001.6082150.60001.75200.70001.90210.80002.05800.90002.21951.00002.3863张昊林做的部分：（1）相关背景知识的介绍（在看书的理解基础上）（1月14）（2）三阶、四阶龙格-库塔思路和程序编写（1月15）（3）截断误差分析（1月16）（4）在理解和认识的基础上对每种方法思路写出的认识和体会（1月17）罗明做的部分：（1）实验目的及内容（2）欧拉法、改进欧拉法思路和程序编写。（1月14）（3）稳定性的研究（1月15）（4）收敛性的研究（1月16）（5）变步长龙格-库塔法的研究（1月17）教师评语指导教师：年月日指导教师：徐红敏2013年1月17日16