常微分方程数值解法的误差分析教材Word格式.docx
《常微分方程数值解法的误差分析教材Word格式.docx》由会员分享，可在线阅读，更多相关《常微分方程数值解法的误差分析教材Word格式.docx（18页珍藏版）》请在冰豆网上搜索。
（一）、欧拉法5
（二）、向后Euler方法6
（三）、0-法7
（四）、改进欧拉法7
（五）Rung*Kutta方法9
2、线性多步法14
总结16
参考文献：
17
自然界中很多事物的运动规律可用微分方程来刻画。
常微分方程是研究自然科学和社会科学中的事物、物体和现象运动、演化和变化规律的最为基本的数学理论和方法。
物理、化学、生物、工程、航空航天、医学、经济和金融领域中的许多原理和规律都可以描述成适当的常微分方程，因此，常微分方程的理论和方法不仅广泛应用丁自然科学，而且越来越多的应用丁社会科学的各个领域。
由丁该问题比较复杂且涉及的面广，使得有些问题的解析解很难求出，而对丁一些典型的微分方程（如线性方程、某些特殊的一阶非线性方程等）可以运用基本方法求出其解析解，并在理论上可以根据初值问题的条件把其中的任意常数完全确定下来。
然而，在生产实际和科学研究中所遇到的微分方程往往很复杂，在很多情况下都不可能给出解的解析表达式，有时即使能求出形式的解，也往往因计算量太大而不实用，而且高次代数方程求根也并不容易，所以用求解析解的方法来计算微分方程的数值解往往是不适宜的。
从实际意义来讲我们更关心的是某些特定的自变量在某一个定义范围内的一系列离散点上的近似值。
本文研究的主要是针对常微分方程各种数值解法的误差进行分析。
一、常微分方程
1、定义
首先，我们在这部分给出所需的一些基本概念和基本知识。
我们已经知道微分方程就是联系着自变量、未知函数以及其导数的关系式。
如果在微分方程中，自变量的个数只有一个，我们称这种微分方程为常微分方程。
方程
变yb曳cy=ftdt2dt
就是常微分方程的例子，这里y是未知函数，t是自变量。
微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶数。
2、常微分方程初值问题描述
在自然科学和经济的许多领域中，常常会遇到一阶常微分方程的初值问题
dy
—=fx,y
』dx'
yj,a4x<
b
（1）
yg）=y（）
这里f（x,y）是充分光滑，即关丁x或y满足李普希茨条件的二元函数，y0是给定的初值，y（xo）=y°
称为初始条件。
3、数值解法的基本思想与途径
一阶微分方程的初值问题
（1）的解y（x）是区间Ia,b］上的连续变量x的函数，因而问题
（1）实际上是一个连续性的问题，求这个问题的数值解，就是要求在区间la,b］上的若干个离散点处的函数近似值，例如：
a《x0<
为<
...<
xn《b,
然后计算出解y（x）的近似值
yx。
yxi,...,yxn.
股常取xo,xi,...,xn为等距离的点，即
为一x。
=x2-为=...=xn—xn」=h
或
称h为步长
xi=ai,h^0,1,...n,
建立数值方法的第1步，就是把连续性问题
（1）通过一定的方法化为在给定的n+1个点上的近似的差分方程的初值问题，称这个过程为离散化。
常用离散化的方法如下：
（1）用差商替代导数
在点x处的导数y^xi）可以近似地表示成差商
y'
（x片
从而把初值问题
（1）化为差分问题
yXo=ys
其中％表小解y（x）在点x处的近似解，即yi=y（x）。
当然，用差商来近似地表示导数，方法不是唯一的，这里所用的是所谓的向前差商。
（二）Taylor展开法
在一点（例如点x）的附近，y（x）的同次数的近似多项式中的Taylor多项式
p
y（Xi+h）女y（x）+hy'
（为）+...+
土ypxip!
为最好。
其中p为一正整数。
通过微分方程y'
=f（x,y）,便可以逐次把各阶导数y；
y"
IH在x处的值表示出来。
（三）数值积分法
对微分方程y'
=f（x,y侄区问仅,X」上求积分，得
xi1
y（x.）—y（xi）=七f（x,y（x）px,j=0,[,“
丁是，初值I可题
（1）便可以近似地化为
yi=y
yxo=y。
xifx,yxdx'
i=0,1川n.
这样，关丁上式右端的积分，可以用数值积分方法计算其近似值。
4、数值解的分类
常微分方程初值问题的数值解法一般分为两大类：
单步法：
所谓单步法是指这类方法在计算yn中时，只用到前一步的值x^4i,xn,yn然后逐步往下计算。
这个算法的代表是龙格---库塔算法，简称R—k方法。
四阶显示Runge-Kutta方法是求解普通常微分方程初值问题数值解法中的重要方法，而隐式Runge-Kutta公式是求解刚性常微分方程初值问题的重要方法。
多步法：
这类方法在计算yn+时，除了用到前一步的值xn书,xn,yn，之外，还要用到
xn*y^p=1,2,HI,k;
k0
这前面k步的值，这个算法的代表就是阿达姆斯（Adams）方法。
5、问题
（1）解的存在惟一性定理
一个常微分方程是不是有特解呢？
如果有，乂有几个呢？
这是微分方程论中一个基本的问题，数学家把它归纳成基本定理，叫做存在和唯一性定理。
因为如果没有解，而我们要去求解，那是没有意义的；
如果有解而乂不是唯一的，那乂不好确定。
因此，存在和唯一性定理对丁微分方程的求解是十分重要的。
这个重要的存在和唯一性就是下面列出的著名的存在惟一性定理。
定理"
如果f（x,y）在带形区域R={（x,y）|a壬x^b,-*<
y<
+=«
}中连续，且关丁y满足Lipchiz条件：
即存在正常数L,使得
f（xy）—f,x2^<
|Li-y|2,y
对所有的xw〔a,b】以及任何y”y2都成立，那么初值问题
（1）存在惟一的连续可微解y=y（x）。
二、几种常用的数值解法及其误差分析
首先来了解下初值问题数值解中最简单的一种方法一一欧拉法。
1、单步法
Euler折线法发生的历史背景。
在微分方程研究之初，瑞士数学家L.Euler（1707.4-1783.9）做出了开创性的工作。
他和其他一些数学家在解决力学、物理学问题的过程中创立了微分方程这门学科。
在常微分方程方面，Euler在1743
年发表的论文中，用代换y=e^扩给出了任意阶常系数线性微分方程的古典解法，最早引入了“通解”和“特解’’的概念。
1768年，Euler在其有关月球运行理论的著作中,创立了广泛用丁求初值问题
（1）的数值解的方法，次年乂把它推广到二阶方程。
欧拉的想法如下：
选择步长h>
0，然后在x0《x£
x0+h情况下用解函数的切线
lx=yx-°
xf,0x0y
代替解函数。
这样对丁点x1=x0+h就可以得到
y1=y0hfx°
y0
在点（xi,y1）重复如上的程序再次计算新的方向，就会得到所谓的递推公式：
xm41=xm+h,ymF=ym+hf（xm,ym）,
这就是Euler方法。
由此，再通过连接所有这些切线得到的函数被称为Euler折线。
如果我们令h->
0,这些折线就会越来越接近解函数。
（一）、欧拉法
Euler方法是最简单的一步法，它是一阶的，精度较差，但公式很简单，即
yn1=ynhyXn,ynn=0,1j。
（3）
Euler方法的几何意义在数值计算思想中已经体现出来了，实际上就是用过已
知点的折线来近似代替过此点的积分曲线。
因此，这种方法乂称为折线法。
在Euler
法中，数值解的误差首先是由差商代替导数引起的，这种近似替代所产生的误差称为截断误差。
另外，计算过程中还会由丁数值的舍入产生另一种误差一一舍入误差。
显然只有当初产生的误差在以后各步的计算中不会无限制扩大时，即当初始误差充分小时，以后各步的误差也可以充分小，Euler法才具有实用价值。
收敛
性、截断误差估计与稳定性闷题是常微分方程各种数值解法研究中必须考虑的基本问题。
显然这些问题在Euler法中是得到验证的，详见下面例子分析。
在xy平面上，微分方程
（1）的解y=y（x）称作它的积分曲线。
积分曲线上一点（x,y）线斜率等丁函数f（x,y）的值，如果按函数f（x,y）在xy平面上建立一个方向场，那么，积分曲线上每一点的切线方向均与方向场在该点的方向相一致，基丁上述几何解释，从初始点Po（x0,yo）出发，先依方向场在该点的方向推进到x=x〔上一点臼，然后再从p1依方向场的方向推进到x=x2上一点p2,循此前进推出一条折线
P0P1P2…，一般地，设已做出该折线的顶点Pn，过Pn（X,yn）依方向场的方向再推进到Pn41（xn+,yn+），显然两个顶点Pn,Pn*的坐标有关系
yn1=ynhfxn,yn
这就是著名的欧拉（Euler）公式。
若初值y°
已知，则依公式⑷可逐步算出
V1=*hfx°
yo,
V2=V"
hf*,乂,
III
例1求解初值问题
y=y_2x/y,（0：
：
x:
:
1）
y（0）二1,
解欧拉公式的具体形式为
取步长h=0.1,计算结果如下表:
表1计算结果对比
初值问题（2.1.2）有解y=5/^芯，按这个解析式子算出的准确值y（xn）同近似值yn一起歹0在表1,两者相比较可以看出欧拉方法的精度很差。
（二）、向后EuIer方法
向后Euler方法和Euler方法差不多，只是把y'
（x"
）用
yXn1-yXn
h
去代替，这时计算公式为
（6）
yn1=yhfxn1,yn1,NJ
（n=0,1,2川）
yx°
=y°
向后Euler方法的总体截断误差也是一阶的，因此向后Euler方法是收敛的。
这里需要指出它与Euler方法的一个很大不同之处，Euler方法是显式方法，即
由xn41,xn,yn,h明显地表示出来了，而向后Euler方法是隐式方法，计算yn41时要解隐式方程（6）。
通常解此方程用迭代法。
因此计算较为麻烦，但比显式Euler方法精度要高。
（三）、8-法
将Euler方法公式与向后Euler方法公式作加权平■均，得到如下公式：
fyJ=yn+h丐f（xn,yn广（1顼）f（x忡，）]。
[刃巾⑺
=0,|,2）
y冷'
）=%,
称式⑺为初值问题⑴的8-法公式。
其中如y（x0）=y0为初值条件。
1-
在此法中当6=-时，即
2h
yn^=yn+2」（xn,Yn）+f（x”甲，Yn*）]（8）
此时的e一法称为梯形公式法。
梯形公式也隐式格式，用起来要进行迭代，其计算公式为
y"
）=mn,yn）+f（XmynQ）k=0,i,2川）⑼
[ynJ°
）=yn+hf（xn,yn）,
这里在应用本迭代法时，是先用Euler方法求初值yn+1的近似值yn+i（°
）即：
yn+i（°
Lyn+hf（Xn,yn）
然后将yn+1（°
陪代梯形公式（8）中的yn+1得到的式（9）。
式⑼乂称为预测校正公式。
换言之，由Euler方法给出预测值，再用梯形法予以校正。
很显然，当步长办取得适当小时，由Euler方法算出的值已是较好的近似。
格式（9）收敛很快，通常只需一两次迭代即可满足精度要求，若需多次迭代,
则应缩小步长h后再行计算。
梯形公式法比用向后Euler方法的迭代步长可以放宽一倍，它的总体截断误差为O（h2），比Euler方法高一阶。
但它每积分一步要计算二次函数值，这说明的精度的提高是以增加计算量为代价的。
（4）
xn中积分，得
（10）
、改进欧拉法
为得到比欧拉法精度高的计算公式，在等式如果对方程
（1）从xn到
xn1
yxn1=yxnft,ytdt
xn
右端积分中若用梯形求积公式近似，并用％代替y（Xn）,yn书代替y（XnG，则得h....
yn1=yn2II」Xn,ynfXn1,y”1（11）
称为改进欧拉法.
改进欧拉方法是隐式单步法，可用迭代法求解.用欧拉方法提供迭代初值，则改进欧拉法的迭代公式为
‘乂拦=人+；
[/（"
吊）+了0心丈2）]
（12）
侬=。
"
成
为了分析迭代过程的收敛性，将（11）式与（10）相减，得
=5仇fm+l）-/（w，照）]M
于是有
Am-欢倍竽*『照
式中L为f（x,y）对y满足Lipschitz常数，如果选取h充分小，使得
hL*
—<
1
2,
则当kT8时有y：
%T外十，这说明迭代过程（12）是收敛的.
例2用改进的欧拉方法求解初值问题
（1）.
解改进的欧拉公式为
*=人+第（L»
）+/g『璀）]
2
（上=。
12…）.
仍取h=0.1，计算结果见下表.同例1中欧拉法的计算结果比较，改进欧拉法明显改善了精度.
表2计算结果对比
%
旭）
所）|
0.1
0.2
0.3
0.4
0.5
1.0959
1.1841
1.2662
1.3434
1.4164
1.0954
1.1832
1.2649
1.3416
1.4142
0.6
0.7
0.8
0.9
1.0
1.4860
1.5525
1.6153
1.6782
1.7379
1.4832
1.5492
1.6125
1.6733
1.7321
（五）Rung「Kutta方法
德国数学家C.D.T.Runge（1851927）是数值方法发展史上具有里程碑
作用的人物。
1895年，他在Hanover发表了关丁微分方程数值解法的经典论文《常
微分方程数值解法》。
此文成为常微分方程RungaKutta方法的发端。
此后，Runge结合教学活动积极投身丁发展一般的数值分析特别是各种实际应用中的Runge-
Kutta方法（严格来说，此方法在Kutta作出工作后才能称作RungbKutta方法）。
RungbKutta方法是一种特殊的单步方法，事实上，这个方法可以看作在（Xm,XmG上取若干条积分曲线的若干个点的切线斜率，再进行一次（或多次）算术（或加权）
平均后产生的新斜率，再按这个斜率从（Xm,ym）出发，以直线带曲线向前推进一步
的过程。
与Taylor展示法相比，RungbKutta方法不用增加微商f（x,y）的次数就可以得到较高的阶。
Runge—Kutta方法除了在微分方程求解中扮演的传统角色外，人们发现相关类型的初值问题可以用Runge一Kutta方法或适合更一般问题的
RungbKutta方法求解，比如RungbKutta方法被应用到了Hamilton系统中。
前面提到的几种数值解法的精度是很低的，下面给出高阶一步法一一Runge-Kutta方法。
它是最常用的一种数值解法，因为它相当精确、稳定、容易编程。
Runge—Kutta方法至今仍然得到广泛地应用。
1、二级二阶RungbKutta方法
由RungbKutta方法的思想，我们得到二阶RungbKutta公式为
=yn+4（Ki+3K2）
Ki=hfXn,yn,
（22）
（13）
K2=hf.Xn+—h,yn+—Ki.
33
yn*=yn+*Ki+K2）
』Ki=hf（Xn,%）,K2=hf（Xn+h,yn+Ki）
（i4）
RungbKutta公式是在计算两次函数值的情况下，局部截断误差的阶最高是3,
式（i3）是允许函数f（x,y）任意变化情况下截断误差最小的二阶方法。
要再提高阶就必须增加计算函数值的次数。
上述式（i4）乂称为欧拉预估一一校正公式。
2、三级三阶RungbKuuta方法
两个常用的三阶RungbKuuta方法分别为：
iyn*=yn+&
h（2Ki+3K2+4K3）
Ki=f（Xn,%）
ii
K2=fXn+，h,yn+，hKi,
I22）
/工3工3）
（i5）
K3=fXn+—h,yn+—hK2,i44
i,、
yn±
=yn+—h（Ki+4K2+K3）,6
Ki=fXn，yn,
Jii）
K2=fXn+—h,yn+—hKi,
I22J
（i6）
K3=fXnh,yn-hKi2hK2.
这两个常用方法在解决实际问题中能够达到较低的精度要求。
但是要更高精度要
求的，我们必须了解更高阶的方法一四级四阶Rung*Kuuta方法。
3、四级四阶RungbKutta方法
这种方法在解决实际问题中常用。
在这里，我们将公式进行导出，以熟悉其
方法的实用过程。
由丁RungesKutta公式对初值问题（i）中的一般f（x,y）都适用,则它必然对特殊的f（x,y）=Xy,f（x,y）=X或f（x,y）=y也适用。
从而可以定出特定
的参数来，因此，这里采用的是一种待定系数法，来导出四级四阶显式
Kutta公式，过程十分简单明了。
为了计算简单，令
=Xnh,
y（xn）近似值为yn,这样求解问题
（1）的四阶显示Runge—Kutta公式为
=忘+顺”+小2+«
+cK）,
Ki=f（Xn，yn）,
』K2=f（Xn+ha2,yn+hb2iKi）,
K3=f（Xn+ha3,yn+hb3iKi+hb43K2）,
K4=f（Xn+ha4,yn+hb4iKi+hb43L）,
这里要求
a2=b2i,a^~b3i炫2,a4=b4i'
b42'
b43-
由丁RungaKutta公式的思想与数值积分
Xni
fX,ydx
Xi
相似。
所以希望四阶显示RungbKutta公式对
y=X3,y0=0
是精确的。
这样，将（i7）式应用于（i9）式，就有
Ki=Xn3,
Runge—
（17）
（18）
（19）
（20）
K2a?
）3=Xn33用端Mx^"
a23,
_3
K4=（Xna4）
K3=（XnW3）3=Xn33Xn2a33x^2a33,=x3■3xn2a4■3xna42■a43,
yni=yn（Ci
i2
-c2a2ca
c2c3c4）ynhc2a2c3&
c4&
yh
%+ca42）*%3-（ca/32Cafc’a）y34h（.4）
6
与Taylor公式
（21）
相比较，有
升级会员 