第九章+常微分方程初值问题数值解法

1第九章常微分方程初值问题的数值解法在自然科学的工程技术的许多领域中，常会遇到常微分方程初值问题，但这种问题大多数情况下不存在初等形式的解析解，只能用近似方法来求解。近似解法主要有两类：一类叫做近似解析方法，它能给出解得近似表达式，例如熟知的级数解法和逐次逼近法等；另一类近似解法称为数值解法，它可以给出解在一些离散点上的近似值。数值方法便于电子计算机求解微分方程。本章讨论初值问题常用的数值解法，并介绍有关的基本理论。假设给定一阶常微分方程的初值问题：(,),()dyfxyaxbdxyaa（9.1）（9.2）其中f为x,y的已知函数，α为给定的初始值。由微分方程的理论可知，如果f(x,y)在区域axb，y内连续，且关于y满足李普希兹(Lipschitz)条件，即存在常数L，使|(,)(,)||fxyfxyLyy对所有的axb及任何y，y均成立，则初值问题（9.1），（9.2）有连续可微的解y(x)存在且唯一。所谓初值问题的数值解，则是问题的解y(x)在一系列点01212...NNaxxxxxb处的值()nyx的近似值ny(n=0,1,…,N)。这里相邻两个节点之间的距离1nnnhxx通常称为步长，通常将步长nh取为常数h。2§9.1欧拉方法与改进的欧拉方法欧拉(Euler)方法是最简单的数值方法，由于它的精确程度较差，已不常用于实际计算。但构造这个方法的基本原理，对于构造一般的数值方法具有普遍意义，因此首先对它进行讨论。9.1-1欧拉方法的构造将方程（9.1）中点nx处的导数()nyx用差商近似地表示为1()()()nnnyxyxyxh即在该点有近似等式1()()(,())nnnnyxyxfxyxh用近似值ny代替()nyx，由上式可以导出其近似值满足的差分方程：1(,()),0,1,...,1nnnnyyhfxyxnN（9.3）（9.3）式称为欧拉方法的计算公式或称为欧拉公式。当初始值0y给定时，利用欧拉公式就可以逐次计算出初值问题的数值解01,,...,Nyyy。9.1-2后退的欧拉公式如果在点1nx处用向后差商近似(9.1)式中的导数1()()()nnnyxyxyxh即有111()()(,())nnnnyxyxhfxyx再用近似值iy代替()iyx(,1)inn，则可导出近似值满足的差分方程3111(,)nnnnyyhfxy（9.6）（9.6）式称为后退的欧拉公式。这是一个两端含有未知数1ny的方程，这样的方法称为隐式方法。相应地，称右端不隐含1ny的欧拉公式（9.3）是一个显式方法。隐式方法求1ny时需解方程，当(,)fxy关于y为非线性时，就需要解非线性方程，通常用迭代法求解。可以证明，隐式的后退欧拉公式（9.6）的局部截断误差为23,()()2nhnhRyxOh（9.7）事实上，设()nnyyx时，由（9.6）式有11111(,)()(,)nnnnnnnyyhfxyyxhfxy（9.8）由中值定理得1111111(,)(,())(,)[()]nnnnynnnfxyfxyxfxyyx（9.9）将(9.9)式代入（9.8）式，并利用微分方程（9.1），得11111()()(,)[()]nnnynnnyyxhyxhfxyyx（9.10）由微分方程解的泰勒展式（9.4）与（9.10）式相减，得11211311()[()()]()(,)2[()]()nnnnnynnnyxyhhyxyxyxhfxyyxOh4以31()()()()nnnyxyxhyxOh代入上式，得1123111()()(,)[()]()2nnnynnnyxyhyxhfxyxyOh移项整理得11231()1[()()]1(,)2nnnynyxyhyxOhhfx注意当h充分小，使1yhf时，有1111(,)...1(,)ynynhfxhfx即得2311()()()2nnnhyxyyxOh9.1-3改进的欧拉方法用1ny与1ny分别表示用欧拉公式与后退的欧拉公式求得的数值解，由式（9.5）与（9.7）可知，欧拉方法与后退的欧拉方法的局部截断误差的主部只是符号相反，于是显然有31111()()()2nnnyxyyOh由此可以构造如下的计算公式111[(,)(,)]2nnnnnnhyyfxyfxy(0,1,2,...1)nN（9.11）（9.11）式称为梯形方法，其局部截断误差为3()Oh，因此它是一个二阶方法，梯形方法也是隐式的，在实际应用中时常与显式的欧拉公式联合使用，构成如下的计算格式：5(0)1(1)()111(,)[(,)(,)]2nnnnkknnnnnnyyhfxyhyyfxyfxy(0,1,2,...)k（9.12）即先用欧拉方法算出初始近似值(0)1ny，然后用（9.12）的第二式进行迭代，反复改进这个近似值，直到(1)()11kknnyy（为所允许的误差）为止。而把(1)1kny取作1()nyx的近似值1ny，类似地计算2ny，3ny，…显然，如果上述迭代序列(0)1ny，(1)1ny，…收敛，其极限便满足方程111[(,)(,)]2nnnnnnhyyfxyfxy即序列的极限是梯形方法所得到的解1ny，容易证明，只要h取得充分小，上述迭代过程必定收敛，事实上，将（9.11）式与（9.12）第二式相减，得(1)()111111[(,)(,)2kknnnnnnhyyfxyfxy由于f满足李普希兹条件，故有(1)()11112kknnnnhLyyyy因此只要h充分小，使12hL，迭代过程就是收敛的，但计算时需要迭代多少次，一般无法估计，故在实用时，在h取得较小的条件下，常常让计算迭代一次就结束，将其一次迭代值取作1ny。这时计算公式为(0)1(0)111(,)[(,)(,)]2nnnnnnnnnnyyhfxyhyyfxyfxy(0,1,2,...)n（9.13）6或直接写成110[(,)(,(,))]2nnnnnnnnhyyfxyfxyhxyya(9.14)(9.13)式或(9.14)式被称为改进的欧拉公式，将(9.14)式展开后与初值问题的解的泰勒展开式比较，可知其局部截断误差仍是3()Oh，即改进的欧拉方法是二阶方法，且是显式方法。通常把(9.13)式中第一式得出的初始近似值(0)1ny称为预测值，第二式是对预测值进行一次校正，因此称这样构造的方法为预测-校正方法。［例1］用欧拉方法与改进的欧拉方法求初值问题223(0)1dyxdxyy在区间[0,1]上取0.1h的数值解。［解］欧拉方法的计算公式为1022()1,0.13nnnnxyyhyhy　　改进的欧拉方法其计算公式为(0)12112(0)2102()322[]233()1,0.1nnnnnnnnnnxyyhyxxhyyyyyh本题的精确解为32()1yxx，可用来检验数值解的精确度，列出计算结果。（请完成）。7§9.2龙格库塔方法龙格-库塔（Runge-Kutta）方法，是间接利用泰勒展开的思想构造的一类数值方法。在讨论龙格－库塔方法之前，我们先介绍泰勒方法。9.2-1泰勒方法假定初值问题(9.1)，(9.2)的解yx及函数,fxy是足够光滑的，利用p阶泰勒公式得2112112!!,,2!,!nnnnpppnnnnnnpppnnhyxyxhyxyxhyxOhphyxhfxyxfxyxhfxyxOhp当h充分小时，略去余项1pOh，则可导出近似值1ny的计算公式：2110,,2!,!nnnnnnppnnhyyhfxyfxyhfxypya(9.18)算式中用到的各阶导数1,iinnnyfxy可以由微分方程的右端函数计算出来。例如82,,,,,,,2nnnnnnnnxnnnnynnnnnxxxyyyyxyyfxyyfxyfxyfxyfxyyfxyfffffff（9.18）式称为p阶泰勒方法。特别地，当1p时（9.18）式就是欧拉公式（9.3）；当2p时，得二阶泰勒方法1,2nnnnxyhffyyhffxy（9.19）［例2］用泰勒方法解223(0)1dyxdxyy分别用二阶、四阶泰勒方法计算点nx＝0.1,0.2,…,1.0处的数值解，并与精确解进行比较。（请完成）9.2-2龙格-库塔方法的基本思想与二阶公式的推导龙格-库塔方法的基本思想是利用,fxy在某些点处的值的线性组合，来构造一类计算公式，使其按泰勒展开后与初值问题的解的泰勒展式比较，在尽可能多的项完全相同以确定其中的参数，从而保证算式有较高的精确度。由于避免了在算式中直接用到,fxy的导数，所以说龙格-库塔方法是基于间接利用泰勒展开的思想。一般龙格-库塔方法的形式为9111221222113331132211,11...(,)(,)(,)...(,...)nnrrnnnnnnrnrnrrrryykkkkhfxykhfxahykkhfxahykkkhfxahykk(9.20)其中，i，ij，i为常数。选取这些常数的原则，是要求(9.20)第一式右端在(,)nnxy处作泰勒展开，且按h的幂次整理得23112311...2!3!nnyyhhh与微分方程解的泰勒展式23111()()...2!3!nnnnnyxyxfhfhfh有尽可能多的项重合，即要求123,,,...nnnvfvfvf通常把（9.20）式叫做r级（或r段）的计算公式。如果(1)ffnvf，对1,2,...,jm成立，而对1jm时不成立，则所得的公式称为m阶的。下面我们以二阶龙格-库塔公式为例，进行具体的推导说明。设想构造如下形式的公式：11122122211(,)(,)nnnnnnyykkkhfxykhfxahyk（9.21）要求适当选取系数1，2，和，使当()nnyyx时，(9.21)式的局部截断误差为3()Oh。为此，将2k在(,)nnxy展开，有22n=[+(+)+O()]xynkhfhfffh将上式代入(9.21)式，整理得231122()(+)+O()nnnxynyyhfhfffh(9.22)10在解()yx的泰勒展式231()()()()+O()2nnnnhyxyxhyxyxh中，()(,())yxfxyx，()(,)(,)(,)xyxyyxfxyfxyfxyfff，即有231()()()+O()2nnnxynhyxyxhffffh(9.23)将(9.22)式与(9.23)式比较，当nny=y(x)时，只需取122211212(9.24)则(9.21)式的局部截断误差为3O()h。(9.24)式中包含有四个未知数，只有三个方程，其中有一个未知数可以任意取值，且显然。(9.24)式的每一组解