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

精品文档 常微分方程初值问题数值解法的比较 数值计算实践课程设计报告课题名称常微分方程初值问题数值解法的比较完成时间2013-1-17姓名班级学号成绩一. 实验目的及内容1实验目的：（1） 了解常微分方程初值问题的理论背景以及初值问题稳定性、收敛性的研究；（2） 熟练掌握欧拉法、改进欧拉法、龙格-库塔法以及截断误差分析；（3） 比较欧拉法、改进欧拉法及龙格-库塔法，能够选择合适的方法进行问题的研究计算； 2实验内容：求微分方程（欧拉法求解）求微分方程（改进欧拉法求解）求微分方程（龙格-库塔求解）根据实验的结果进行分析，了解一般方法的的优缺点，稳定性，收敛性以及截断误差的分析，针对相应问题拿出有效方法得出最优的结果。二相关背景知识介绍以及初值问题稳定性的研究： 在科学与工程问题中，常微分方程表述物理量的变化规律，应用非常广泛，比如，天体运动的轨迹，机器人控制，化学反应过程的描述和控制以及电路瞬态过程分析等。这些问题中要求解随时间变化的物理量，即位置函数表示时间，而微分方程描述了未知函数与它的一阶或高阶导数之间的关系。 考虑一阶常微分方程的初值问题 ，如果存在实数使得则称f关于y满足利普希茨条件，L称为利普希茨常数。 对于常微分方程初值问题，考虑初值的扰动是问题的解发生偏差的情形。若时的偏差被控制在有界范围内，则称该初值问题是稳定的，否则该初值问题不稳定的。特别地，若时的偏差收敛于零，则称该初值问题是渐进稳定的。 对于初值问题稳定性的研究，易知其准确解为，假定初值经过扰动后变为，对于扰动后的解为因此带来的扰动误差为，因此考虑时的值，它取决于。易知，若 ，则原问题是稳定的；若 ，则原问题是不稳定的；若 ，则原问题是渐进稳定的。实际遇到的大多数常微分方程初值问题都是稳定的，因此在后面的讨论数值解法时这常常是默认条件。 1 欧拉法：依据：积分曲线上一点的切线斜率等于函数值。方法：推进法，初始点出发，依照方向场在改点的方向推进到 向前欧拉法的得到:(1)将在处泰勒展开取h的线性部分，得（2）将初值问题中得导数用向前差商来代替有 ，因此（3）将两边同时对x的区间上积分 对右端用左矩形公式得，此方法称向前欧拉法，也叫显示欧拉法。（4）对右端用右矩形公式得，也叫隐式欧拉法。误差分析：1.称为计算时的局部截断误差； 2.如果数值方法的局部截断误差为，那么称这种数值方法的阶数是p，其实p为非负整数。通常情况下，步长h越小，p越高，则局部截断误差越小，计算精初泰勒展开有则有可见欧拉方法是一阶方法，精度不是很高。2 改进欧拉方法： 梯形公式：对右端用梯形公式得+显然梯形公式是隐式公式。 改进欧拉公式：先用欧拉公式求的一个初步的近似值，成为预测值，预测值的精度可能达不到要求，在用梯形公式将他校正一次，记为，这个结果成为校正值。 预测： 校正： 误差分析：记为改进欧拉公式在处的截断误差， 记因此有，表示在出的局部截断误差。由此得，梯形公式的局部截断误差为，因此改进欧拉的截断误差为，可见改进欧拉的方法是二阶方法，改进欧拉方法优于欧拉方法。3 龙格库塔法: 根据拉格朗日微分中值定理，记得到，这样，只要给出一种计算的算法，就可以得到相应的计算公式。欧拉公式可以写为 改进欧拉公式可以写成因此推出一般的推出广式称为p阶龙格-库塔方法，简称p阶R-K方法。因为 这里的均为常数。因为给定的系数不唯一，因此这里的常数有无穷多个解，下面是特殊情况下和一般情况下得结果（一阶龙格-库塔）当r=1时，这就是欧拉法。（二阶龙格-库塔）当r=2时，这就是改进欧拉法。三阶和四阶龙格-库塔也只是在一般情况下得结果。三阶四阶其局部截断误差是： 三程序代码 欧拉法代码如下：（1）向前欧拉法：function y=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）向后欧拉法：function y=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; end z0=z1; end y(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;x(1)=x0;for i=1:n x(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;13.龙格-库塔代码如下：（三阶龙格-库塔）function R=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; %初值赋予for j=1:n k1=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;（四阶龙格-库塔）function R=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; %初值赋予for j=1:n k1=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;四数值结果：输入：定义M文件并保存ffx.m Euler1(ffx,0,1,1,10)结果：T = 0 1.0000 0.1000 1.1000 0.2000 1.2100 0.3000 1.3310 0.4000 1.4641 0.5000 1.6105 0.6000 1.7716 0.7000 1.9487 0.8000 2.1436 0.9000 2.3579 1.0000 2.5937ans = Columns 1 through 7 1.0000 1.1000 1.2100 1.3310 1.4641 1.6105 1.7716 Columns 8 through 11 1.9487 2.1436 2.3579 2.5937向前欧拉公式结果：输入：Euler2(ffx,0,1,1,10)结果T = 0 1.0000 0.1000 1.1110 0.2000 1.2344 0.3000 1.3716 0.4000 1.5240 0.5000 1.6933 0.6000 1.8814 0.7000 2.0904 0.8000 2.3227 0.9000 2.5807 1.0000 2.8674ans = Columns 1 through 7 1.0000 1.1110 1.2344 1.3716 1.5240 1.6933 1.8814 Columns 8 through 112.0904 2.3227 2.5807 2.8674改进欧拉公式结果： 输入：x,y=Gaijineuler(f,0,1,1,0.1)结果：x = 0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000y = 1.0000 1.0959 1.1841 1.2662 1.3434 1.4164 1.4860 1.5525 1.6165 1.6782 1.7379龙格-库塔计算结果：（三阶）输入：Longgekuta3(ff,0,1,1,0.1) 结果：ans =Longgekuta3(ff,0,1,1,0.1)ans = 0 1.0000 0.1000 1.1048 0.2000 1.2188 0.3000 1.3411 0.4000 1.4711 0.5000 1.6082 0.6000 1.7520 0.7000 1.9021 0.8000 2.0580 0.9000 2.2195 1.0000 2.3863（四阶）输入：Longgekuta4(ff,0,1,1,0.1) 结果：ans = 0 1.0000 0.1000 1.1048 0.2000 1.2188 0.3000 1.3411 0.4000 1.4711 0.5000 1.6082 0.6000 1.7520 0.7000 1.9021 0.8000 2.0580 0.9000 2.2195 1.0000 2.3863五计算结果分析:方法显示欧拉简单精度低隐式欧拉稳定性最好精度低，计算量大梯形公式精度提高计算量大中点公式精度提高，显式多一个初值，可能影响精度1.收敛性：若某算法对于任意固定的 x = xi = x0 + i h，当 h0 ( 同时 i ) 时有 yi y( xi )，则称该算法是收敛的。以下讨论的都是单步法（指在计算时只用到它前一步的信息）欧拉法，龙格-库塔法都是单步法的例子；例：就初值问题 考察欧拉显式格式的收敛性。解：该问题的精确解为欧拉公式为 对于固定的 x = xi = x0 + i h，有 设单步法具有p阶精度，满足利普希茨条件，其整体的截断误差：判断单步法的收敛性，归结为验证增量函数能否满足利普希茨条件。对于欧拉法，由于其增量函数就是f（x，y），因此当f（x，y）关于y满足利普希茨条件时候它就是收敛的。对于改进欧拉法和龙格-库塔方法，找到增量函数，满足利普希茨条件的时候，找到利普希茨常数，这些方法都是收敛的。 2.稳定性： 例：考察初值问题在区间0,0.5上的解，分别用欧拉显，隐式格式，改进欧拉格式计算数值解：节点欧拉显式欧拉隐式改进欧拉法精确解01.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)若某算法在计算过程中任一步产生的误差在以后的计算中都逐步衰减，则称该算法是绝对稳定的.考察显示欧拉法：，由此可见，要保证初始误差以后逐步衰减必须满足：考察后退欧拉法：由此可见，据对绝对稳定绝对稳定绝对稳定的区域为梯形公式：模型，绝对稳定区域为相对稳定区间为改进欧拉公式：相对稳定区域为稳定区间为四阶龙格-库塔：相对稳定区域：相对稳定区间：考察初值问题在区间0,1上的解节点改进欧拉法四节龙格-库塔法准确解01110.21.1866671.1832291.1832160.41.3483121.3416671.3416410.61.4937041.4832811.4832400.81.6278611.6125141.61245211.7542051.7321421.732051 变步长的龙格-库塔法：实际计算中，步长h的选择是一个十分重要的问题，若步长h取的太小，将增加许多不必要的的计算量，若步长取的太大，其计算竞速很难保证，因此，与数值积分一样，微分方程的数值解也有一个选择步长的问题。由四阶龙格-库塔法的性质，其局部截断误差为这里的c与在内的值有关系。然后将步长折半，即取为步长，从跨两步到，在求一个近似值，每跨一步的误差为，因此，比较之后误差的大约减少。六计算时出现的问题，解决方法及体会： 对于初值问题，在用不同方法计算离散点的函数值的时候，有精确值给定的情况下，不同方法跟精确值的差别很大，也就是不是一种好方法，这与实际问题的稳定性有很大的关系，只要所求的初值问题，能够落在使用方法的稳定区间内，所得到的结果和精确值的离散点的函数值相差不会很大，并且在稳定区间内，改进欧拉法优于欧拉法，三阶之后的龙格-库塔方法优于改进欧拉法，离散点的函数值更加接近精确值。所谓的欧拉法就是一阶的龙格-库塔方法，改进欧拉法也是在显示欧拉法预测和梯形公式校正的基础上作出的改进，也是二阶龙格-库塔的特殊情况。所谓的三阶，四阶甚至n阶龙格-库塔方法，也是在待定系数有无数解得情况下的一般情况。 一阶常微分初值问题的收敛性，只要其增量函数关于y能满足利普希茨条件，这就收敛。而起在计算过程过有避免不了的误差，欧拉法精度相对来说比较低，龙格-库塔的阶数越高，精度就越高，所得的结果和精确值越接近，截断误差越小。同一个问题，想要减少截断误差，我们可以采用变步长的的方法，将步长变为，这样得出的结果，根据不同的方法，截断误差会有相应的改变。而在每种方法的程序设计方面，主要思想还是根据每种方法的公式，找到步长，分点数，用循环语句，得出相应的结果，只是在龙格-库塔方法里面的三阶，四阶的公式里面的待定系数有无数个解，我们给出的也是一般情况的三阶，四阶龙格-库塔方法，精度，误差相对于欧拉法，改进欧拉法，梯形公式都优于他们。数值计算实践课程设计原始数据记录实验名称：常微分方程初值问题数值解法研究 实验时间： 2013 年 1 月 14 日 输入：定义M文件并保存ffx.m Euler1(ffx,0,1,1,10)结果：T = 0 1.0000 0.1000 1.1000 0.2000 1.2100 0.3000 1.3310 0.4000 1.4641 0.5000 1.6105 0.6000 1.7716 0.7000 1.9487 0.8000 2.1436 0.9000 2.3579 1.0000 2.5937ans = Columns 1 through 7 1.0000 1.1000 1.2100 1.3310 1.4641 1.6105 1.7716 Columns 8 through 11 1.9487 2.1436 2.3579 2.5937向前欧拉公式结果：输入：Euler2(ffx,0,1,1,10)结果T = 0 1.0000 0.1000 1.1110 0.2000 1.2344 0.3000 1.3716 0.4000 1.5240 0.5000 1.6933 0.6000 1.8814 0.7000 2.0904 0.8000 2.3227 0.9000 2.5807 1.0000 2.8674ans = Columns 1 through 7 1.0000 1.1110 1.2344 1.3716 1.5240 1.6933 1.8814 Columns 8 through 112.0904 2.3227 2.5807 2.8674改进欧拉公式结果： 输入：x,y=Gaijineuler(f,0,1,1,0.1)结果：x = 0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000y = 1.0000 1.0959 1.1841 1.2662 1