您好,欢迎访问三七文档
当前位置:首页 > 办公文档 > 会议纪要 > 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是精确的。则),()()(phOxykhykh=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.梯形方法:1ky00)(),(),(yxyyxfyxf2hyykk1k1kk1k;3.改进的欧拉方法:预测值:),(^kkk1kyxfhyy校正值:1ky),(),(^kk1k1kkyxfyxf2hy。或改写为:)(21),(),(11cpkpkkckkkpyyyyxfhyyyxfhyy4、梯形方法与改进欧拉方法的截断误差是O(h3),具有二阶精度。5、龙格-库塔法的思想1).二阶龙格-库塔法计算公式:))((),(),(211121KK1hyyhKpyhpxfKyxfKkkkkkk当:21p时,得一簇龙格-库塔公式,其局部截断误差均为O(h3)都是二阶精度。特别取1p,21,就是改进欧拉公式。取1p,21,得二阶龙格—库塔法为:321121)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方法:][211kkkkff3hyy;三阶显式Adams方法:][2211kkkkkf5f16f231hyy;四阶显式Adams方法:]5559379[241231kkkkkkffffhyy.2).隐式Adams方法二阶隐式Adams方法:][211kkkkff3hyy;三阶隐式Adams方法:]8[2111kkkkkfff51hyy;四阶隐式Adams方法:)9195(241121kkkkkkffffhyy3).Adams预报-校正系统:先用显式格式作为预测值,再用隐式格式来校正。预测值:)9375955(243211kkkkkkffffhyy校正值:),9195(2411121kkkkkkkyxffffhyy4).改进的Adams预报-校正系统:预测:)9375955(243211kkkkkkffffhyp改进:)(kk1k1kpc270251pm校正值:)519,9(2421111kkkkkkkfffmxfhyc改进:)(1k1k1k1kpc27019cy7.收敛与稳定性对于固定的khxxk0,如果数值解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
本文标题:98常微分方程数值解法
链接地址:https://www.777doc.com/doc-4409028 .html