数值分析常微分方程数值解.pdf
《数值分析常微分方程数值解.pdf》由会员分享,可在线阅读,更多相关《数值分析常微分方程数值解.pdf(13页珍藏版)》请在得力文库 - 分享文档赚钱的网站上搜索。
1、许多实际问题的数学模型是微分方程或微分方程的定解问题。如物体运动、电路振荡、化学反映及生物群体的变化等。常微分方程可分为线性、非线性、高阶方程与方程组等类;线性方程包含于非线性类中,高阶方程可化为一阶方程组。若方程组中的所有未知量视作一个向量,则方程组可写成向量形式的单个方程。因此研究一阶微分方程的初值问题 0)(),(yaybxayxfdxdy,(9-1)的数值解法具有典型性。常微分方程的解能用初等函数、特殊函数或它们的级数与积分表达的很少。用解析方法只能求出线性常系数等特殊类型的方程的解。对非线性方程来说,解析方法一般是无能为力的,即使某些解具有解析表达式,这个表达式也可能非常复杂而不便计
2、算。因此研究微分方程的数值解法是非常必要的。只有保证问题(9-1)的解存在唯一的前提下,研究其数值解法或者说寻求其数值解才有意义。由常微分方程的理论知,如果(9-1)中的),(yxf满足条件(1)),(yxf在区域 ),(ybxayxD,上连续;(2)),(yxf在上关于满足 Lipschitz 条件,即存在常数,使得 yyLyxfyxf),(),(则初值问题(9-1)在区间,ba上存在惟一的连续解)(xyy。在下面的讨论中,我们总假定方程满足以上两个条件。所谓数值解法,就是求问题(9-1)的解)(xyy 在若干点 bxxxxaN210 处的近似值),2,1(Nnyn的方法。),2,1(Nny
3、n称为问题(9-1)的数值解,nnxxh1称为由到1nx的步长。今后如无特别说明,我们总假定步长为常量。建立数值解法,首先要将微分方程离散化,一般采用以下几种方法:(1)用差商近似导数 在问题(9-1)中,若用向前差商hxyxynn)()(1代替)(nxy,则得)1,1,0()(,()()(1Nnxyxfhxyxynnnnn)(nxy用其近似值代替,所得结果作为)(1nxy的近似值,记为1ny,则有 1(,)(0,1,1)nnnnyyhf xynN 这样,问题(9-1)的近似解可通过求解下述问题 100(,)(0,1,1)()nnnnyyhf xynNyy x(9-2)得到,按式(9-2)由初
4、值经过步迭代,可逐次算出Nyyy,21。此方程称为差分方程。需要说明的是,用不同的差商近似导数,将得到不同的计算公式。(2)用数值积分法 将问题(9-1)中的微分方程在区间,1nnxx上两边积分,可得)1,1,0()(,()()(11Nndxxyxfxyxynnxxnn(9-3)用1ny,分别代替)(1nxy,)(nxy,若对右端积分采用取左端点的矩形公式,即),()(,(1nnxxyxhfdxxyxfnn 同样可得出显式公式(9-2)。类似地,对右端积分采用其它数值积分方法,又可得到不同的计算公式。(3)用 Taylor 多项式近似。把1()ny x在点处 Taylor 展开,取一次多项式近
5、似,则得 2121()()()()2!()(,()(),2!nnnnnnnnhy xy xhy xyhy xhf xy xyxx 设1h,略去余项,并以代替()ny x,便得 1(,)nnnnyyhf xy 以上三种方法都是将微分方程离散化的常用方法,每一类方法又可导出不同形式的计算公式。其中 Taylor 展开法,不仅可以得到求数值解的公式,而且容易估计截断误差。上面我们给出了求解初值问题(9-1)的一种最简单的数值公式(9-2)。虽然它的精度比较低,实践中很少采用,但它的导出过程能较清楚地说明构造数值解公式的基本思想,且几何意义明确,因此它在理论上仍占有一定的地位。1 简单的数值方法和基本
6、概念 1.1 Euler 法与向后 Euler 法 一、Euler 法 Euler 方法就是用差分方程初值问题 10(,)(0,1,1)()nnnnyyhf xynNyy a(9-4)的解来近似微分方程初值问题(9-1)的解,即由公式(9-4)依次算出()ny x的近似值(1,2,)nyn。从几何上看,微分方程(,)yf x y 在xoy平面上确定了一个向量场:点(,)x y处的方向斜率为(,)f x y。问题(9-1)的解()yy x代表一条过点00(,)xy的曲线,称为积分曲线,且此曲线上每点的切向都与向量场在这点的方向一致。从点000(,)P xy出发,以00(,)f xy为斜率作一直线
7、段,与直线1xx交于点111(,)P x y,显然有1000(,)yyhf xy,再从出发,以11(,)f x y为斜率作直线段推进到2xx上一点222(,)P xy,其余类推,这样得到解曲线的一条近似曲线,它就是折线012P PP。因此 Euler方法又称为 Euler 折线法。二、向后 Euler法 在微分方程离散化时,用向后差商代替导数,即11()()()nnny xy xy xh,则得到如下差分方程 11100(,)(0,1,1)()nnnnyyhf xynNyy x(9-5)用这组公式求问题(9-1)的数值解称为向后 Euler 法。向后 Euler法与 Euler法形式上相似,但实
8、际计算时却复杂得多。Euler 法计算1ny的公式中不含有1ny,这样的公式称为显式公式;向后 Euler 法计算1ny的公式中含有1ny,称为隐式公式。显式公式与隐式公式各有特点。显式公式的优点是使用方便,计算简单,效率高。其缺点是计算精度低,稳定性差;隐式公式正好与它相反,它具有计算精度高,稳定性好等优点,但求解过程很复杂,一般采用迭代法。为了结合各自的优点,通常将显式公式与隐式公式配合使用,由显式公式提供迭代初值,再经隐式公式迭代校正。上面隐式公式中,在求解1ny时,1,nnyx为已知,1ny是方程111(,)nnnnyyhf xy的根。一般说来,这是一个非线性方程,因此我们通过构造简单
9、迭代法来求解。迭代格式为(0)1(1)()111(,)(,)(0,1,2,)nnnnkknnnnyyhf xyyyhf xyk 由于(,)f x y满足 Lipschitz 条件,所以(1)()()11111111(,)(,)kkknnnnnnnnyyh f xyf xyhL yy 由此可知,只要1hL,迭代法就收敛到解1ny。1.2 梯形公式 利用数值积分方法将微分方程离散化时,若用梯形公式计算式(9-3)中右端积分,即 111(,()(,()(,()2nnxnnnnxhf x y x dxf xy xf xy x 并用1,nnyy代替1(),()nny xy x,则得计算公式 111(,)
10、(,)2nnnnnnhyyf xyf xy(9-6)这就是求解初值问题(9-1)的梯形公式。梯形公式也是隐式格式,一般需用迭代法求解,迭代公式为(0)1(1)()111(,)(,)(,)(0,1,)2nnnnkknnnnnnyyhf xyhyyf xyf xyk 由于函数(,)f x y关于满足 Lipschitz 条件,所以(1)()()11111111(,)(,)22kkknnnnnnnnhhLyyf xyf xyyy 其中为 Lipschitz 常数。因此,当012hL时,迭代法就收敛到解1ny。1.3 局部截断误差与方法的精度 为了刻画近似解的准确程度,引入局部截断误差与方法精度的概念
11、。定义 假设在某一步的近似解是准确的,即()nnyy x(这个假设称为局部化假设)。在此前提下,用某公式推算所得1ny,我们称 111()nnnRy xy 为该公式(即该方法)的局部截断误差。定义 9.2 如果某种方法的局部截断误差是 1111()()pnnnRy xyO h 则称该方法是阶方法,或具有阶精度。显然越大,方法的精度越高。1)Euler 法的截断误差 假设问题的解)(xy充分光滑,且前步计算结果是精确的,即)(iixyy,()(,()iiiy xf x y x)(ni 于是 Euler法的截断误差是 11111()()(,)()()()nnnnnnnnnnRy xyy xyhf
12、xyy xy xhy x 23()()2nhyxO h(9-7)这里2()2nhyx称为局部截断误差主项。显然21()nRO h。2)向后 Euler 法的截断误差。计算公式是),(111nnnnyxhfyy 将(,)f x y对用微分中值定理,有 1111111(,)(,()(,)()nnnnynnnf xyf xy xfxyy x(在1ny与1()ny x之间)将11(,()nnf xy x在处 Taylor 展开 2111(,()()()()()nnnnnf xy xy xy xhyxO h 于是 231111()()()(,)()()nnnnynnnyy xhy xh y xhfxyy
13、 xO h 将方程的解作 Taylor 展开 231()()()()()2nnnnhy xy xhy xyxO h 因此 2311111()()(,)()()22nnnynnnhhy xyyxfxyy xO h 故 2311112311()()()1(,)21(,)()()2nnnnynynnhRy xyyxO hhfxhhfxyxO h 23()()2nhy xO h(9-8)3)梯形法的计算公式是 111(,)(,)2nnnnnnhyyf xyf xy 将(,)f x y对用微分中值定理,有 1111111(,)(,()(,)()nnnnynnnf xyf xy xfxyy x(在1ny与
14、1()ny x之间)将11(,()nnf xy x在处 Taylor 展开 23111(,()()()()()()2nnnnnnhf xy xy xy xhyxyxO h 1Eny 1Bny 1Tny 图 9-3 于是 2341111()()()()(,)()()242nnnnnynnnhhhyy xhy xyxyxfxyy xO h 将方程的解作 Taylor 展开 2341()()()()()()26nnnnnhhy xy xhy xyxyxO h 因此 3411111()()(,)()()122nnnynnnhhy xyyxfxyy xO h 故 341111341341()()()1(
15、,)121(,)()()12()()12nnnnynynnnhRy xyyxO hhfxhhfxyxO hhyxO h(9-9)所以梯形法是二阶方法。1.4 改进的 Euler 法 我们看到,虽然梯形方法提高了精度,但其算法复杂,在应用迭代公式进行实际计算时,每迭代一次,都要重新计算函数(,)f x y的值,而迭代又要反复进行若干次,计算量很大。为了控制计算量,通常只迭代一两次就转入下一步的计算,这就简化了算法。具体地说,我们先用欧拉公式求得一个初步的近似值1ny,称之为预测值,预测值1ny的精度可能很差,再用梯形公式将它校正一次得1ny,称为校正值。这样的预测校正系统通常称为改进的欧拉公式。
16、即 预测:1(,)nnnnyyhf xy 校正:111(,)(,)2nnnnnnhyyf xyf xy 为了编制程序上机,将上式改写成 1(,)(,)1()2pnnnqnnpnpqyyhf xyyyhf xh yyyy(9-10)算法(1)输入,(,)a b f x y,整数,初值(2)置0,0,bahnxa yyN 输出(,)x y(3)计算(,)(,)pqpyyhf x yyyhf xh y 1(),2pqyyy xhx 输出(,)x y(4)若1nN,置1nn,转 3;否则停机。改进欧拉法的截断误差。计算公式是),(,(),(211nnnnnnnnyxhfyxfyxfhyy 在)(,(n
17、nxyx处作 Taylor 展开式,注意到)(nnxyy,有)()(2)()()()(,()(,()(,()(,(2)(,(2)(3221hOxyhxyhxyhOxyxfxyxhfxyxhfxyxfhxyxfhxyynnnnnnnynnxnnnnnn 于是 3111()()nnnRy xyO h(9-11)所以改进 Euler法是二阶方法。2 龙格-库塔(Runge-Kutta)法 2.1 Runge-Kutta法的基本思想及一般形式 设初值问题(9-1)的解1(),yy xC a b,按微分中值定理,必存在1,nnxx,使 1*()()()(,()nnnny xy xhyyhfyyhK(9-
18、12)其中叫做()y x在1,nnxx上的平均斜率。只要对平均斜率提供一种算法,公式(9-12)就给出一种数值解公式。例如,用1(,)nnKf xy代替,就得到前进 Euler 公式;用211(,)nnKf xy替代,就得到后退Euler 公式;如果用12,K K的算术平均值替代,则可得到 2阶精度的梯形公式,如图 9-3。可以设想,如果在1,nnxx上能多预测几个点的斜率值,用它们的加权平均值替代,就有望得到具有较高精度的数值解公式,这就是 Runge-Kutta 法的基本思想。Runge-Kutta 公式的一般形式是 11111(,)(2,)(,)rnniiinniininijjjyyhc
19、 KKf xyirKf xh yhK(9-13)其中是()yy x在 (01)niixh点的斜率预测值。,iiijc 均为常数。选取这些常数的原则是使公式(9-13)具有尽可能高的精度。公式(9-13)叫做级显式 Runge-Kutta 公式,简称 RK方法。2.2 二阶 RK法的推导 2r 的 RK方法的近似公式为 11122122211()(,)(,)nnnnnnyyh c Kc KKf xyKf xh yhK(9-14)这里12221,c c 均为待定常数,我们希望适当选取这些系数,使公式阶数尽量高。先展开,按照二元函数 Taylor级数:01(,)(,)!knkf as btstf a
20、 bkxy,得 22222122222221211(,)(,)(,)(,)(,)2!2(,)(,)(,)(,)nnxnnynnnnxxnnxynnnnyynnnnKf xyhfxyhfxyf xyh fxyh fxyf xyh fxyfxy(9-15)为了叙述方便,把(,)nnf xy及其偏导数中的,nnxy省略不写,将式(9-15)代入(9-14)的第一式得 2112222132222222121()()(2)2nnxyxxxyyyyyh ccfh cff fchfffff 再展开1()ny x,注意到公式 2(,)2()xyxxxyyyyxyyf x yyff fyffffffff f 得
21、 231232()()()()()2!3!()(2)26nnnnnnxyxxxyyyyxyhhy xy xhy xyxyxhhyhfff fffffff ff f 于是,局部截断误差为 111212222213222222221221()11(1)2211111()62362nnnxyxxxyyyyxyRy xyhccfhcfcf fhcfcffcfffff f (9-16)要使式(9-16)的局部截断误差为3()O h,则应要求,xyfff f的系数为零,于是有 122222110.50.5cccc(9-17)方程有 4个未知数,3个方程,所以有无穷多组解,它的每组解代入式(9-14)得到的
22、近似公式,局部截断误差均为3()O h,故这些方法统称为二阶方法。例如,取120.5cc,得2211,近似公式为 1121210.5()(,)(,)nnnnnnyyh KKKf xyKf xh yhK(9-18)这就是改进 Euler 公式。如果取120,1cc,有22112,近似公式为 12121(,)(0.5,0.5)nnnnnnyyhKKf xyKf xh yhK(9-19)这也是常用的二阶公式,称为中点公式。注:对于一般函数(,)f x y,由于()0yxyfff f,所以不论参数如何选取,只能有 31nRO h。这说明(9-14)至多是二阶的方法。类似地,对3r 和4r 的情形,通过
23、更复杂的计算,可以导出三阶和四阶 RK公式,其中常用的三阶和四阶RK公式为 1123121312(4)6(,)(,)22(,2)nnnnnnnnhyyKKKKf xyhhKf xyKKf xh yhKhK(9-20)和 112341213243(22)6(,)(,)22(,)22(,)nnnnnnnnnnhyyKKKKKf xyhhKf xyKhhKf xyKKf xh yhK(9-21)式(9-21)称为经典的四阶 RK公式,通常说四阶 RK方法就是指用式(9-21)求解。算法(1)输入0,(,),a b f x yN y(2)置0,1,bahnxaN(3)计算 01002001300240
24、0301234(,)(,)22(,)22(,)(22)6xxhKf xyhhKf xyKhhKf xyKKf xh yhKhyyKKKK 输出(,)x y(4)若nN,则置001,nn xxyy 转 3;否则停机。例 4 用 Euler法(0.025h),改进 Euler 法(0.05h)和经典 R-K方法(0.1h)计算初值问题 1 (00.5)(0)0yyxy .其结果如下表 9-1 表 9-1 Euler 方法(0.025h)改进 Euler法(0.05h)经典 R-K法(0.1h)()ny x 0 0 0 0 0 注:1)通过比较 RK四阶经典公式和前面用 Euler法、改进 Eule
25、r 法(它是一种二阶 RK方法)的计算结果,我们发现,四阶 RK方法的精度高得多。就计算量来说,虽然 Euler法、改进 Euler 法每步只需计算一个或二个函数值,而四阶 RK方法每步需计算四个函数值,但由于放大了步长,计算量几乎相同。可见,选择有效的算法是非常重要的。2)RK方法的导出基于 Taylor 展开,故它要求问题的解具有较高的光滑度。当解充分光滑时,四阶 RK方法确实优于改进 Euler 法。如果解的光滑性较差,则用四阶 RK方法求得的数值解,其精度可能反而比改进 Euler 法的差。因此,在实际计算时,需根据问题的具体情况选择合适的算法。3)当1,2,3,4r 时,RK公式的最
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 微分方程
限制150内