常微分方程数值解法的误差分析.docx
《常微分方程数值解法的误差分析.docx》由会员分享，可在线阅读，更多相关《常微分方程数值解法的误差分析.docx（15页珍藏版）》请在冰豆网上搜索。
常微分方程数值解法的误差分析
师大学
2021届学士学位论文
常微分方程数值解法的误差分析
学院、专业数学科学学院数学与应用数学
研究方向
学生姓名翅
学号
指导教师」
指导教师职称讲师
年月日
常微分方程数值解法的误差分析
娜
〔师大学数学科学学院，，235000J
摘要
自然界与工程技术中的很多现象，往往归结为常微分方程定解问题。
许多偏微分方程问题也可以化为常微分方程问题来近似求解。
因此，研究常微分方程的数值解法是有实际应用意义的。
数值解法是一种离散化的数学方法，可以求出函数的准确解在自变量一系列离散点处的近似值。
随着计算机计算能力的增强以及数值计算方法的开展，常微分方程的数值求解方法越来越多，比拟成熟的有Euler法、后退Euler法、梯形方法、Runge-Kutta方法、投影法和多步法，等等.本文将对这些解的误差进展分析，以求能够得到求解常微分数值解的精度更好的方法。
关键词：
常微分方程，数值解法，单步法，线性多步法，局部截断误差
ErrorAnalysisofNumericalMethodforSolvingthe
OrdinaryDifferentialEquation
LiNa
（SchoolofMathematicalScienee,HuaibeiNormalUniversity,Huaibei,235000）
Abstract
Innatureandengineeringhavemanyphenomena,definitesolutionoftheproblemoftenboilsdowntoordinarydifferentialequations.Sostudythenumericalsolutionofordinarydifferentialequationsispracticalsignificanee.Thenumericalmethodisadiscretemathematicalmethods,andexactsolutionofthefunctioncanbeobtainedintheapproximationofaseriesofdiscretepointsoftheargument.Withtheenhancedputingpowerandthedevelopmentofnumericalmethods,ordinarydifferentialequationshavemoreandmorenumericalsolution，therearesomematuremethods.SuchasEulermethod,backwardEulermethod,trapezoidalmethod,Runge-Kuttamethod,project!
onmethodandmulti-stepmethodandsoon.Therefore,numericalsolutionofdifferentialequationisofgreatpracticalsignificanee.Throughthispaper,errorofthesesolutionswillbeanalyzedinordertogetatheaccuracybetterwaytosolvethenumericalsolutionofordinarydifferential.
Keywords:
Ordinarydifferentialequations,numericalsolutionmethods,singlestepmethods,linearmulti-stepmethods,localtruncationerror
引言1
一、常微分方程1
1、定义1
2、常微分方程初值问题描述2
3、数值解法的根本思想与途径2
4、数值解的分类3
5、问题⑴解的存在惟一性定理4
二、几种常用的数值解法及其误差分析4
1、单步法4
（一）、欧拉法5
（二）、向后Euler方法6
（三）、法7
（四）、改良欧拉法8
（五）Runge—Kutta方法9
2、线性多步法14
总结17
参考文献：
17
引言
自然界中很多事物的运动规律可用微分方程来刻画。
常微分方程是研究自然科学和社会科学中的事物、物体和现象运动、演化和变化规律的最为根本的数学理论和方法。
物理、化学、生物、工程、航空航天、医学、经济和金融领域中的许多原理和规律都可以描述成适当的常微分方程，因此，常微分方程的理论和方法不仅广泛应用于自然科学，而且越来越多的应用于社会科学的各个领域。
由于该问题比拟复杂且涉及的面广，使得有些问题的解析解很难求岀，而对于一些典型的微分方程（如线性方程、某些特殊的一阶非线性方程等）可以运用根本方法求岀其解析解，并在理论上可以根据初值问题的条件把其中的任意常数完全确定下来。
然而，在生产实际和科学研究中所遇到的微分方程往往很复杂，在很多情况下都不可能给出解的解析表达式，有时即使能求岀形式的解，也往往因计算量太大而不实用，而且高次代数方程求根也并不容易，所以用求解析解的方法来计算微分方程的数值解往往是不适宜的。
从实际意义来讲我们更关心的是某些特定的自变量在某一个定义围的一系列离散点上的近似值。
本文研究的主要是针对常微分方程各种数值解法的误差进展分析。
一、常微分方程
K定义
首先，我们在这局部给出所需的一些根本概念和根本知识。
我们已经知道微分方程就是联系着自变量、未知函数以及其导数的关系式。
如果在微分方程中，自变量的个数只有一个，我们称这种微分方程为常微分方程。
方程d2y
就是常微分方程的例子，这里y是未知函数，t是自变量。
微分方程中岀现的未知函数最高阶导数的阶数称为微分方程的阶数。
7（xo）=^
其中X表示解）心）在点兀处的近似解，即乃=）0。
当然，用差商来近似地表示导数，方法不是唯一的，这里所用的是所谓的向前差商。
（-）Taylor展开法
在一点（例如点兀）的附近，y（x）的同次数的近似多项式中的Taylor多项式y（xi+h）^y（xi）+hy\xi）+...+
为最好。
其中p为一正整数。
通过微分方程yf=f（^y）,便可以逐次把各阶导数*在兀处的值表示岀来。
（三）数值积分法
对微分方程y=f（^y）在区间山,和]上求积分，得
）心+J-yM=f（x，y⑴亦<•=o,i.
于是，初值问题
（1）便可以近似地化为X+]=X+J「/（x，y（x））dx,.
<也2=0丄••讥
）Oo）=）b
这样，关于上式右端的积分，可以用数值积分方法计算其近似值。
久数值解的分类
常微分方程初值问题的数值解法一般分为两大类:
单步法：
所谓单步法是指这类方法在计算儿+1时，只用到前一步的值x“+i，£，y“，然后逐步往下计算。
这个算法的代表是龙格…库塔算法，简称r—k方法。
四阶显示Runge-Kutta方法是求解普通常微分方程初值问题数值解法中的重要方法，而隐式Runge-Kutta公式是求解刚性常微分方程初值问题的重要方法。
多步法：
这类方法在计算儿昭时，除了用到前一步的值£+|，兀,儿,之外，还要用到
X"片卡（/2=1,2….**>0）
这前面R步的值，这个算法的代表就是阿达姆斯（Adams）方法。
5、问题⑴解的存在惟一性定理
一个常微分方程是不是有特解呢?
如果有，又有几个呢?
这是微分方程论中一个根本的问题，数学家把它归纳成根本定理，叫做存在和唯一性定理。
因为如果没有解，而我们要去求解，那是没有意义的；如果有解而又不是唯一的，那又不好确定。
因此，存在和唯一性定理对于微分方程的求解是十分重要的。
这个重要的存在和唯一性就是下面列出的著名的存在惟一性定理。
定理山如果在带形区域R=｛（x,y）|^-oo即存在正常数L,使得
|/（匕必）一/（32）|"卜］一力|，
对所有的x^\a,b\以及任何必，）転都成立，那么初值问题
（1）存在惟一的连续可微解y二）‘（兀）。
二、几种常用的数值解法及其误差分析
首先来了解下初值问题数值解中最简单的一种方法——欧拉法。
K单步法
Euler折线法发生的历史背景。
在微分方程研究之初，瑞士数学家L.Euler（1707.4-1783.9）做岀了开创性的工作。
他和其他一些数学家在解决力学、物理学问题的过程中创立了微分方程这门学科。
在常微分方程方面，Euler在1743年发表的论文中，用代换y=占扩给出了任意阶常系数线性微分方程的古典解法，最早引入了“通解"和“特解''的概念。
1768年，Euler在其有关月球运行理论的著作中，创立了广泛用于求初值问题
（1）的数值解的方法，次年又把它推广到二阶方程。
欧拉的想法如下：
选择步长力>0,然后在x0心）=北+（兀一勺）/（心北）
代替解函数。
这样对于点可=珀）+力就可以得到
）1=儿+妙'（兀，％）・
在点（勺，X）重复如上的程序再次计算新的方向，就会得到所谓的递推公式：
=兀”+力，儿屮=）1+/?
广（兀"儿
这就是Euler方法。
由此，再通过连接所有这些切线得到的函数被称为Euler折线。
如果我们令"TO,这些折线就会越来越接近解函数。
（一）、欧拉法
Euler方法是最简单的一步法，它是一阶的，精度较差，但公式很简单，即儿+1=几+力）'（禺，儿）W=0,1,2...）（3）
Euler方法的几何意义在数值计算思想中已经表达岀来了，实际上就是用过点的折线来近似代替过此点的积分曲线。
因此，这种方法又称为折线法。
在Euler法中，数值解的误差首先是由差商代替导数引起的，这种近似替代所产生的误差称为截断误差。
另外，计算过程中还会由于数值的舍入产生另一种误差一一舍入误差。
显然只有当初产生的误差在以后各步的计算中不会无限制扩大时，即当初始误差充分小时，以后各步的误差也可以充分小，Euler法才具有实用价值。
收敛性、截断误差估计与稳定性闷题是常微分方程各种数值解法研究中必须考虑的根本问题。
显然这些问题在Euler法中是得到验证的，详见下面例子分析。
在Q平面上，微分方程⑴的解尸y（x淋作它的积分曲线。
积分曲线上一点（忑刃线斜率等于函数的值，如果按函数/（兀刃在Q平面上建立一个方向场，那么，积分曲线上每一点的切线方向均与方向场在该点的方向相一致，基于上述几何解释，从初始点PuUoOo）出发，先依方向场在该点的方向推进到尤=西上—点、Pl,然后再从门依方向场的方向推进到X=X2±—点E，循此前进推岀一条折线
…,—般地，设已做出该折线的顶点几，过几（儿，几）依方向场的方向再推进到几+】（£+],儿+J,显然两个顶点几，几利的坐标有关系
生d=儿）
£+1一&
即
片+1=几+/犷（兀,儿）（4）
这就是著名的欧拉（Euler）公式。
假设初值儿，那么依公式（4）可逐步算出
”=儿+妙'（兀,儿），
例1求解初值问题
（5）
y'=y-2x/y,（0J（0）=h
解欧拉公式的具体形式为
几“二几+如几一2耳Dj
取步长//=0.1,计算结果如下表:
农1计算结果比照
©
咻）
©
心）
0.1
1.1000
1.0954
0.6
1.5090
1.4832
0.2
1.1918
1.1832
0.7
1.5803
1.5492
0.3
1.2774
1.2649
0.8
1.6498
1.6125
0.4
1.3582
1.3416
0.9
1.7178
1.6733
0.5
1.4351
1.4142
1.0
1.7848
1.7321
初值问题（2.1.2）有解y=gx,按这个解析式子算出的准确值y（xn）同近似值几一起列在表X两者相比拟可以看出欧拉方法的精度很差。
（二人向后Euler方法向后Euler方法和Euler方法差不多，只是把y（x/j+1）用y（x“+J-y（x“）
h
去代替，这时计算公式为卜：
了M和％J，（归丄2…）〔6〕
[yM=y^
向后Euler方法的总体截断误差也是一阶的，因此向后Euler方法是收敛的。
这里需要指出它与Euler方法的一个很大不同之处，Euler方法是显式方法，即儿小由兀小，®，儿/明显地表示出来了，而向后Euler方法是隐式方法，计算儿利时要解隐式方程（6）。
通常解此方程用迭代法。
因此计算较为麻烦，但比显式Euler方法精度要咼。
（三）、於法
将Euler方法公式与向后Euler方法公式作加权平均，得到如下公式：
）和=儿+力［&/（兀讥）+（1-&）/（畑,畑）］,…、「力
（〃=0丄2…）I/」
丿（兀）=弘
称式⑺为初值问题
（1）的&-法公式。
其中如y（x°）=儿为初值条件。
在此法中当0=^时，即
儿+】=片+£［/（x”，y”）+/g，y”+J］⑻
此时的法称为梯形公式法。
梯形公式也隐式格式，用起来要进展迭代，其计算公式为
X”，儿）+/（£+］，
这里在应用本迭代法时，是先用Euler方法求初值几+】的近似值儿+F）即：
几+$"=儿+妙（9,几）
然后将儿+F）替代梯形公式（8）中的儿+】得到的式（9）。
式（9）又称为预测校正公式。
换言之，由Euler方法给出预测值，再用梯形法予以校正。
很显然，当步长办取得适当小时，由Euler方法算岀的值已是较好的近似。
格式（9）收敛很快，通常只需一两次迭代即可满足精度要求，假设需屡次迭代,那么应缩小步长力后再行计算。
梯形公式法比用向后Euler方法的迭代步长可以放宽一倍，它的总体截断误差为O（胪），比Euler方法高一阶。
但它每积分一步要计算二次函数值，这说明的精度的提高是以增加计算量为代价的。
（四人改艮欧拉法
为得到比欧拉法精度高的计算公式，在等式如果对方程⑴从兀到心积分，得y（暫+J=yM+广M，y⑴W（1o）
右端积分中假设用梯形求积公式近似，并用儿代替y（xj,儿/弋替y（A^）,那么得
儿+】=y”+£[/（x”，儿）+/g，y”J]⑴）
称为改良欧拉法.
改良欧拉方法是隐式单步法，可用迭代法求解.用欧拉方法提供迭代初值，那么改良欧拉法的迭代公式为
或上儿十幼（S几）；
（12）
<或严=儿+£[/（忑必）+了氏+L，闌）]（—0,12…）.
为了分析迭代过程的收敛性，将⑴）式与（10）相减，得
儿+1-必T）=*1/（心+1,儿+1）-于（心+1,期）]
于是有
式中厶为/（X」）对y满足Lipschitz常数,如果选取h充分小,使得hL“
——<1
2
那么当k-g时有弗T畑，这说明迭代过程（12）是收敛的•例2用改良的欧拉方法求解初值问题
（1）・解改良的欧拉公式为
必?
1二几+幼（S几）；
■或严=儿+级（心必）+/（%闌）.
（—0,12…）.
仍取力=0.1,计算结果见下表.同例1中欧拉法的计算结果比拟，改良欧拉法明显改善了精度.
表2计算结果比照
©
心）
©
心）
0.1
1.0959
1.0954
0.6
1.4860
1.4832
0.2
1.1841
1.1832
0.7
1.5525
1.5492
0.3
1.2662
1.2649
0.8
1.6153
1.6125
0.4
1.3434
1.3416
0.9
1.6782
1.6733
0.5
1.4164
1.4142
1.0
1.7379
1.7321
（五）Runge—Kutta方法
德国数学家C.D.T.Runge（1856——1927）是数值方法开展史上具有里程碑作用的人物。
1895年，他在Hanover发表了关于微分方程数值解法的经典论文?
常微分方程数值解法?
。
此文成为常微分方程Runge-Kutta方法的发端。
此后，Runge结合教学活动积极投身于开展一般的数值分析特别是各种实际应用中的Runge—Kutta方法（严格来说，此方法在Kutta作岀工作后才能称作Runge—Kutta方法）。
Runge-Kutta方法是一种特殊的单步方法,事实上，这个方法可以看作在（x„_,xw+1）上取假设干条积分曲线的假设干个点的切线斜率，再进展一次（或屡次）算术（或加权）平均后产生的新斜率，再按这个斜率从（心，儿）出发，以直线带曲线向前推进一步的过程。
与Taylor展示法相比，Runge—Kutta方法不用增加微商/（俎刃的次数就可以得到较高的阶oRunge—Kutta方法除了在微分方程求解中扮演的传统角色外，人们发现相关类型的初值问题可以用Runge—Kutta方法或适合更一般问题的Runge—Kutta方法求解，比方Runge—Kutta方法被应用到了Hamilton系统中。
前面提到的几种数值解法的精度是很低的，下面给出高阶一步法——Runge-
Kutta方法。
它是最常用的一种数值解法，因为它相当准确、稳定、容易编程。
Runge—Kutta方法至今仍然得到广泛地应用。
1、二级二阶Runge—Kutta方法
由Runge—Kutta方法的思想，我们得到二阶Runge—Kutta公式为
儿+产儿+扌（心+3瓦），
<儿），
22
Ki=hf耳+刃h儿+厅心・
133丿〔13〕
>；+1=儿+£（«|+心），*K严财（兀，儿），笛=炒（£+九儿+«）.
（14）
Runge-Kutta公式是在计算两次函数值的情况下，局部截断误差的阶最高是3,式（13）是允许函数f（x,y）任意变化情况下截断误差最小的二阶方法。
要再提高阶就必须增加计算函数值的次数。
上述式（14）又称为欧拉预估——校正公式。
2、三级三阶Runge—Kuuta方法两个常用的三阶Runge—Kuuta方法分别为:
儿+产儿+討（2匕+3氐+4岛）,
K严.©，儿），
儿+严儿+》（（+%+KJ,
O
K]（兀，儿），
心=fx,t+尹儿+2必J，
K3=f（xn+h,yn-hKl+2hK2）.〔⑹
这两个常用方法在解决实际问题中能够到达较低的精度要求。
但是要更高精度要求的，我们必须了解更高阶的方法一四级四阶Runge-Kuuta方法。
3、四级四阶Runge—Kutta方法
这种方法在解决实际问题中常用。
在这里，我们将公式进展导岀，以熟悉其方法的实用过程。
由于Runge—Kutta公式对初值问题⑴中的一般/（俎刃都适用，那么它必然对特殊的/（x,y）=xy,f（x,y）=x或/（x,y）=y也适用。
从而可以定岀特定的参数来，因此，这里采用的是一种待定系数法，来导出四级四阶显式Runge-Kutta公式，过程十分简单明了。
为了计算简单，令
y（&）近似值为儿，这样求解问题⑴的四阶显示Runge—Kutta公式为儿+1=儿+〃（冰|+“心+C3K3+C4K4）,
&=/（血,儿）,
=f（叫+1叫,儿+/%KJ,K3=/*（£+加3'儿+肋31心+地3K2）,
£=/（召+屁4,儿+肋41&+肪43心）,〔［7〕
这里要求
〔18〕
=*户6=$]+仅2,4=Ql+22+“43・
由于Runge—Kutta公式的思想与数值积分
f1f（x,y）dx人八丿〔19〕
相似。
所以希望四阶显示Runge—Kutta公式对y'=x‘，y（o）=o〔20〕
是准确的。
这样，将〔17〕式应用于〔19〕式，就有
K",
K2=（x+a2）3
=x/+3x：
Clr+3x«72+d,,
K3=（XM+a3y=兀；+3说+3xna^+君,
=（xn+c/J=x„3+3x”咕4+3xna^+dj,
儿+】=儿+（C]+c2+C3+q）y»+（C02+eg+c4a4）y"h2
+c.a3-+c4a43）yn（4'h\
++C3«32+c46/42）+i（c26/23
与Taylor公式
〔21〕
=扣2石+c3a/+c4a4-=1,
儿+i=儿+儿％++y：
h?
+iy：
f+右儿⑷胪
相比拟，有
Cx+
2+6+C
4=1,叱+*+%4=-
C2a2+°3^3+C4a4
又由于Runge—Kutta公式（17）中既要反映x的变化，又要反映y的变化，所以简单的选取/（x,y）为一对线性函数，这样将（13）式用于
y=D，y（O）=O〔22〕
得
儿
K2=K^+a^h+cqKJr
K严K\+3；/+（必32必“+G%）/,+（居仇2心兀+。
2%32丈"，+心3仇2心"，
忑=&+a4ynh+[（a2b42+a3b4i）yxn+a^Kx\lf+
[（alb42K{+a2bi2b42yflxfl+如盲KJ兀+a4（a2bA2+a3b4i）yji\〔23〕
儿+严儿+h（cK+c2K2+C3K3+C4KJ
=儿（q+c2+c3+q）儿+（c2a2+c3a3+c4a4）yh2
+[（qd；+c3af+c4a；）yn+（c3a2by2+c4a2b42+c4^2Z?
43）x'xJ/?
3
+{[56。
032+c4a4（a2b42+a3b43）]yn+
（C/詁32+C/；$2詛43）几亠+。
4心亠2曲
对〔22〕式求导数有
yn=2儿+呂儿，儿=3儿+召儿
将上两式代入（23）式并与Taylor公式（21）比拟有：
即2伏2+4也2+。
讪3=右，即2也2+56（也2+。
血）=右,
5。
2‘毎2+5。
2%42+46%相=占,C02仪2^3=£，
q22
c2a^+c3a3+c4a4
于是总结有如下四阶Runge—Kutta-般显格式中13个参数所满足的11个参数方程
①=2"
h=S\+@2,4=g+“42+“43，
c}+c2+C
3+c4=l9c2a2+c3a3+c4a4=i9
=|，G右+6冶+c/J=土,
勺心2+C4a2h42+C4a3h43=舟，6“2“0324-54（^2^42+^43）=|*彳色‘纽？
+c4a22b42+c4afb43=,c4a2b32b43=右，卫厶〕于是经典的四阶显示Runge-Kutta公式为
儿+严儿+：
（心+2心+2心+心）,O
K严妙（心,儿），
K2=,Jf兀+彳,儿+*
升级会员 