计算方法―-常微分方程的数值解法3.16

2020/1/181第九章常微分方程的数值解法9.1、引言9.2、初值问题的数值解法--单步法9.3、龙格-库塔方法9.4、收敛性与稳定性主要内容2020/1/182主要内容研究的问题数值解法的意义9.1引言2020/1/183现实世界中大多数事物内部联系非常复杂找出其状态和状态变化规律之间的相互联系，也即一个或一些函数与他们的导数之间的关系此种关系的数学表达就为微分方程1.什么是微分方程?其状态随着时间、地点、条件的不同而不同2020/1/184如何利用数值方法求解微分方程（组）的问题。2.数值求解微分方程的意义如何建立数学模型？各类微分方程本身和他们的解所具有的特性已在常微分方程及数学物理方程中得以解释，本章专门讨论2020/1/185常微分方程表示方法:xy-2y=4xy=f(x,y)=y(1)=-3(,)()2y4x2y4xdy2yfxy4dxxy13例方程令：且给出初值就得到一阶常微分方程的初值问题：2020/1/18600(,)(1.1)()(1.2)yfxyaxbyxy本章主要讨论一阶常微方程的初值问题基本知识其中f(x,y)是已知函数，(1.2)是定解条件也称为初值条件。各种数值解法2020/1/187则称f(x,y)对y满足李普希兹条件，L称为Lipschitz常数.常微分方程的理论指出：当f(x,y)定义在区域G=(a≤x≤b,｜y｜＜∞)若存在正的常数L使：1212|(,)(,)|||(1.3)fxyfxyLyy12[,],xabyy使得对任意的及都成立就可保证方程解的存在唯一性(Lipschitz)条件2020/1/188可以证明,如果函数在带形区域R=a≤x≤b,-∞＜y＜∞}内连续，且关于y满足李普希兹(Lipschitz)条件，即存在常数L(它与x,y无关)使2121),(),(yyLyxfyxf对R内任意两个都成立,则方程(7.1)的解在a,b上存在且唯一。21,yy)(xyy2020/1/189寻找解析解的过程称为求解微分方程组。一个或一组具有所要求阶连续导数的解析函数，将它代入微分方程（组），恰使其所有条件都得到满足的解称为解析解(或古典解),称为真解或解。3.什么是微分方程(组)的解析解?在高等数学中，对于常微分方程的求解，给出了一些典型方程求解析解的基本方法，如可分离变量法、常系数齐次线性方程的解法、常系数非齐次线性方程的解法等。2020/1/1810'0022220222(,)()11221(0)11cos1cossin(0)111sinxyfxyyxyxydyxdxycxcyyyxcyyxyxdyxdxxcyyyyx本章重点讨论一阶常微分方程初值问题dydxdydx当x=0时,y=1,可得c=1特解当x=0时,y=1,可得c=-1特解两边积分通解常微分方程解法回顾-2020/1/1811但能求解的常微分方程仍然是有限的，大多数的常微分方程是不可能给出解析解。譬如22yxy这个一阶微分方程就不能用初等函数及其积分来表达它的解。3.什么是微分方程(组)的解析解?再如，方程1)0(yyy的解,虽然有表可查,但对于表上没有给出的值,仍需插值方法来计算xeyxe2020/1/1812在微分方程中,自变量的个数只有一个,称为常微分方程.自变量的个数为两个或两个以上的微分方程叫偏微分方程。222222222y=4xy-(y)+12xy=000ydydytydtdtTTTxyz一阶线性二阶非线性一阶非线性二阶偏微分方程2020/1/18134.什么是微分方程的数值解?虽然求解微分方程有许多解析方法，但解析方法只能够求解一些特殊类型的方程，从实际意义上来讲，我们更关心的是某些特定的自变量在某一个定义范围内的一系列离散点上的近似值。寻找数值解的过程称为数值求解微分方程。把这样一组近似解称为微分方程在该范围内的数值解2020/1/1814在大量的实际方程中出现的函数起码的连续性都无法保证，更何况要求阶的导数求解数值解很多微分方程根本求不到问题的解析解！重要手段。2020/1/1815常微分方程的数值解法常用来求近似解根据提供的算法通过计算机便捷地实现5.常微分方程数值解法的特点数值解法得到的近似解(含误差)是一个离散的函数表.2020/1/1816我们以下的讨论，都在满足上述条件下进行。若f(x,y)在区域G连续，关于y一阶常微分方程的初值问题的解存在且唯一。满足李普希兹条件一阶常微分方程组常表述为：1111(,,,)()(,,,)mmmmyfxyyaxbyfxyy1010()()mmyxyx方程组初值条件2020/1/1817写成向量形式为(1)()00000(,)(),(,,)mTyfxyaxbyxyxxx高阶常微分方程定解问题如二阶定解问题：(,,)()()yfxyyaxbyaya2020/1/1818也就解决了高阶方程的定解问题。这些解法都可以写成向量形式用于一阶常微分方程组的初值问题。2020/1/1819简单的数值方法与基本概念1.简单欧拉法（Euler）2．后退的欧拉法3．梯形法4．改进Euler法9.2、初值问题的数值解法―单步法2020/1/18201.简单的欧拉(Euler)方法考虑模型：00(,)(1.1)()(1.2)yfxyaxbyxy在精度要求不高时使用。通过欧拉方法的讨论弄清常微方程初值问题数值解法的一些基本概念和构造方法的思路.欧拉方法最简单而直观实用方法2020/1/18212.欧拉方法的导出把区间［a,b］分为n个小区间步长为要计算出解函数y(x)在一系列节点()iiyyx)/()(iiixaihhhban,一般取即等距节点处的近似值01naxxxb1(-)iiihxxN等分2020/1/182200(,)(1.1)()(1.2)yfxyaxbyxy对微分方程(1.1)两端从1nnxx到进行积分11(,())nnnnxxxxydxfxyxdx11()()(,())nnxnnxyxyxfxyxdx2020/1/1823右端积分用左矩形数值求积公式2()()()()()2baggxdxbagaba()(,())gxfxyx令1)1(,())((,())nnnnxxnnfxyxnnyyfxyxh得1()()()(,)nnnnnnyxyxhyxyhfxy或用向前差商近似导数1()()()nnnyxyxyxhx0x12020/1/1824)1,...,0(),(1niyxfhyyiiii亦称为欧拉折线法/*Euler’spolygonalarcmethod*/100021111(,)(,)(,)nnnnyyhfxyyyhfxyyyhfxy依上述公式逐次计算可得：故也称Euler为单步法。ny公式右端只含有已知项所以又称为显格式的单步法。每步计算只用到1nyny2020/1/18253.欧拉公式有明显的几何意义依此类推得到一折线000011122(,)()(,)()(,)(,)xyyxxyyxxxxyxy过点的曲线是解在作的切线与直线交于再作切线交于2020/1/1826就是用这条折线近似地代替曲线欧拉方法()yx()yyxxynx1nxnp1np1npx从几何意义上得知，由Euler法所得的折线明显偏离了积分曲线，可见此方法非常粗糙。2020/1/1827利用Euler方法求初值问题解此时的Euler公式为例20,21122xyxy0)0(y的数值解.此问题的精确解是2()/(1).yxxx)211(221nnnnyxhyy00,0,1,2yn221(,)2102fxyyxx分别取步长h=0.2,0.1,0.05,计算结果如下2020/1/1828hxnyny(xn)y(xn)-ynh=0.20.000.400.801.201.602.000.000000.376310.542280.527090.466320.406820.000000.344830.487800.491800.449440.400000.00000-0.03148-0.05448-0.03529-0.01689-0.00682h=0.10.000.400.801.201.602.000.000000.360850.513710.509610.458720.404190.000000.344830.487800.491800.449440.400000.00000-0.01603-0.02590-0.01781-0.00928-0.00419h=0.050.000.400.801.201.602.000.000000.352870.500490.500730.454250.402270.000000.344830.487800.491800.449440.400000.00000-0.00804-0.01268-0.00892-0.00481-0.002272020/1/18294.欧拉法的局部截断误差：定义在假设yi=y(xi)，即第i步计算是精确的前提下，考虑的截断误差Ri=y(xi+1)yi+1称为局部截断误差/*localtruncationerror*/。定义若某算法的局部截断误差为O(hp+1)，则称该算法有p阶精度。如果单步差分公式的局部截断误差为O(hp+1)，则称该公式为p阶方法。这里p为非负整数。显然,阶数越高，方法的精度越高。2020/1/1830Ri的主项/*leadingterm*/欧拉法的局部截断误差：232()()hiyxOh21132()[()()()()][(,)]iiihiiiiiiRyxyyxhyxyxOhyhfxy欧拉法具有1阶精度。Taylor展开式,一元函数的Taylor展开式为:321!3)(!2)()()()()(hxyhxyhxyxyhxyxynnnnnn2020/1/18315.欧拉公式的改进：隐式欧拉法/*implicitEulermethod*/向后差商近似导数hxyxyxy)()()(011x0x1))(,()(1101xyxfhyxy)1,...,0(),(111niyxfhyyiiii由于未知数yi+1同时出现在等式的两边，不能直接得到，故称为隐式/*implicit*/欧拉公式，而前者称为显式/*explicit*/欧拉公式。一般先用显式计算一个初值，再迭代求解。隐式欧拉法的局部截断误差：11)(iiiyxyR)()(322hOxyih即隐式欧拉公式具有1阶精度。见P2852020/1/18326.梯形公式/*trapezoidformula*/—显、隐式两种算法的平均)1,...,0()],(),([2111niyxfyxfhyyiiiiii注：的确有局部截断误差，即梯形公式具有2阶精度，比欧拉方法有了进步。但注意到该公式是隐式公式，计算时不得不用到迭代法，其迭代收敛性与欧拉公式相似。)()(311hOyxyRiii中点欧拉公式/*midpointformula*/中心差商近似导数hxyxyxy2)()()(021x0x2x1))(,(2)()(1102xyxfhxyxy1,...,1),(211niyxfhyyiiii假设，则可以导出即中点公式具有2阶精度。)(),(11iiiixyyxyy)()(311hO