常微分方程数值解法

第八章 常微分方程的数值解法一．内容要点考虑一阶常微分方程初值问题:⎪⎩⎪⎨⎧==00)(),(y x y y x f dx dy微分方程的数值解：设微分方程的解y (x )的存在区间是[a,b ]，在[a,b ]内取一系列节点a= x 0< x 1<…< x n =b ,其中h k =x k+1-x k ；(一般采用等距节点，h=(b-a)/n 称为步长)。

在每个节点x k 求解函数y(x)的近似值:y k ≈y(x k )，这样y 0 , y 1 ,...,y n 称为微分方程的数值解。

用数值方法，求得f(x k )的近似值y k ，再用插值或拟合方法就求得y(x)的近似函数。

(一)常微分方程处置问题解得存在唯一性定理对于常微分方程初值问题:⎪⎩⎪⎨⎧==00)(),(y x y y x f dx dy如果：(1) 在B y y A x x 00≤-≤≤,的矩形内),(y x f 是一个二元连续函数。

(2) ),(y x f 对于y 满足利普希茨条件，即2121y y L y x f y x f -≤-),(),(则在C x x 0≤≤上方程⎪⎩⎪⎨⎧==00)(),(y x y y x f dxdy的解存在且唯一，这里C=min((A-x 0),x 0+B/L),L 是利普希茨常数。

定义：任何一个一步方法可以写为),,(h y x h y y k k k 1k Φ+=+，其中),,(h y x k k Φ称为算法的增量函数。

收敛性定理：若一步方法满足： (1)是p 解的.(2) 增量函数),,(h y x k k Φ对于y 满足利普希茨条件.(3) 初始值y 0是精确的。

则),()()(p h O x y kh y =-kh =x -x 0,也就是有0x y y lim k x x kh 0h 0=--=→)((一)、主要算法 1．局部截断误差局部截断误差：当y(x k )是精确解时，由y(x k )按照数值方法计算出来的1~+k y 的误差y (x k+1)- 1~+k y 称为局部截断误差。

注意：y k+1和1~+k y 的区别。

因而局部截断误差与误差e k +1=y (x k +1) -y k +1不同。

如果局部截断误差是O (h p+1),我们就说该数值方法具有p 阶精度。

1. 显式欧拉格式：⎩⎨⎧=+=+001)(),(y x y y x hf y y k k k k k =1,2,⋯n-1.显式欧拉格式的特点：(1)、单步方法； (2)、显式格式；(3)、局部截断误差)(2h O 因而是一阶精度 。

隐式欧拉格式⎩⎨⎧=+=+++00111)(),(y x y y x hf y y k k k k2. 两步欧拉格式：⎩⎨⎧=+=-+001)(),(y x y y x hf 2y y k k 1k k k =1,2,⋯n-1.两步欧拉格式的局部截断误差)(3h O ，因而是二阶精度 3. 梯形方法:=+1k y []⎪⎩⎪⎨⎧=+⋅+=+++00)(),(),(y x y y x f y x f 2h y y k k 1k 1k k 1k ； 3.改进的欧拉方法：预测值: ),(^k k k 1k y x f h y y ⋅+=+校正值: =+1k y ⎥⎦⎤⎢⎣⎡+⋅+++),(),(^k k 1k 1k k y x f y x f 2h y 。

或改写为:⎪⎪⎩⎪⎪⎨⎧+=+=+=++)(21),(),(11c p k p k k c k k k p y y y y x f h y y y x f h y y4、梯形方法与改进欧拉方法的截断误差是O(h 3),具有二阶精度。

5、龙格-库塔法的思想1). 二阶龙格-库塔法计算公式： ⎪⎩⎪⎨⎧λ+λ-+=++==+))((),(),(211121K K 1h y y hK p y h p x f K y x f K k k k k k k当:21=λp 时，得一簇龙格-库塔公式，其局部截断误差均为O (h 3)都是二阶精度。

特别取1p ==λ,21,就是改进欧拉公式。

取1p =λ=,21,得二阶龙格—库塔法为：⎪⎩⎪⎨⎧+=++==+21121)2,2(),(hK y y K h y h x f K y x f K k k kk k k ，称为二阶中点格式。

2)、经典龙格－库塔格式(也称为标准龙格－库塔方法)：⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧++++=++=++=++==+)22(6),()2,2()2,2(),(432113423121K K K K h y y hK y h x f K K h y h x f K K h y h x f K y x f K kk k k kk kk k k 四阶龙格－库塔方法的截断误差为()5h O ，具有四阶精度。

一般一阶常微分方程初值问题用四阶龙格－库塔方法计算，其精度均满足实际问题精度要求。

3).变步长龙格－库塔方法： 从节点x k 出发，以步长h 据四阶龙格－库塔方法求出一个近似值)(1h k y +,然后以步长h/2求出一个近似值)2/(1h k y +，得误差事后估计式：)()(1)2/(1)2/(1h k h k h k 1k y y 151y y ++++-≈- 根据)(1)2/(1h k h k y y ++-=∆来选取步长h 。

4). RKF 格式：变步长龙格－库塔方法，因频繁加倍或折半步长会浪费计算量。

Felhberg 改进了传统龙格－库塔方法，得到RKF 格式，较好解决了步长的确定，而且提高了精度与稳定性，为Matlab 等许多数值计算软件采用。

4/5阶RKF 格式：由4阶龙格－库塔方法与5阶龙格－库塔方法结合而成。

⎪⎭⎫⎝⎛+-+++=⎪⎭⎫⎝⎛-+++=++654311^5431155250956430285611282566561351651410121972565140821625K K K K K h y y K K K K h y y n n n n⎪⎪⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎪⎪⎨⎧-+-+-+=-+-++=+-++=+++=++==)401141041859256535442278,2()410484551336808216439,()219772962197720021971932,1312()329323,83()4,4(),(443216432153214213121K K K K K y h x f K K K K K y h x f K K K K y h x f K K K y h x f K K h y h x f K y x f K n n n n n nn n n n n nFelhberg 得到的最佳步长hs ，其中h 为当前步长，4/111^⎪⎪⎪⎪⎭⎫⎝⎛-ε=++n n y y h s .ε为精度要求，若s<0.75,步长折半；若s>1.5步长加倍。

6.亚当斯方法(Adams) 1).显式Adams 方法：记：),(k k k y x f f =；二阶显式Adams 方法：][211-+-+=k k k k f f 3hy y ;三阶显式Adams 方法：][2211--++-+=k k k k k f 5f 16f 231hy y ;四阶显式Adams 方法：]5559379[241231k k k k k k f f f f hy y +-+-+=---+.2). 隐式Adams 方法二阶隐式Adams 方法：][211++-+=k k k k f f 3hy y ;三阶隐式Adams 方法：]8[2111-++-++=k k k k k f f f 51hy y ;四阶隐式Adams 方法：)9195(241121+--+++-+=k k k k k k f f f f hy y3）.Adams 预报-校正系统：先用显式格式作为预测值，再用隐式格式来校正。

预测值： )9375955(243211---+-+-+=k k k k k k f f f f hy y校正值：()),9195(2411121++--+++-+=k k k k k k k y x f f f f hy y 4).改进的Adams 预报-校正系统：预测：)9375955(243211---+-+-+=k k k k k k f f f f hy p改进：)(k k 1k 1k p c 270251p m -+=++校正值：())519,9(2421111--++++-++=k k k k k k k f f f m x f hy c 改进：)(1k 1k 1k 1k p c 27019c y ++++--= 7.收敛与稳定性对于固定的kh x x k +=0，如果数值解y k 当,0h →（同时+∞→n ）时趋向于准确解y (x k )，则称该数值方法方法是收敛的。

如果一种数值方法，在节点值y k 上大小为δ的扰动，于以后各节点值y m (m >k )上产生的偏差均不超过δ，则称数值方法是稳定的.一般，隐式格式较显式格式有较好的稳定性。

8.刚性方程组：考虑n 阶常微分方程组：)(x Ay y φ+='，⋯⋯⋯（*）若矩阵A 的特征值λ1,λ2,⋯λn 的实部Re(λi )<0，i=1,2,⋯,n.，则称)Re()Re(i ni 1i n i 1min max S λλ=≤≤≤≤为*式的刚性比，当刚性比很大时称*式为刚性方程组。

刚性方程组需采用稳定性较好的方法。

二. 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-Moulton PECE)多步方法。

刚性系统：ode15s(Gear 方法)多步方法ode23s(二阶modified Rosenbroack formula)单步。

ode23t(trapezoidal rule)solve DAEs. 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;for n=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;for n=1:length(x)-1y(n+1)=iter(dyfun,x(n+1),y(n),h);endx=x';y=y';function y=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;while abs(y-y1)>ey1=y;y=y0+h*feval(dyfun,x,y);k=k+1;if k>K,error('迭代发散');endend改进Euler格式function[x,y]=naeulerg(dyfun,xspan,y0,h)x=xspan(1):h:xspan(2);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格式function[x,y]=nak4(dyfun,xspan,y0,h)x=xspan(1):h:xspan(2);y(1)=y0;for n=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)if nargin<5,h=(xspan(2)-xspan(1))/10;endn=1;x(n)=xspan(1);y(n)=y0;[y1,y2]=comput(dyfun,x(n),y(n),h);while x(n)<xspan(2)-eps;if abs(y2-y1)/10>ewhile abs(y2-y1)/10>eh=h/2;[y1,y2]=comput(dyfun,x(n),y(n),h);end elsewhile abs(y2-y1)/10<=e h=2*h;[y1,y2]=comput(dyfun,x(n),y(n),h); endh=h/2;h=min(h,xspan(2)-x(n)); [y1,y2]=comput(dyfun,x(n),y(n),h); endn=n+1;x(n)=x(n-1)+h;y(n)=y2;[y1,y2]=comput(dyfun,x(n),y(n),h); endx=x';y=y';function [y1,y2]=comput(dyfun,x,y,h); y1=rk4(dyfun,x,y,h); y21=rk4(dyfun,x,y,h/2);y2=rk4(dyfun,x+h/2,y21,h/2); function y=rk4(dyfun,x,y,h); K1=feval(dyfun,x,y);K2=feval(dyfun,x+h/2,y+h/2*K1); K3=feval(dyfun,x+h/2,y+h/2*K2); K4=feval(dyfun,x+h,y+h*K3); y=y+h*(K1+2*K2+2*K3+K4)/6;五.试验例题例1 不同方法的精度比较用多种方法解下述初值问题，并与其准确解x e y x +=-进行比较。