《常微分方程的数值解法分析解析》由会员分享，可在线阅读，更多相关《常微分方程的数值解法分析解析（32页珍藏版）》请在金锄头文库上搜索。
1、第八章常微分方程的数值解法内容要点考虑一阶常微分方程初值问题：dxf (x, y)y(xo) = yo微分方程的数值解：设微分方程的解y(x)的存在区间是a,b,在a,b内取一系列节 点a= xo xi xn =其中hk=xk+i-xk ；(一般采用等距节点，h=(b-a)/n称为步长)。在 每个节点xk求解函数y(x)的近似值：y- y(xk),这样yo, yi,.,yn称为微分方程的数值解。用数值方法，求得f(xk)的近似值y.再用插值或拟合方法就求得 y(x)的近似函数。(一)常微分方程处置问题解得存在唯一性定理对于常微分方程初值问题:dXH(x, y)如果：y(xo).。(1)在x。，
2、人,卜%的矩形内f (x, y)是-个二元连续函数。(2) f (x, y)对于y满足利普希茨条件，即f (x, yi) (x, y2 )|L|yi21则在x。|x|c上方程.=f (x, y)的解存在且唯-,这 y(xo)3y0里C=min(A-xo),xo+B/L),L是利普希茨常数。定义：任何一个一步方法可以写为y|Jyk(xk,yk,h),其中同(xk,yk,h)称为算法 的增量函数。收敛性定理：若一步方法满足：是p解的.(2)增量函数图x,y.h)对于y满足利普希茨条件.(3)初始值yo是精确的。则|y(kh)，y(x)倡O(hp), kh=x-xo，也就是有hlim，k.y(x)号
3、kh -x - xQ(一)、主要算法1.局部截断误差局部截断误差：当y(xk)是精确解时，由y(xk)按照数值方法计算出来的丫息的误差y(xk+i)- yk称为局部截断误差。注意：yk+i和y届i区别。因而局部截断误差与误差ek+i=y(xk+i)-yk+i不同。如果局部截断误差是O(hp+i),我们就说该数值方法具有p阶精度1 .显式欧拉格式:，k 1=yk hf (xk, yjk=1,2, n-1.y(xo) = yo显式欧拉格式的特点：(1)、单步方法；(2)、显式格式；(3)、局部截断误差O(h2)因而是一阶精度隐式欧拉格式yk hf(Xk i,yk i) y(%) = y。2 .两步
4、欧拉格式:k 1 =%2hf(Xk，yk)k=1,2, n-1.y(X。) = y。两步欧拉格式的局部截断误差 o(h3),因而是二阶精度3 .梯形方法：yk 1,ykyk h f (Xk 13.改进的欧拉方法:预测值：yky(xo) = yoykh f (Xk, yk)校正值:ykykhyp BykhDf(Xk,yk)f(Xk i,yk i) . f (Xk,yk)或改写为ICBykhBf(xlllyp)1 .、yk i 二2(yp yJ4、梯形方法与改进欧拉方法的截断误差是 O(h3),具有二阶精度。5、龙格-库塔法的思想1).二阶龙格-库塔法计算公式:K1. (Xk, yk)心(XkHp
5、h, yk 3phK1)yk i时，得一簇龙格-库塔公式，其局部截断误差均为 O(h3)都是二阶精度。特别 2取厦1, P.1，就是改进欧拉公式。2取p =1,得二阶龙格一库塔法为:Ki = f (x.yk)K2Bf(xkBh,ykHhK1),称为二阶中点格式。yk hK22)、经典龙格-库塔格式(也称为标准龙格-库塔方法)：Ki = f(Xk, yk)h hK2 = f (Xk, yk Ki)22K3 = f (Xk h, yk h.) 22K4 = f (xh, yh.)yk h(K1 2K2 2K3 K4)6四阶龙格一库塔方法的截断误差为h5 I,具有四阶精度。值问题用四阶龙格-库塔方法
6、计算，其精度均满足实际问题精度要求般一阶常微分方程初3) .变步长龙格库塔方法：从节点Xk出发，以步长h据四阶龙格库塔方法求出一个近似值ykj|然后以步长h/2求出一个近似值y|,得误差事后估计式：(h/2)1 (h/2)(h)yk - yk 1 .(yk i -丫)根yllli y|来选取步长h。4) . RKF格式：变步长龙格库塔方法，因频繁加倍或折半步长会浪费计算量。Felhberg 改进了传统龙格-库塔方法，得到 RKF格式，较好解决了步长的确定，而且提高了精 度与稳定性，为Matlab等许多数值计算软件采用。4/5阶RKFB式：由4阶龙格-库塔方法与5阶龙格-库塔方法结合而成。252
7、16K1 1408 K321972565410116135ayy6656.1282528561 K3K56430K4 -555J50K1 - f(xn,yn )K2 =f(Xn Lyn 1) 44K3=f(xn 8h，yn 329K132K2)12h1932K4 = f(Xn ITyn -7200 K12197K2 K2197K5 = f (Xn h,yn 4393680K18K2 巫0216513f(xn h,yn-_8Kl 2K2 - 3544 K 322725658454104K3)K61859410411K4-40K4)1/4为精度要求,y二阶显式Adams方法:三阶显式Adams方法
8、:四阶显式Adams方法:hyk 1 = yk 23fk - fk,h .yk 1 = yk 12 23 fk-16fk5fk/;yk 1 = y2h4-9 fkj3 37 fk/ 一59 fk55 fk.Felhberg得到的最佳步长hs,其中h为当前步长,若s1.5步长加倍6.亚当斯方法(Adams)1) .显式 AdamspJ法:记：fkHf(xk, yk);三阶隐式Adams方法:yk 23fk - fk 1;yk .5fk 1 8fk - fk四阶隐式Adams方法:yk 1 = y24(fk -5f19fk 9fk1)2) .隐式Adams方法二阶隐式Adams方法:3) .Ada
9、ms预报-校正系统：先用显式格式作为预测值，再用隐式格式来校正。预测化 yk1 = yk ，(55fk -59fk37f_ -9f-)校正值:ykLyk 24(f-5f19fk 94).改进的Adams预报-校正系统:改进：m预测：Pk校正值:改进：yk 1(ck ii)0,(同时fi 一 Pk i)时趋向于准确解y(x.，(55fk .59fk4 37 fy _9f.),m251 ,、(ck - pk)270Ck 1 = y,9 f xk 119ck 12707.收敛与稳定性对于固定的xk.x0 Bkh ,如果数值解yk当19fk -5 f则称该数值方法方法是收敛的。ym(mk)上产生的如果
10、一种数值方法，在节点值 yk上大小为事扰动，于以后各节点值 偏差均不超过、则称数值方法是稳定的.8.刚性方程组：考虑n阶常微分方程组:)n的实部Re(i)eyi=y;y=y0+h*feval(dyfun,x,y);k=k+1;if kK,error(迭代发散);end end 改进Euler格式functionx,y=naeulerg(dyfun,xspan,y0,h) x=xspan:h:xspan； y(1)=y0;for n=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格式functionx,y=nak4(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2);y(1)=y0;for n=1:length
