常微分方程数值解法教材.docx
- 文档编号:28270646
- 上传时间:2023-07-10
- 格式:DOCX
- 页数:39
- 大小:143.89KB
常微分方程数值解法教材.docx
《常微分方程数值解法教材.docx》由会员分享,可在线阅读,更多相关《常微分方程数值解法教材.docx(39页珍藏版)》请在冰豆网上搜索。
常微分方程数值解法教材
常微分方程数值解法
科学技术中的许多问题,在数学中往往归结为微分方程的求解问题。
例如天文学中研究星体运动,空间技术中研究物体飞行等,都需要求解常微分方程初值问题。
最简单的常微分方程初值问题是
座=f(x,y),a£x£b,
dx'"(*)
y(xo)=y°.
定理如果f(x,y)在区域D={(x,y)a《x〈b,y<°°}内连续,且关于y满足Lipschitz条件,那么初值问题(*)的解存在且惟一。
因为数值解是求微分方程解y(x)的近似值,所以总假定微分方程的解存在惟一,即初值
问题(*)中的f(x,y)满足定理的条件。
除特殊情形外,初值问题(*)一般求不出解析解,即使有的能求出解析解,其函数表示式
也比较复杂,计算量比较大,况且实际问题往往只要求在某一时刻解的函数值,因此,有必要研究初值问题(*)的数值解法。
所谓初值问题(*)的数值解法,就是寻求解y(x)在区间[a,b]上的一系列点
x〔:
:
:
x^:
:
x3:
:
:
..•:
:
:
xn:
:
:
■■-
上的近似值yi,y2,…,yn,….记hj=x—x口(i=1,2,|||)表示相邻两个节点的间距,称为步长。
求微分方程数值解的主要问题:
(1)如何将微分方程也=f(x,y)离散化,并建立求其数值解的递推公式;
dx
(2)递推公式的局部截断误差、数值数yn与精确解y(xn)的误差估计;
(3)递推公式的稳定性与收敛性.
、Euler方法
Euler法是最早的解决一阶常微分方程初值问题的一种数值方法,虽然它不够精确,很少被采
用,但是它在某种程度上反映了数值方法的基本思想和特征。
考虑初值问题
%f(x,y),(1.1)
dx
yMi。
,(1.2)
为了求得它在等距离散点x1 设 h=Xi—x—(i=1,2,|||),将式(1.1)离散化的办法有Taylor展开法、数值微分法及数值积分法。 如果在点xn将y(xn+^y(x^+h)作Taylor展开,得 h2 (1.3) y(Xn.1)=y(Xn)hy(Xn);Y(n),n(Xn,Xn1) 那么当h充分小时,略去误差项}y”(知),用yn近似替代y(xn)、yn41近似替代y(xn书), (1.4) 并注意到y'(Xn)=f(Xn,y(Xn)),便得y(o=y(x()),yn1=ynhf(Xn,yn),n=0,1,山,N-1b-a 其中Xn=x。 +nh,h=z「.用(1.4)求解(1.1)的万法称为Euler万法。 如果利用差商近似替代微商,那么可得 y(Xn1)-yg)y(Xn)=f(Xn,y(Xn)).(1.5)h 在(1.5)中若用yn近似替代y(xn)、yn中近似替代y(xn41),同样得到递推公式(1.4). 如果在[xn,xn书]上对y「=f(x,y(x))积分,得 y(Xn*)—y(Xn)=L: +f(X,y(x))dx.(1.6) 那么对上式右端积分用左矩形求积公式,并用yn近似替代y(xn)、yn虫近似替代y(xn虫),也 可得到递推公式(1.4). 我们知道,在xOy平面上,微分方程业=f(x,y)的解称为积分曲线,积分曲线上一点dx (x,y)的切线斜率等于函数f(x,y)的值。 如果在D中每一点,都画上一条以f(x,y)在这点的值为斜率并指向x增加方向的有向线段,即在D上作出了一个由方程d^=f(x,y)确定的方向场,那么方程的解f=y(x).从几何上看,就是位于此方向场中的曲线,它在所经过的每 一点都与方向场的该点的方向相一致。 y=y(x),斜率为f(x0,y0).设在x=x0 y=y°+f(X0,y°)(x—x°).当x=x〔时, y1,这就是得到x=x1时计算y(x1)的近 从初始点P0(x0,y0)出发,过这点的积分曲线为附近y(x)可用过P0点的切线近似表示,切线方程为y(x1)的近似值为y0+f(x0,y0)(x1—x0),并记为 似公式 y=y〔f(*,山)(*一为). 当x=x2时,y(x2)的近似值为y1+f(x1,y1)(x2—x1),并记为y2.于是就得到当x=x2时讨 算y(x2)的近似公式 V2=*f(X1,y〔)(x2—x) 重复上面方法,一般可得当x=xn土的计算y(xn*)的近似公式 Yn^Ynf(Xn,yn)(Xni*. 如果h=x—X-(i=1,2,III),则上面公式就是(1.4).将P°,Pi,…,Pn连续起来,就得到一条 折线,所以Euler方法又称为折线法。 由公式(1.4)看出,已知y0便可算出Y1.已知y[,便可算出y2,如此继续下去,这种只用前一步的值yk便可计算出yk*的递推公式称为单步法。 Yn替代y(Xn),Yhi替 若在(1.6)中,其右边的积分由数值积分的右矩形公式近似,并用 代y(xnq,则可得到 (1.7) §0=)(处), Yn*=Yn+hf(Xn为Yn+), 并称(1.7)为后退Euler公式。 Euler公式(1.4)是关于yn*的一个直接计算公式,是显式的;而公式(1.7)右端是含有yn中 的一个函数方程,因此是隐式的,也是单步法。 图9-1 若在公式(1.6)中,其右边积分用数值梯形积分公式近似,并用yn替代y(xn),y^H替代f(Xn*),则可得到梯形方法公式 hYn1=Yn: [f(Xn,Yn)f(Xn1,Yn帛(1.8)2 梯形方法同后退Euler方法一样都是隐式单步法.对于隐式方法,通常采用迭代法。 对后退Euler方法,先Euler方法计算y^,并将它作为初值y: %,即y;%=Yn+hf(Xn,Yn), 再把它代入(1.7)的右端,便得到后退Euler方法的迭代公式为 (1.9) (1.10) jYn%=Ynff(Xn,yn),(/项=Yn+hf(Xn*,yn%),k=0,1,…. 同样地,仍用Euler方法提供初始值,梯形方法的迭代公式为 yn°1=Ynhf(Xn,Yn), 'yn")=Yn+*f(Xn,Yn)+f(XmYnW'k=0,1, 二、误差分析 首先讨论Euler方法的误差.不妨假定在计算yn书时用到前面一步的值是准确值y(Xn),即yn=y(Xn).利用算法由y(Xn)计算出的y(Xn+)的近似值为y^H,则尸(2书)一贝书称为局部截断误差,也就是计算一步所产生的误差,记为Tn书. 将y(Xn+)在Xn处作Taylor展开 h2 y(Xn.1)=y(Xnh)=y(Xn)hf(Xn,y(Xn))wy(),Xn: : : ,: Xn4. 由Euler方法(1.4)得 ~Ii_r/\/\Ii_r//\\ Vn4=Ynhf(Xn,yn)=y(Xn)hf(^,y(Xn)). 上面两式相减得 Tn^=yg*)-"书=? 广(令).(1.11) 若v(x)在解的存在区间[a,b]上充分光滑,且记M2=maXy"(X),则Euler方法局部截X^a,b]J 断误差为 h22 Tn1£M25=°(h)(1.12) 由于计算y^1时(n=0除外),一般用到的不是y(Xn),而是近似值yn,因为每一步计算除局部截断误差外,还应考虑前一步不准确而引起的误差。 我们称这种误差为总体截断误差 或方法误差。 体截断误差有一个积累过程,它与计算步数有关。 把第n步的总体截断误差记为en=y(Xn)-yn.对任意的第n+1步有 en+=y(Xn+)—yn]Wy(Xnf)-Vn/+Vn*-y"由(二11)知yg*)-]丽,故只须估计反尚一尸"方. yn1-yn1=y(Xn)hf(Xn,y(Xn))-y”-hf(Xn,y”)Wy(Xn)-yn|+hf(Xn,y(Xn))-f(Xn,yn). 因为f(x,v)关于y满足Lip条件,所以 |,如+-yn^K+hL*=(1+hL)en 这里L是Lip常数。 综上所述 |eJg+(1+hL)斜(1.13) (1.13)给出了第n+1步总体截断误差与第n步总体截断误差之间的关系,它对一对n都成立。 特别当n=0时,eJ eo=y(x0)—y0=0,因此 en (1hL)e~M|Tn|+(1+hL)[|Tn』+(1+hL)|en项=Tn(1hL)TnJ(1hL)2en/-HI〈|Tn|+(1+hL)Tn」+(1+hL)2Tu+H|+(1+hL尸T1 斜|苴£(1+hL)kTj=O(h2)! (1+hL)k=(: **): 1O(h2)=O(h).(1.14) 心1+hL—1 这说明在区间[x0,xN]上用Euler方法求初值问题(1.1)的数值解,若取步长h=Xn一X0,则在 N Xn(n=1,…,N)各点上的总体截断误差为en=O(h). 这表明当hT0时,ent0,也就是说当h充分小时近似值yn能和y(xn)充分接近,即 数值解是收敛的。 对后退Euler方法,类似上面的讨论可知,其局部截断误差为 h23 y(Xn1)-yn『"y(Xn)O(h), 即局部截断误差关于h是二阶的,整体截断误差关于h是一阶的. 梯形方法公式是将(1.6)中右边积分用梯形求积公式得到的,而梯形方法公式的截断误 h3 差为—一广(豪),Xn<豪<Xz,因此梯形方法公式(1.8)的局部截断误差关于h是三阶的, 12 总体截断误差关于h是二阶的。 可以证明: 对后退Euler方法的迭代公式(1.9),如果h充分小,使得0<hL<1,那么迭代过程是收敛的;对梯形方法的迭代公式(1.10),如果h充分小,使得0〈里<1,那么迭代过 2 程也是收敛的(这里L是Lip常数)。 三、改进的Euler方法 梯形方法的迭代公式(1.10)比Euler方法精度高,但其计算较复杂,在应用公式(1.10)进行计算时,每迭代一次,都要重新计算函数f(x,y)的值,且还要判断何时可以终止或转下一步计算. 为了控制计算量和简化计算法,通常只迭代一次就转入下一步计算.具体地说,我们先用Euler公式求得一个初步的近似值Vn+,称之为预测值,然后用公式(1.10)作一次迭代得yn4i,即将yn书校正一次.这样建立的预测一校正方法称为改进的Euler方法: 预测: yn^ynhf(Xn,yn), h 校正: yn1=yn2[f(Xn,yn)f(Xn1,yn1)].(1.15) 这个计算公式也可以表示为 r yp=yn+hf(Xn,yn), yc=ynhf(Xn1,yp), 1,,、 =」yp+y) 2 例1取步长h=0.1,分别用Euler方法及改进的Euler方法求解初值问题 手=一y(1xy),0£x5,y(0)=1. 解这个初值问题的准确解为y(x)=1/(2eX—x-1).根据题设知f(x,y)=—y(1+xy). (1)Euler方法的计算式为 Vn1二尾一°.1[yn(14忘)], 由y°=y(°)=1,得 y1=1-0.1[1(101)]=0.9, y2=0.9-0.1[0.9(10.10.9)]=0.8019, 这样继续计算下去,其结果列于表9.1. ⑵改进的Euler方法的计算式为 yp=yn—0.1K[yn(1+Xnyn)], Vg=yn-。 .1[yp(1Xndyp)], '1,,、 yn+=;(yp成), 2 由y0=y(0)=1,得 yp=1-0.1尺[1乂(1+0乂1)]=0.9,yc=1-0.1[0.9(10.10.9)]=0.9019, 1. y〔=? (0.9+0.9019)=0.90095 f yp=0.90095-0.1乂[0.90095乂(1+0.1又0.90095)]=0.80274, yc=0.90095-0.1[0.80274(V0.20.80274)]=0.80779, 1工 y2=—(0.80274+0.80779)=0.80526 2 这样继续计算下去,其结果列于表 Euler方法 改进的Euler方法 准确值 Xn yn yn y(Xn) 0.1 0.9000000 0.9009500 0.9006235 0.2 0.8019000 0.8052632 0.8046311 0.3 0.7088491 0.7153279 0.7144298 0.4 0.6228902 0.6325651 0.6314529 0.5 0.5450815 0.5576153 0.5563460 0.6 0.4757177 0.4905510 0.4891800 0.7 0.4145675 0.4310681 0.4296445 0.8 0.3610801 0.3786397 0.3772045 0.9 0.3145418 0.3326278 0.3312129 1.0 0.2741833 0.2923593 0.2909884 从表9.1可以看出,Euler方法的计算结果只有2位有效数字,而改进的Euler方法确有3位 有效数字,这表明改进的Euler方法的精度比Euler方法高. §2Runge-Kutta方法 上一节介绍的Euler方法、梯形方法及改进的Euler方法都是单步法,即计算yn+i时,只 2、 用到yn.Euler方法及后退的Euler方法的局部截断快差是0(h),总体截断快差是0(h),梯形方法的总体截断误差是0(h2).一个方法的总体截断误差若为0(hp),贝U我们称它为p阶方法。 体截断误差和局部截断误差的关系是: 总体截断误差=03。 乂局部截断误差 一般地,方法的总体截断误差阶越高,该方法的精度也越高。 在上一节我们利用y(Xnt)在点X。 的Taylor展开的前两项导出Euler公式.很自然地,我们 想到: 若将y(xn+)在点x0的Taylor展开多取几项,则有希望获得高阶方法。 但是直接利用 Taylor展开取多项的办法,需要计算f(x,y)的高阶导数,计算量较大。 因此在建立Runge-Kutta方法时不采用求高阶导数的办法,而是用计算不同点上的函数值,然后对这些值作线性组合, 构造近似公式,把近似公式和解的Taylor展开相比较,使它们在前面的若干项相同,从而使 近似公式达到一定的阶。 二阶Runge-Kutta方法 我们已经知道Euler方法每一步只计算一次f(x,y)的值,总体截断误差为0(h),它的计 算式可以写成 ;y2=yn+h&(y。 已知) 0(h2),它的计算公式 &=f(Bn) 而改进的Euler方法,每迭代一步要计算两个函数的值,其总体误差为可改写成 hhyn*=yn+二名+匚k2,22』临=f(xn,yn),k2=f(xnh,ynhki). 以计算两次函数值为例,设一般计算式为 j_yn1=ynh(c1k1C2k2), *1=f(xn,yn),(2.1) “2=f(xn+ph,yn+qhki). 适当选择参数G,C2,p,q的值,使得在y(xn)=yn的假设下,截断误差y(xz)-y^H阶尽可 能高。 为此,将yn*=yn+h(Ciki+C2k2)在点(xn,yn)作Taylor展开,因为 ki=f(xn,y(xn))=y(xn), k2的展开式为 所以 yn1=y(Xn)(CC2)hf(Xn,y”) C2ph2[fx(Xn,yn)Qfy(Xn,yn)f(Xn,yn)]O(h3). P 再将y(X")在x=Xn作Taylor展开 y(Xn1)=y(Xn)hy(Xn)}[fx(Xn,yn)fy(Xn,Yn)f(Xn,Yn)]O(h3). 若要求局部截断误差达到 O(h3),则通过比较上面两式知,参数Ci,C2,p,q必须满足 一一一_1…1 c1'C2=1,c2p=—,c2q=22 、……-•、…,-、一一'“、I…一1 这是四个未知量,三个方程式的方程组,因此解不惟一。 当我们选取G=c2=」,p=q 2 (2.2)和(2,.3)的前三项完全一致。 将参数值代入(2.1)便得到二阶Runge-Kutta R-K方法)的一种迭代方法 hh,yn*=yn十二k〔十二k2,22*=f(Xn,yn),k2=f(Xn+h,yn+hk〔), (2.2) (2.3) (2.4) 它的局部截断误差为O(h2),(2.4)实际上就是改进的Euler方法(1.15). …1 右取q=0,c2=1,p=q,则碍到一阶R-K方法的又一种迭代公式 2 yn*=yn+hk2, 』k1=f(Xn,yn), k2=f(Xn+£,yn+%), L22 或写成 …h1…、、 yn1=ynhf(Xn匚,yn: hf(Xn,yn)).(2.5) 22 公式(2.4)、(2.5)都是显式单步二阶公式. 可以证明了不论这四个参数如何选择,都不能使局部截断误差达到O(h4).这说明在计算 两次函数值的情况下,局部截断误差的阶最高为O(h3),要再提高阶就必须增加计算函数值 的次数。 二、三阶及四阶Runge-Kutta方法 为了提高方法的精度,考虑每步计算三次函数f(x,y)值。 根据两次计算函数值的做法,很 自然地取据+的形式为 Yn1=ynh(Ciki弘2块), J=f(Xn,Yn),k2=f(Xna〔h,Ynb21hk〔),k3=f(Xna2h,Ynb31hk1bazhk? ). 适当选择参数c、G、G、a、&、a3、质、Qi、bj2,使截断误差y(xn*)—y忡的阶尽可能高。 由于在推导公式时只考虑局部截断误差,故设yn=y(xn)。 类似前面二阶r-k公式的推 导方法,将yz在(xn,yn)处作Taylor展开,然后再将y(xnG在X=xn处作Taylor展开(这里不详细写出展开式),只要两个展开式的前四项相同,便有y(xn噂)-yn41=0(h4),而要两 个展开式的前四项相同,参数必须满足 c1c2c3~1,1 职2+a2%=—, 2 这是8个未知数6个方程的方程组,解不是惟一的,可以得到很多公式。 当我们取方程组的 1一1一一-141..„.组解为a1=—,a2=1,b21=—,b31=-1,b32=2,c1—,c2=,c3=-时,就碍到一 22666 种常用的三阶Runge-Kutta公式 hyn*=yn+二(崎+4k2+k3),6 (2.7) 、=f(Xn,yn), k2=f(Xn+}yn+? 峪), k3=f(Xn+h,yn-hk1+2hk2). 从上面导出公式(2.7)的过程知,它的局部截断误差为O(h4). 如果每步计算四次函数f(x,y)的值,完全类似的,可以导出局部截断误差为O(h5)的四 公式(2.8)就是一种常用的四阶Runge-Kutta公式,也称为标准的(或经典的)四阶R-K方法。 三阶、四阶Runge-Kutta公式都是单步显式公式。 例2试用Euler方法、改进的Euler方法及四阶经典R-K方法在不同步长下计算初值问 y(0)=1 在0.2、0.4、0.8、1.0处的近似值,并比较它们的数值结果. 解对上述三种方法,每执行一步所需计算f(x,y)=—y(1+xy)的次数分别为1、2、4。 为了公正起见,上述三种方法的步长之此应为1: 2: 4。 因此,在用Euler方法、改进的Euler 方法及四阶经典R-K方法计算0。 2、0。 4、0。 8、1。 0处的近似值时,它们的步长应分别取为0。 05、0。 1、0。 2,以使三种方法的计算量大致相等。 Euler方法的计算格式为 yn1=yn-0.05[yn(1xnyn)]. 改进的Eluer方法的计算格式为 yp=yn—0.W[yn(1+xnyn)], «yc=yn—0.1K[yp(1+xn”p)], 1,,、 yn*=^(yp+yc). 'J二 四阶经典R-K方法的计算格式为 ynT=yn*&*(k1+2k2+2k3*k4), k〔=—yn(1+xnyn),.02.k2=—(yn+°-Xk1)[1+(xn2k3=-(Yn+°;2*2)[1+M2k4=-(yn0.2k3)[1(xn0.2)(yn0.2k3)] 初始值均为y0=y(0)=1,将计算结果列于表9.2. 表9.2 Euler方法 (步长h=0.05) 改进的Euler方法 (步长h=0.1) 四阶经典R-K方法 (步长h=0.2) 准确解 xn yn yn yn y(xn) 0.2 0.8031866 0.8052632 0.8046363 0.8046311 0.4 0.6271777 0.6325651 0.6314653 0.6314529 0.6 0.4825586 0.4905510 0.4891979 0.4891800 0.8 0.3693036 0.3786397 0.3772249 0.3772045 1.0 0.2827482 0.2923593 0.2910086 0.2909884 从表9.2可以看出,在计算量大致相等的情况下,Euler方法计算的结果只有2位有效数字, 改进的Euler方法计算的结果有3位有效数字,而四阶经典R-K方法计算的结果却有5位有效数字,这与理论分析是一致的。 例1和例2的计算结果说明,在解决实际问题时,选择恰当的算法是非常必要的。 需要指出的是Runge-Kutta方法的基于Taylor展开法,因而要求解具有足够的
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 微分方程 数值 解法 教材