数值分析：第九章 常微分方程数值解

1第九章

常微分方程数值解2常微分方程初始问题

许多科学技术问题，例如天文学中的星体运动，空间技术中的物体飞行，自动控制中的系统分析，力学中的振动，工程问题中的电路分析等，都可归结为常微分方程的初值问题。

常微分方程中往往只有少数较简单和典型的常微分方程（线性常系数常微分方程等）可求出其解析解。本章考虑的是一阶微分方程（组）的数值解。高阶的常微分方程可以转化成一阶微分方程组。3举例例：单摆运动问题OAm可以转化成一阶微分方程组4常微分方程解的存在唯一性

定理

考虑一阶常微分方程的初值问题对任意定义在上的都成立，则上述方程存在唯一解。只要在上连续,且关于y满足Lipschitz条件，即存在与无关的常数L使5常微分方程的数值解

数值解的基本思想满足微分方程

常微分方程等价于下列积分方程

常微分方程初值问题的数值解法，是要寻求解函数y(x)在一系列节点

处的近似值

关键在于导数如何近似。

关键在于积分如何近似。6初始问题的欧拉方法

初值问题的数值解法一般按节点从左至右的顺序依次求出的近似，所以称为步进法

单步法

从初值开始,依次求出,后一步的值只依靠前一步的来计算典型的单步法是Euler(欧拉)方法,其计算格式是:

多步法在计算时需要前面若干步值来计算。7欧拉方法的推导

以Euler法为例说明初值问题数值方法的三种基本途径

Taylor展开法

忽略高阶项,取近似值可得到Euler公式

数值微分法,用差商代替导数

8欧拉方法的推导则

数值积分法

将区间积分右端采用左矩形公式数值积分，得9几个简单的数值方法

显式欧拉法

隐式欧拉法方法上是向后差商近似导数，或右矩形数值积分得到

由于未知数yn+1同时出现在等式的两边，不能直接得到，故称为隐式/*implicit*/欧拉公式.

一般先用显式计算一个初值，再迭代求解。10几个简单的数值方法

梯形公式

显、隐式两种算法的平均或梯形公式数值积分得到

改进的欧拉公式

先用显式欧拉公式作预测，算出),(nnnyxfhy+=1+ny

再将代入隐式梯形公式的右边作校正，得到1+ny)],(),([211++++=nnnnnxfyxfhyy1+ny11单步法的误差分析

整体截断误差

数值解和精确解之差

整体截断误差除与步计算有关外,还与的计算有关，所以一般很难分析清楚。

局部截断误差一般单步法有增量形式

在假设yn=y(xn)，即第n步计算是精确的前提下，考虑截断误差Tn+1=y(xn+1)

yn+1称为局部截断误差

12算法的精度若某算法的局部截断误差为O(hp+1),则称该算法有p阶精度。显式欧拉，隐式欧拉为1阶精度

梯形公式和改进欧拉公式为2阶精度

显式欧拉局部截断误差梯形公式局部截断误差13例例：考察初值问题在区间[0,0.5]上的解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。0.00.10.20.30.40.5精确解改进欧拉法

欧拉隐式欧拉显式

节点xi1.0000

2.00004.0000

8.00001.6000101

3.2000101

1.00002.5000101

6.25001021.56251023.90631039.76561041.00002.50006.25001.56261013.90631019.76561011.00004.97871022.47881031.23411046.14421063.0590107从表中可以看出，这几种方法计算误差比较大，说明这些方法精度不够，需寻找更高精度的方法14高精度方法

基于Taylor展开式可以建立任意精度的单步法

而所以，可以构造格式

这种格式使用到了各阶偏导数，使用不便。15Runge-Kutta方法

Runge-Kutta方法的基本思想

只要对平均斜率提供一种算法，便相应地得到一种微分方程的数值计算公式。16Runge-Kutta方法

考察改进的欧拉法可以将其改写为：

改进欧拉公式却是利用了xn与xn-1两个点的斜率值k1=f(xn,yn)与k2=f(xn+1,yn+hk1)取算术平均作为平均斜率的近似值。17Runge-Kutta方法启示设法在[xk,xk+1]上多预报几个点的斜率值，然后将它们加权平均作为平均斜率的近似值，则有可能构造出更高精度的计算公式推广其中

由f(x,y)在一些点的值加权平均所得这种方法称为Runge-Kutta方法,简记为R-K公式.若

由r个f(x,y)的值加权平均所得的公式称为r级R-K公式

18二级Runge-Kutta方法

二级Runge-Kutta方法),(),(][12122111phKyphxfKyxfKKKhyyiiiiii++==++=+ll

希望能确定系数

1、

2、p，使得到的算法格式有2阶精度，即在的前提假设下，使得

将K2在(xi,yi)点作Taylor展开19二级Runge-Kutta方法

将K2代入第1式，得到

将yi+1与y(xi+1)在xi点的泰勒展开作比较20二级Runge-Kutta方法要求，则必须有：

这里有3个未知数，2个方程,解不唯一

就是改进的欧拉法。

就是中点公式。

21三阶与四阶显式R-K方法

三级Runge-Kutta方法其中及均为待定参数.

局部截断误差为

只要将按二元函数泰勒展开，使，22三阶与四阶显式R-K方法可得待定参数满足方程

这是8个未知数6个方程的方程组，解也不是惟一的.23三阶与四阶显式R-K方法一个特殊的方法

同样的方法，可以得到经典4阶R-K方法24注更一般R-K方法

龙格-库塔法的主要运算在于计算的值，即计算f的值。Butcher于1965年给出了计算量与可达到的最高精度阶数的关系：25注753可达到的最高精度642每步须算Ki的个数

由于龙格-库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好采用低阶算法而将步长h

取小。26

例

解设取步长，从到用四阶龙格-库塔方法求解初值问题这里，经典的四阶龙格-库塔公式为2728F='y-2*x/y';a=0;b=1;h=0.1;n=(b-a)/h;X=a:h:b;Y=zeros(1,n+1);Y(1)=1;fori=1:nx=X(i);y=Y(i);

K1=h*eval(F);

x=x+h/2;y=y+K1/2;K2=h*eval(F);x=x;y=Y(i)+K2/2;K3=h*eval(F);x=X(i)+h;y=Y(i)+K3;K4=h*eval(F);Y(i+1)=Y(i)+(K1+2*K2+2*K3+K4)/6;end2930

求微分方程（组）的解析解命令:dsolve(‘方程1’,‘方程2’,…‘方程n’,‘初始条件’,‘自变量’)

求微分方程（组）的数值解命令

[t，x]=solver（’f’,ts,x0,options）ode45ode23ode113ode15sode23s由待解方程写成的m-文件名ts=[t0，tf]，t0、tf为自变量的初值和终值函数的初值函数值用于设定误差限(缺省时设定相对误差10-3,绝对误差10-6),命令为：options=odeset（’reltol’,rt,’abstol’,at）,rt，at：分别为设定的相对误差和绝对误差.31单步法的收敛性定义若某算法对于任意固定的x=x0+nh，当h0

(即n

)时有yn

y(xn)，则称该算法是收敛的。

例

就初值问题考察欧拉显式格式的收敛性。解：该问题的精确解为欧拉公式为对任意固定的x=xn=nh

，有

32定理:对初值问题的单步法,若局部截断误差为,且函数

对y满足Lipschitz条件,即存在L>0,使得对一切成立,则该方法收敛,且有单步法的收敛性收敛性的条件欧拉，改进欧拉以及Runge-Kutta方法都是收敛的

整体截断误差比局部截断误差低一阶33

对于欧拉方法，由于其增量函数就是，故当关于满足利普希茨条件时它是收敛的.

改进的欧拉方法，其增量函数由欧拉方法的收敛性

改进欧拉方法的收敛性

假设关于满足利普希茨条件，记利普希茨常数为,34单步法的稳定性

考虑收敛性时，当h0时有yn

y(xn)，此时yn为理论数值解，也就没有考虑到计算过程中的舍入误差

当考虑计算过程中舍入误差对最后结果是否有影响，所以引入数值算法绝对稳定的概念定义

若一种数值方法在节点值上大小为扰动，于以后各节点值上产生的偏差均不超过，则称该方法是绝对稳定的.由于f(x,y)比较复杂，给稳定性研究带来困难，所以考虑试验方程35单步法的稳定性

试验方程

该试验方程为线性方程，其精确解是稳定的。若对该方程数值算法都稳定，则对其它方程也靠不住。对一般非线性方程，局部可以化成线性的来近似

单步法的稳定性

单步法用于试验方程，若得到的解为满足

则是绝对稳定的。36单步法的稳定性记，，的范围叫绝对稳定区域

hlh=h绝对稳定区域与实轴之交叫绝对稳定区间

方法A比方法B的绝对稳定区域大称A比B稳定显式欧拉法的稳定性0-1-2ReImgh在复平面上是以(-1,0)为圆心的单位圆绝对稳定区间为【-2,0】37隐式欧拉法的稳定性210ReImg所以绝对稳定区域为左半平面绝对稳定区间隐式欧拉比显式欧拉稳定性好38梯形公式的稳定性所以绝对稳定区域为左半平面绝对稳定区间39经典4阶龙格-库塔方法的稳定性40

研究一般方程的稳定性时，相当于模型方程中的λ

41线性多步法

一个特例对于初值问题若用中心差商离散导数称为二步法42线性多步法

推广当

10时，为隐式公式;

1=0则为显式公式。

43线性多步方法的构造

基于Taylor展开的构造方法: