常微分方程数值解法

第八章常微分方程数值解法教学目的1.掌握解常微分方程的单步法：Euler方法、Taylor方法和Runge-Kutta方法；2.掌握解常微分方程的多步法：Adams步法、Simpson方法和Milne方法等；3.了解单步法的收敛性、相容性与稳定性；多步法的稳定性。教学重点及难点重点是解常微分方程的单步法：Euler方法、Taylor方法和Runge-Kutta方法和解常微分方程的多步法：Adams步法、Simpson方法和Milne方法等；难点是理解单步法的收敛性、相容性与稳定性及多步法的稳定性。教学时数20学时教学过程§1基本概念1．1常微分方程初值问题的一般提法常微分方程初值问题的一般提法是求函数bxaxy),(，满足)2.1()()1.1(),,(aybxayxfdxdy其中),(yxf是已知函数，是已知值。假设),(yxf在区域},),{(ybxayxD上满足条件：（1）),(yxf在D上连续；（2）),(yxf在D上关于变量y满足Lipschitz条件：2121),(),(yyLyxfyxf，21,,yybxa（1.3）其中常数L称为Lipschitz常数。我们简称条件（1）、（2）的基本条件。由常微分方程的基本理论，我们有：定理1当),(yxf在D上满足基本条件时，一阶常微分方程初值问题（1.1）、（1.2）对任意给定存在唯一解)(xy在],[ba上连续可微。定义1方程（1.1）、（1.2）的解)(xy称为适定的，若存在常数0和0K，对任意满足条件及)(x的和)(x，常微分方程初值问题aazbxaxzxfdxdz)(),(),(（1.4）存在唯一解)(xz，且}.{)()(Kxzxy适定问题的解)(xy连续依赖于（1.1）右端的),(yxf和初值。由常微分方程的基本理论，还有：定理2当),(yxf在D上满足基本条件时，微分方程（1.1）、（1.2）的解)(xy是适定的。我们在本章中假设),(yxf在D上满足基本条件，从而（1.1）、（1.2）的解)(xy存在且适定。一般的一阶常微分方程组初值问题是求解niaybxaniyyxfydxdiinii,,2,1,)(,,,2,1),,,,(1（1.5）（1．5）的向量形式是)(),,(aybxayxFydxd（1.5）其中.),,,(,)),(,),,((),(,))(,),(()(2111TnTnTnyxfyxfyxFxyxyxy记},,2,1,,),,,{(1niybxayyxDin。类似于定理1和定理2，我们有：定理3若映射),(yxF满足条件(1)),(yxF在D上是从1nR到nR上的连续映射;(2)),(yxF在D上关于y满足Lipschits条件;212121,,),(),(yybayyLyxFyxF任意。则常微分方程组初值问题（1.5）存在的唯一的连续可微解),(xy而且解)(xy是适定的。高阶常微分方程初值问题一般为1,,1,0,)(),,,,,(111niaaydxddxaydxddxdyyxfydxdiiinnnn（1.6）其中),,,,(uuyxf是给定多元函数，naa,1为给定值。引进新的变量函数nkbxaxydxdxykkk,,2,1,),()(11（1.7.）则初值问题（1.6）化成了一阶常微分方程组初值问题niayyyxfdxdyyydxdbxayydxdiinnnn,2,1,)(),,,()8.1(1121通过求解（1.8）得到(1.6)的解)()(1xyxy。1.2初值问题数值解基本概念初值问题的数值解法，是通过微分方程离散化而给出解在某些节点上的近似值。在ba,上引入节点),,1(,:1100nkxxhbxxxaxkkknnkk称为步长。在多数情况下，采用等步长，即),1,0(,nkkhaxnabhk。记（1.1）,(1.2)的为准确解为)(xy,记)(kxy的近似值为ky,记),(kkyxf为kf.。求值问题数值解的方法是步进法，即在计算出kiyi,后计算1ky。数值的方法有单步与单步法之分。单步法在计算1ky时只利用ky而多步法在计算1ky时不仅要利用ky还要利用前面已算出的若干个1,,2,1,ljyjk。我们称要用到11,,,lkkkyyy的多步法为l步方法。单步法可以看作多步法，但两者有很大差别。l步方法只能用于lkyk,的计算，110,,,lyyy要用其它的方法计算；而且在稳定性上单性法比1l的多步法容易分析；此外单步法容易改变步长。单步法和多步法又都有显式方法和稳式方法之分。单步显式法的计算公式可写成),,(1hyxhyykkkk（1.9）隐式单步法的计算公式可写成),,,(11hyyxhyykkkkk（1.10）在(1.10)中右端项显含1ky。从而（1.10）是1ky的方程式，要通过解方程求出1ky。显式多步法计算公式为),,,(11hyyxhyykkkkk（1.11）而隐式多步法计算公式为),,,,,)(111hyyyxhyylkkkkkk(1.12)右端项含1ky。多步法中一类常用方法是线性多步法1,11101lkfihyayliikiliikik（1.13）其中101110,,,,.,,lla是独立于k和f的常数。01时（1.13）是显式的，01时是隐式的。数值解法及到方法构造、误差分析、稳定性分析等内容。一些概念和定义在后面的论述中逐步引入。§2Euler的方法Euler方法是常微分方程初值问题数值方法中最简单的方法。Euler方法精度低，较少有直接使用。但我们通过Euler方法介绍离散化途径、数值解法中的基本概念、术语和加速方法等。2.1显式Euler方法设节点为bxxxan10。初值问题（1.1）、（1.2）的显式Euler方法为1,,1,0,10nkfhyyaykkkk（2.1）其中),(,1kkkkkkyxffxxh。显式Euler方法可以用多种途径导出。导出方法1：Taylor展开法。将)(1kxy在kxx点进行Taylor展开,得121,,!2)())(,()()(kkkkkkkkkkxxhyxyxfhxyxy（2.2）忽略2kh这一阶项,分别用kkkkkyxffyy,,,1近似)(kxy，)(1kxy和))(,(kkxyxf，得kkkkfhyy1。结合初值条件)0(y即得（2.1）。导出方法2：向前差分近似微分法。用向前差分kkkhxyxy)()(1近似微分)('kxy，得))(,()()(1kkkkkxyxfhxyxy（2.3）将近似号改作等号，用kkkfyy,,1近似)(kxy，)(1kxy))(,(kkxyxf，并结合初值条件即得（2.1）。方法3：左矩数值积分法。将（1.1）两边从kx到1kx积分得dxxyxfxyxykkxxkk1))(,()()(1（2.4）用ky，1ky近似)(kxy、)(1kxy，数值积分采用左矩公式得),(1kkkkkyxfhyy，从而亦得（2.1）Euler方法有几何意义，如图8-1，式（1.1）,(1.2)的解曲线)(xy过点),(000yxP，且具斜率0f。从0P出发以0f为斜率作直线段，交1xx于),(111yxP，显然0001fhyy。式（1.1）过),(111yxP的解曲线具有斜率1f从1P出发以1f为斜率作直线要交2xx于),(222yxP，余类推。这样我们得到了一条折线nPPP10，它在点kP的右侧具有斜率kf，与（1.1）过kP的解曲线相切。我们取折线nPPP10，作为（1.1）、（1.2）解曲线)(xyy的近似曲线，所以Euler方法又称折线法。2.2隐式Euler方法和梯形方法若将)(kxy在1kx展开11()()(kxkkxfhxyxy、,)(!!21))(21kkkhyxy1kkkxx忽略2h项，用1,kkyy和),(111kkkyxff分别近似)(kxy，)(1kxy及))(,(11kkxyxf，可以得另一计算公式1,,1,0),,(111nkyxfhyykkkkk（2.5）（2．5）式称为隐式Enler方法。隐式Euler方法也可以利用向后差分近似微分或用右矩数值求积公式来建立。读者可自行推导。隐式Euler方法(2.5)给出了1ky要满足的方程，要通过解方程才能得到1ky。在显式和稳式Euler方法中，忽略的项都是2h项，为了得到更高精确度的方法，我们可将1211,)(21))(,()()(kkkkkkkkkkxxhyxyxfhxyxy12111,)(21))(,()()(kkkkkkkkkkxxhyxyxfhxyxy取平均得)]()([4))](,())(,([2)()(2111kkkkkkkkkkyyhxyxfxyxfhxyxy当)(xy三次连续可微时,)()()(kkkhOyy。忽略)(3khO项，用1,kkyy分别近似)(),(1kkxyxy，得)],(),([2111kkkkkkkyxfyxfhyy（2.6）（2.6）称为梯形方法。取这个名称的原因是利用梯形求积公式1))(,([2)(,(kkxxkkkxyxfhdxxyxfxxyxfDhxyxfxkkk))(,(12))](,(2311其中xD表示关于x的全微分，忽略数值求积余项也可建立（2.6）。梯形方法也是隐式方法，要通过解（2.6）来得到1ky。与（1.10）式中单步法公式相对应，显式Euler方法取),(),,(kkkkkyxfhyx隐式Euler方法取),(),,(111kkkkkyxfhyx梯形方法取21),(21),,(1kkkkkyxfhyx,(1kxf)1ky当),(yxf在D上满足基本条件，),(yxf关于y的Lipschitz常数为L时，只要)5.2(,1Lhk确定了唯一的1ky；同样，只要,2kLh（2.6）确定了唯一的1ky。以（2.6）为例，当,2kLh以y为变量的函数)],(),([21yxfyxfhykkkkk在y上关于y满足Lipschitz条件，且Lipschitz常数为,12khL从而由第七章压缩不动点定理得方程)],(),([21yxfyxfhyykkkkk有唯一点不动,1yk而且从任意)0(1ky出发，迭代)],,(),([2)(11)1(1ikkkkkkikyxfyxfhyy1,,1,0ni（2.7）都收敛到1ky在实际计算中总希望有较好的)0(1ky，用较少的迭代步，取得有足够精度的1ky。2．3预估–校正Euler方法在实际计算中，),()(11ikkyxf的计算量比较大，往往取)1()(1mymk作为1ky来用。我们称)(1mky为)0(1ky的m次迭代改进。最常用的方法之一是先用显式Euler方法所得的1ky为)0(1ky，再用梯形方法改进一次1,,1,0)],,(),([2),(1111nkyxfyxfhyyyxfhyykkkkkkkkkkkx（2.8）方法（2.8）称为预估-校正Euler方法，或改进Euler方