第六章 常微分方程数值解法简介-2011-04

科学计算与数学建模中南大学数学科学与计算技术学院第7章常微分方程数值解法简介简单的数值方法与基本概念2数值线性多步法3非线性高阶单步法——Runge-Kutte法45一阶方程组和高阶方程的初值问题6常微分方程边值问题的数值解法1实际问题的微分方程模型第七章常微分方程数值解法简介微分方程在科学和工程技术中有很广泛的应用。许多实际问题的数学模型都可以用微分方程来描述，归结为常微分方程的定解问题。很多偏微分方程问题，也可以化为常微分方程问题来近似求解，但是求出所需的解绝非易事。实际上，除了极特殊情形外，人们不可能求出微分方程的解析解，只能用各种近似方法得到满足一定精度的近似解。在常微分方程中已经熟悉了级数解法和Picard逐步逼近法，这些方法可以给出解的近似表达式，称为近似解析方法。另一类方法只给出解在一些离散点上的值，称为数值方法。数值方法应用范围更广，特别适合用计算机计算，本章主要介绍常用的常微分方程数值解法。函数是事物的内部联系在数量方面的反映,如何寻找变量之间的函数关系,在实际应用中具有重要意义.在许多实际问题中，往往不能直接找出变量之间的函数关系,但是有时却容易找出变量的改变量之间的关系，从而建立描述问题的微分方程模型.§1实际问题的微分方程模型例7.1将初始温度的一碗汤放置于环境温度保持在的桌上，10分钟后测得汤的温度为1000C。如果汤的温度低于550C才可以喝，问再过20分钟后这碗汤能喝了吗？00150uCau024C设物体在时刻的温度为，从，温度从，注意到热量总是从温度高的物体向温度低的物体传导，因而，所以温度差恒正，又因物体将随时间而逐渐冷却，则温度的改变量为两边除以，并令得温度变化速度为解：为了解决这一问题，需要了解有关热力学的一些基本规律.热量总是从温度高的物体向温度低的物体传导的；在一定的温度范围内，一个物体的温度变化速度与这个物体的温度和其所在介质温度的差值成正比。t()uutttt()()ututt0auuauu()()(())auuttutkuttut()adukuudtt0t其中是比例常数.从而得出描述物体冷却过程的微分方程模型为容易求出这个一阶微分方程初值问题的解为根据问题所给的条件知,当时，，代入，得0()(0)adukuudtuu0k(7.1.1)0()ktaauuuue(7.1.2)10t01100uuC1010()kaauuuue00150uC024auC0111lnln1.660.0511010aauukuu从而得到这碗汤的温度随时间变化的函数关系为于是，将代入计算得到再过20min汤的温度，这说明再过20min后这碗汤能喝了.不过，并不是所有的微分方程模型都可求出解析解。例如，看似简单的微分方程，自德国数学家WilhelmvonLeibniz提出100多年后才被法国数学家JosephLiouville证明它没有解析解，只能借助于数值的方法求数值解.0.0510.0510()()24126ttaauutuuuee(7.1.3)30t0251uC22dyxydx例7.2某地区发现一种有免疫性的传染病，为了控制疫情扩散对该地人群进行隔离处理.为了分析受感染人数的变化规律，需要建立描述传染病传播过程的数学模型.解设该地区的总人数为常数，任意时刻病人、健康人和病人治愈后移出感染系统的移出者的比例分别为，病人的日接触率，日治愈率，则容易得出从时刻，病人和健康人的改变量为每个方程两边除以，并令，化简后得(),(),()itstrtttt[()()]()()()[()()]()()NittitNstittNittNsttstNstitt00(0),(0)disiidtdssidtiiss(7.1.4)t0t其中（对任意的t）.式（7.1.4）就是描述病人和健康人的比例和随时间变化的微分方程模型，这是一个微分方程组的初值问题.但是，这一初值问和的解析式来分析和解决问题。()()()1stitrt()it()st()it()st在数学建模课程中学到的大量数学模型都是用微分方程形式给出的，各类微分方程本身和它们的解所具有的特性在常微分方程及数学物理方程中有所解释.虽然，求解微分方程有许多解析方法，但解析方法只能够求解一些特殊类型的方程，在实际应用中人们更关心的是某些特定的自变量在某一个定义范围内的一系列离散点上的近似值.这样一组近似解称为微分方程在该范围内的数值解，寻找微分方程数值解的过程称为微分方程的数值解法。§2简单的数值方法与基本概念设在区域上连续，求满足其中是已知常数,这就是一阶常微分方程的初值问题.为使问题（7.2.1）的解存在、唯一且连续依赖初值，即初值问题（7.2.1）适定，还必须对右端项加以适当限制，通常要求关于是已知函数,且满足Lipschitz条件，即存在常数L，使7.2.1常微分方程初值问题(,)fxy:,Gaxby()yyx0(,),()dyfxyaxbdxyay(7.2.1)0y0y(,)fxy(,)fxyy1212(,)(,)fxyfxyLyy(7.2.2)对所有及本章总假定满足条件（7.2.2）。1.Euler方法的导出与几何意义最简单的数值解法是Euler将区间作N称为步长，点列称为节点，。由已知初值，可算出在的导数。[,]xab12,(,)yy7.2.2Euler法及改进的Euler法,(0)abTba/hTNixaih(0,1,...,)iN0xa00()yxy()yx0xx'00000()(,())(,)yxfxyxfxy其中，并略去二阶小量，得下面用3种方法导出Euler法.本章用表示函数在点表示的近似值。就是的近似值。利用可算出，如此下去可算出在所有节点上的值，它的一般递推公式为()iyx()yxixiy()iyx23''''''100000000()()()()()()26(,)hhyxyxhyxhyxyxyyhfxyR(7.2.3)01(,)xx0R1000(,)yyhfxy1）幂级数展开法利用Taylor展式1()yx1y1y2y()yx1(,)nnnnyyhfxy0,1,...,1nN(7.2.4)这就是Euler法。实际上，初值问题（7.2.1）的解是xy平面上过点的一条积分曲线。按Euler法，过初始点作经过此点的积分曲线的切线（斜率为（按式（7.2.4）计算）作为的近似.然后，过（斜率为），沿切线取点（按式（7.2.4）计算）作为的近似.如此下去，即得一条以为顶点的折线，这就是用Euler法得到的近似积分曲线，如图7-1所示.从几何图形上看，越小，此折线逼近积分曲线越好，因此也称Euler法为Euler折线法00(,)xy00(,)xy00(,)fxy11(,)xy1y11(,())xyx11(,)xy11(,)fxy22(,)xy22(,())xyx(,)nnxy2yhEuler法有明显的几何意义：2)数值微分法利用向前差商近似导数从而得出Euler法的一般递推公式1()()()nnnyxyxyxh1()()()(,)nnnnnnyxyxhyxyhfxy1(,)nnnnyyhfxy0,1,...,1nN3)数值积分法将初值问题（7.2.1）写成等价的积分形式：取，得用左矩形公式作为右端积分的近似，并用替代从而也可得出Euler法的一般递推公式为00()()(,())xxyxyxftytdt1xx1010()()(,())xxyxyxftytdt1000(,)yyhfxy1(,)nnnnyyhfxy0,1,...,1nN1y1()yx(7.2.5)2.改进的Euler由Euler方法的数值积分导出法可知只要给出式（7.2.5）右端定积分的一种近似计算方法，就可得出初值问题（7.2.1）的一种如果用右矩形公式作为式（7.2.51011(,)yyhfxy从而也可得出一般递推公式为111(,)nnnnyyhfxy0,1,...,1nN称式（7.2.6）为后退Euler法(7.2.6)显然，改进的Euler法比Euler法精度更高。后退Euler法和改进的Euler法，由于未知数同时出现在等式的两边，故称为隐式算法；如果未知数由已知量直接计算（即不出现在等式右端），则称为显式算法。对于隐式算法，每步计算需要解关于的方程，而这样的方程往往是非线性的，通常将初值取为，用迭代法求解，一般只需迭代几步即可收敛。一般先用显式公式计算一个初值，再用隐式公式迭代求解。如果用梯形公式近似替代式（7.2.5100011[(,)(,)]2hyyfxyfxy111[(,)(,)]2nnnnnnhyyfxyfxy0,1,...,1nN称式（7.2.7）为改进的Euler法(7.2.7)1ny1ny1ny[0]1nnyy如果先用显式Euler公式作预测，算出再将代入隐式梯形公式的右边作校正，得到1ny111[(,)(,)]2nnnnnnhyyfxyfxy1(,)nnnnyyhfxy从而可得11(,),(,)2nnnnnnnnhyyfxyfxyhfxy0,1,...,1nN这种方法称为预估－校正法。可以看到它是显示格式，迭代求解过程比隐式公式的简单。后面将看到，它的稳定性高于显式Euler法。如果在区间11[,]nnxx上对初值问题（7.2.1）的方程两边积分，则有1111()()(,())nnxnnxyxyxfxyxdx并用中矩形求积公式近似替代右端的定积分，则得出一般递推公式为112(,)nnnnyyhfxy0,1,...,1nN称式（7.2.8）为Euler中点公式。(7.2.8)和，这样的方法称为双步法（或1()nyx1nyny0y1nyny1ny1ny如果计算的近似值时只用到前一节点的值，则从初值样的方法称为单步法；而Euler中点公式计算到前两个节点的值多步法.多步法附加初值才能逐一计算出以后各节点的值。出发可逐一计算出以后各节点的值，这时需要用二步法）；计算时需要用到前面多个节点值的方法称为7.2.3截断误差与算法精度的阶从Euler法的几何意义得知，由Euler法所得的折线明显偏离了积分曲线，可见此方法非常粗糙（即误差太大）。现在分析一下求解初值问题（7.2.1）的数值方法误差的来源。为使问题简化，不考虑因计算机字长限制引起的摄入误差。在假设第步计算是精确的前提下），第的截断误差称为局部截断误差。若某算法的局部截断误差为的同阶无穷小，则称该算法有阶精度。若局部截断误差n1n()nnyyx1()nyx11()nnnRyxy1()pOh1php1211()(,)()PpnnnnnRyxyxyhOh则称为误差主项，为误差主项系数。1(,)Pnnxyh(,)nnxy1.Euler法的局部截断误差由Euler法的一般递推公式和的Taylor展式，得所以，Euler法的局部截断误差为，即Euler法为1阶精度算法，其误差主项为1()nyx11()()[(,)]nnnnnnnRyxyyxhyhfxy2()Oh2()2nhyx23()()[()()][()]2!3!nnnnnnyxyxyxyxhhhyhyx2233()()()()2!3!2nnnyxyxhhhyxOh2.后退Euler法的局部截断误差同理，由后退Euler法的一般递推公式，的Taylor展式，得所以，后退Euler法的局部截断误差也是，即后退Euler法也是11