98常微分方程数值解法

1第八章常微分方程的数值解法一．内容要点考虑一阶常微分方程初值问题:00)(),(yxyyxfdxdy微分方程的数值解：设微分方程的解y(x)的存在区间是[a,b]，在[a,b]内取一系列节点a=x0x1…xn=b,其中hk=xk+1-xk；(一般采用等距节点，h=(b-a)/n称为步长)。在每个节点xk求解函数y(x)的近似值:yk≈y(xk)，这样y0,y1,...,yn称为微分方程的数值解。用数值方法，求得f(xk)的近似值yk，再用插值或拟合方法就求得y(x)的近似函数。(一)常微分方程处置问题解得存在唯一性定理对于常微分方程初值问题:00)(),(yxyyxfdxdy如果：(1)在ByyAxx00,的矩形内),(yxf是一个二元连续函数。(2)),(yxf对于y满足利普希茨条件，即2121yyLyxfyxf),(),(则在Cxx0上方程00)(),(yxyyxfdxdy的解存在且唯一，这里C=min((A-x0),x0+B/L),L是利普希茨常数。定义：任何一个一步方法可以写为),,(hyxhyykkk1k，其中),,(hyxkk称为算法的增量函数。收敛性定理：若一步方法满足：(1)是p解的.(2)增量函数),,(hyxkk对于y满足利普希茨条件.(3)初始值y0是精确的。则),()()(phOxykhykh=x-x0,也就是有0xyylimkxxkh0h0)((一)、主要算法1．局部截断误差局部截断误差：当y(xk)是精确解时，由y(xk)按照数值方法计算出来的1~ky的误差y(xk+1)-1~ky称为局部截断误差。注意：yk+1和1~ky的区别。因而局部截断误差与误差ek+1=y(xk+1)-yk+1不同。如果局部截断误差是O(hp+1),我们就说该数值方法具有p阶精度。21.显式欧拉格式：001)(),(yxyyxhfyykkkkk=1,2,n-1.显式欧拉格式的特点：(1)、单步方法；(2)、显式格式；(3)、局部截断误差)(2hO因而是一阶精度。隐式欧拉格式00111)(),(yxyyxhfyykkkk2.两步欧拉格式：001)(),(yxyyxhf2yykk1kkk=1,2,n-1.两步欧拉格式的局部截断误差)(3hO，因而是二阶精度3.梯形方法:1ky00)(),(),(yxyyxfyxf2hyykk1k1kk1k；3.改进的欧拉方法：预测值:),(^kkk1kyxfhyy校正值:1ky),(),(^kk1k1kkyxfyxf2hy。或改写为:)(21),(),(11cpkpkkckkkpyyyyxfhyyyxfhyy4、梯形方法与改进欧拉方法的截断误差是O(h3),具有二阶精度。5、龙格-库塔法的思想1).二阶龙格-库塔法计算公式：))((),(),(211121KK1hyyhKpyhpxfKyxfKkkkkkk当:21p时，得一簇龙格-库塔公式，其局部截断误差均为O(h3)都是二阶精度。特别取1p,21,就是改进欧拉公式。取1p,21,得二阶龙格—库塔法为：321121)2,2(),(hKyyKhyhxfKyxfKkkkkkk，称为二阶中点格式。2)、经典龙格－库塔格式(也称为标准龙格－库塔方法)：)22(6),()2,2()2,2(),(432113423121KKKKhyyhKyhxfKKhyhxfKKhyhxfKyxfKkkkkkkkkkk四阶龙格－库塔方法的截断误差为5h，具有四阶精度。一般一阶常微分方程初值问题用四阶龙格－库塔方法计算，其精度均满足实际问题精度要求。3).变步长龙格－库塔方法：从节点xk出发，以步长h据四阶龙格－库塔方法求出一个近似值)(1hky,然后以步长h/2求出一个近似值)2/(1hky，得误差事后估计式：)()(1)2/(1)2/(1hkhkhk1kyy151yy根据)(1)2/(1hkhkyy来选取步长h。4).RKF格式：变步长龙格－库塔方法，因频繁加倍或折半步长会浪费计算量。Felhberg改进了传统龙格－库塔方法，得到RKF格式，较好解决了步长的确定，而且提高了精度与稳定性，为Matlab等许多数值计算软件采用。4/5阶RKF格式：由4阶龙格－库塔方法与5阶龙格－库塔方法结合而成。654311^5431155250956430285611282566561351651410121972565140821625KKKKKhyyKKKKhyynnnn)401141041859256535442278,2()410484551336808216439,()219772962197720021971932,1312()329323,83()4,4(),(443216432153214213121KKKKKyhxfKKKKKyhxfKKKKyhxfKKKyhxfKKhyhxfKyxfKnnnnnnnnnnnn4Felhberg得到的最佳步长hs，其中h为当前步长，4/111^nnyyhs.为精度要求，若s0.75,步长折半；若s1.5步长加倍。6.亚当斯方法(Adams)1).显式Adams方法：记：),(kkkyxff；二阶显式Adams方法：][211kkkkff3hyy;三阶显式Adams方法：][2211kkkkkf5f16f231hyy;四阶显式Adams方法：]5559379[241231kkkkkkffffhyy.2).隐式Adams方法二阶隐式Adams方法：][211kkkkff3hyy;三阶隐式Adams方法：]8[2111kkkkkfff51hyy;四阶隐式Adams方法：)9195(241121kkkkkkffffhyy3）.Adams预报-校正系统：先用显式格式作为预测值，再用隐式格式来校正。预测值：)9375955(243211kkkkkkffffhyy校正值：),9195(2411121kkkkkkkyxffffhyy4).改进的Adams预报-校正系统：预测：)9375955(243211kkkkkkffffhyp改进：)(kk1k1kpc270251pm校正值：)519,9(2421111kkkkkkkfffmxfhyc改进：)(1k1k1k1kpc27019cy7.收敛与稳定性对于固定的khxxk0，如果数值解yk当,0h（同时n）时趋向于准确解y(xk)，则称该数值方法方法是收敛的。如果一种数值方法，在节点值yk上大小为的扰动，于以后各节点值ym(mk)上产生的偏差均不超过，则称数值方法是稳定的.5一般，隐式格式较显式格式有较好的稳定性。8.刚性方程组：考虑n阶常微分方程组：)(xAyy，（*）若矩阵A的特征值1,2,n的实部Re(i)0，i=1,2,,n.，则称)Re()Re(ini1ini1minmaxS为*式的刚性比，当刚性比很大时称*式为刚性方程组。刚性方程组需采用稳定性较好的方法。二.MATLAB的有关命令[t，x]=solver（’f’,ts,x0,options）t:输出的自变量值,x:输出的函数值；solver:包括ode45、ode23、ode113、ode15、sode23s.非刚性系统：ode23：用2阶(及3阶)龙格-库塔算法。ode45：用4阶(及5阶)龙格-库塔-算法。ode113(Adams-Bashforth-MoultonPECE)多步方法。刚性系统：ode15s(Gear方法)多步方法ode23s(二阶modifiedRosenbroackformula)单步。ode23t(trapezoidalrule)solveDAEs.ode23tb(TR-BDF2)低精度算法。f：由待解方程写成的m-文件名ts=[t0，tf]，t0、tf为自变量的初值和终值x0：函数的初值options=odeset(’reltol’,rt,’abstol’,at),rt，at：分别为设定的相对误差和绝对误差.(缺省时设定相对误差10-3,绝对误差10-6)注意:1、在解n个未知函数的方程组时，x0和x均为n维向量，m-文件中的待解方程组应以x的分量形式写成.2、使用Matlab软件求数值解时，高阶微分方程必须等价地变换成一阶微分方程组.四.重点算法Matlab程序Euler格式、function[x,y]=naeuler(dyfun,xspan,y0,h)x=xspan(1):h:xspan(2);y(1)=y0;forn=1:length(x)-1y(n+1)=y(n)+h*feval(dyfun,x(n),y(n));endx=x';y=y';隐式Euler格式function[x,y]=naeulerb(dyfun,xspan,y0,h)x=xspan(1):h:xspan(2);y(1)=y0;6forn=1:length(x)-1y(n+1)=iter(dyfun,x(n+1),y(n),h);endx=x';y=y';functiony=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;whileabs(y-y1)ey1=y;y=y0+h*feval(dyfun,x,y);k=k+1;ifkK,error('迭代发散');endend改进Euler格式function[x,y]=naeulerg(dyfun,xspan,y0,h)x=xspan(1):h:xspan(2);y(1)=y0;forn=1:length(x)-1K1=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;endx=x';y=y';4阶经典Runge-Kutta格式function[x,y]=nak4(dyfun,xspan,y0,h)x=xspan(1):h:xspan(2);y(1)=y0;forn=1:length(x)-1K1=feval(dyfun,x(n),y(n));K2=feval(dyfun,x(n)+h/2,y(n)+h/2*K1);K3=feval(dyfun,x(n)+h/2,y(n)+h/2*K2);K4=feval(dyfun,x(n+1),y(n)+h*K3);y(n+1)=y(n)+h*(K1+2*K2+2*K3+K4)/6;endx=x';y=y';变步长4阶经典Runge-Kutta格式function[x,y]=nark4v(dyfun,xspan,y0,e,h)ifnargin5,h=(xspan(2)-xspan(1))/10;endn=1;x(n)=xspan(1);y(n)=y0;[y1,y2]=comput(dyfun,x(n),y(n),h);whilex(n)xspan(2)-eps;ifabs(y2-y1)/10ewhileabs(y2-y1)/10eh=h/2;[y1,y2]=comput(dyf