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

该【常微分方程初值问题数值解法的比较 】是由【wawa】上传分享，文档一共【16】页，该文档可以免费在线阅读，需要了解更多关于【常微分方程初值问题数值解法的比较 】的内容，可以使用淘豆网的站内搜索功能，选择自己适合的文档，以下文字是截取该文章内的部分文字，如需要获得完整电子版，请下载此文档到您的设备，方便您编辑和打印。2常微分方程初值问题数值解法的比较数值计算实践—课程设计报告课题名称常微分方程初值问题数值解法的比较完成时间2013-1-17姓名班级学号成果试验目的及内容1试验目的:(1)了解常微分方程初值问题的理论背景以及初值问题稳定性、收敛性的探讨;娴熟驾驭欧拉法、改进欧拉法、龙格-库塔法以及截断误差分析;比较欧拉法、改进欧拉法及龙格-库塔法,能够选择合适的方法进行问题的探讨计算;2试验内容:求微分方程(欧拉法求解)求微分方程(改进欧拉法求解)求微分方程(龙格-库塔求解)依据试验的结果进行分析,了解一般方法的的优缺点,稳定性,收敛性以及截断误差的分析,针对相应问题拿出有效方法得出最优的结果。:在科学与工程问题中,常微分方程表述物理量的变更规律,应用特别广泛,比如,天体运动的轨迹,机器人限制,化学反应过程的描述和限制以及电路瞬态过程分析等。这些问题中要求解随时间变更的物理量,即位置函数表示时间,而微分方程描述了未知函数与它的一阶或高阶导数之间的关系。考虑一阶常微分方程的初值问题,假如存在实数使得则称f关于y满意利普希茨条件,L称为利普希茨常数3。对于常微分方程初值问题,考虑初值的扰动是问题的解发生偏差的情形。若时的偏差被限制在有界范围内,则称该初值问题是稳定的,否则该初值问题不稳定的。特殊地,若时的偏差收敛于零,则称该初值问题是渐进稳定的。对于初值问题稳定性的探讨,易知其精确解为,假定初值经过扰动后变为,对于扰动后的解为因此带来的扰动误差为,因此考虑时的值,它取决于。易知,若,则原问题是稳定的;若,则原问题是不稳定的;若,则原问题是渐进稳定的。实际遇到的大多数常微分方程初值问题都是稳定的,因此在后面的探讨数值解法时这经常是默认条件。欧拉法:依据:积分曲线上一点的切线斜率等于函数值。方法:推动法,初始点动身,依照方向场在改点的方向推动到向前欧拉法的得到:(1)将在处泰勒绽开取h的线性部分,得(2)将初值问题中得导数用向前差商来代替有,因此(3)将两边同时对x的区间上积分对右端用左矩形公式得,此方法称向前欧拉法,也叫显示欧拉法。(4)对右端用右矩形公式得,也叫3隐式欧拉法。误差分析:;,那么称这种数值方法的阶数是p,其实p为非负整数。通常状况下,步长h越小,p越高,则局部截断误差越小,计算精初泰勒绽开有则有可见欧拉方法是一阶方法,精度不是很高。改进欧拉方法:梯形公式:对右端用梯形公式得+明显梯形公式是隐式公式。改进欧拉公式:先用欧拉公式求的一个初步的近似值,成为预料值,预料值的精度可能达不到要求,在用梯形公式将他校正一次,记为,这个结果成为校正值。预料:校正:误差分析:记为改进欧拉公式在处的截断误差,记因此有,表示在出的局部截断误差。由此得,梯形公式的局部截断误差为,因此改进欧拉的截断误差为,可见改进欧拉的方法是二阶方法,改进欧拉方法优于欧拉方法。龙格—库塔法:依据拉格朗日微分中值定理,,记得到,这样,只要给出一种计算的算法,就可以得到相应的计算公式。欧拉公式可以写为4改进欧拉公式可以写成因此推出一般的推出广式称为p阶龙格-库塔方法,简称p阶R-K方法。因为这里的均为常数。因为给定的系数不唯一,因此这里的常数有无穷多个解,下面是特殊状况下和一般状况下得结果(一阶龙格-库塔)当r=1时,这就是欧拉法。(二阶龙格-库塔)当r=2时,,这就是改进欧拉法。三阶和四阶龙格-库塔也只是在一般状况下得结果。三阶四阶其局部截断误差是::(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;x(1)=x0;fori=1:n6x(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';-库塔代码如下:(三阶龙格-库塔)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'];:输入:('ffx',0,1,1,10)7结果:T==:输入:Euler2('ffx',0,1,1,10)结果T==:输入:[x,y]=Gaijineuler('f',0,1,1,)结果:x==-库塔计算结果:(三阶)输入:Longgekuta3('ff',0,1,1,)结果:ans=Longgekuta3('ff',0,1,1,)ans=(四阶)输入:Longgekuta4('ff',0,1,1,)结果:ans=:方法显示欧拉简洁精度低隐式欧拉稳定性最好精度低,计算量大梯形公式精度提高计算量大中点公式精度提高,显式多一个初值,:若某算法对于随意固定的x=xi=x0+ih,当h?0(同时i?￥)时有yi?y(xi),则称该算法是收敛的。以下探讨的都是单步法(指在计算时只用到它前一步的信息)欧拉法,龙格-库塔法都是单步法的例子;例:就初值问题考察欧拉显式格式的收敛性。解:该问题的精确解为欧拉公式为?对于固定的x=xi=x0+ih,有设单步法具有p阶精度,满意利普希茨条件,其整体的截断误差:推断单步法的收敛性,归结为验证增量函数能否满意利普希茨条件。对于欧拉法,由于其增量函数就是f(x,y),因此当f(x,y)关于y满意利普希茨条件时候它就是收敛的。对于改进欧拉法和龙格-库塔方法,找到增量函数,满意利普希茨条件的时候,找到利普希茨常数,这些方法都是收敛的。:例:考察初值问题在区间[0,]上的解,分别用欧拉显,隐式格式,改进欧拉格式计算数值解:节点欧拉显式欧拉隐式改进欧拉法精确解