常微分方程初值问题数值解法1
常微分方程的初值问题及其解法
常微分方程的初值问题及其解法常微分方程是自然界中各种变化的基础模型，广泛应用于物理、工程、生物、经济学等领域。
初值问题是其中最基本的问题之一。
本文将从初值问题的意义入手，介绍几种不同的数值解法，并评价其优缺点。
1. 初值问题的意义首先，我们来看一个简单的例子。
假设有一个人从一楼的窗户往下跳，忽略空气阻力，我们可以列出他下落的物理规律：$$\frac{d^2h}{dt^2}=g$$其中$h$是跳下来后距离地面的高度，$t$是时间，$g$是常数，表示重力加速度。
上面这条式子就是一个二阶常微分方程。
我们的问题是，如果知道了他的初速度$v_0$和起始高度$h_0$，能否求得他下落到地面时的时间和高度。
这个例子中，$h$和$t$都是连续的量，但是我们并不能解析地求出$h(t)$的解析式，因此需要用数值方法去近似求解。
这就是初值问题的意义。
通常，初值问题是指某一初始时刻$t_0$的初值：$$y'(t_0)=f(y(t_0),t_0),\ y(t_0)=y_0$$其中$y$是未知函数，而$f$则是已知函数。
对于一阶常微分方程，这个条件是充分的，可以唯一地决定一个解。
但是对于更高阶的常微分方程，则需要多个初始条件才能确定一个解。
然而，这已经超出了本文的范畴，这里只讨论一阶常微分方程的初值问题。
2. 数值解法下面将介绍几种常见的数值解法。
2.1. 欧拉法欧拉法是最简单的数值解法之一，其思路是将初值问题离散化。
具体来说，我们可以将时间$t$分成若干个小段，每段的长度为$\Delta t$。
于是，我们可以将初始时刻$t_0$的初始值$y(t_0)=y_0$，并通过欧拉法近似计算下一个时间点$t_0+\Delta t$的值$y_1$：$$y_1=y_0+f(y_0,t_0)\Delta t$$同理，我们可以通过已知的$y_1$和$t_1=t_0+\Delta t$，计算下一个时间点$t_2=t_0+2\Delta t$的值$y_2$：$$y_2=y_1+f(y_1,t_1)\Delta t$$依此类推，直到我们得到一个目标时间$t_m$的值$y_m$。
浅谈常微分方程初值问题数值解法
浅谈常微分方程初值问题数值解法在自然科学、工程技术、甚至社会科学的一些领域中，常常会遇见一阶常微分方程的求解问题：()上述问题，寻求解的具体表达式十分困难，仅对一些特殊形式的才有可能找到解的解析表达式，在大多情况下，初值问题的解不能用初等函数表示出来即使可写出解的解析表达式，但因为这些表达式过于复杂，要计算它在某些点上的函数值也异常困难。
在实际问题中，经常需要的恰是解在某些点上的函数值，因此研究初值问题的数值解法十分必要。
1 常微分方程初值问题的数值解法常微分方程的近似解法大体可分成三大类：一类是图解法和器械法；第二类是解的近似法；第三类是数值解法，即通过离散化的方法直接求出函数在某些点上的近似值，此数值解仅为精确解的近似解。
其基本原理为：一阶常微分方程的初值问题的解是上变量的连续函数，因此求上述问题的数值解，就是在区间上的若干离散点上用离散化的方法将初值问题化成离散变量的相应问题，从而相应问题的解可作为初值问题理论解的近似值。
由常微分方程的理论可知，只要在区域内连续，且关于满足林普希兹条件，则方程的解存在且唯一。
初值问题的数值解法通常采取“步进法”，而“步进法”又可分为“单步法”和“多步法”两类。
（1）单步法。
所谓“单步法”是指在计算时，只用到前一步的有关信息。
其一般形式为：，主要包括下面三种方法：Euler方法，改进的Euler公式-梯形公式和Runge-Kutta法。
（2）线性多步法。
单步法没有用到前几步计算得到的信息，因此为了提高精度，需重新计算多个点处的函数数值，如RK方法，故计算量较大。
线性多步法的基本思想是充分利用前面的已知信息来构造精度高且计算量小的算法来计算。
多步法常用方法是线性多步法，求解公式为：构造的常用方法是Taylor展开和数值积分方法。
常用的线性多步公式有：四阶Adams显式公式：四阶Adams隐式公式：四阶Milne显式公式：三阶Hamming公式：（隐式公式）预测校正系统和预测校正修正法：一般地，同阶的隐式法比显式法精确，而且数值稳定性好，但隐式公式中的求解较难，需要用到迭代法，这就增加了计算量。
常微分方程初值问题数值解法
数值解法的必要性
实际应用需求
许多实际问题需要求解常微分方程初值问题，如物理、 化学、生物、工程等领域。
解析解的局限性
对于复杂问题，解析解难以求得或不存在，因此需要 采用数值方法近似求解。
数值解法的优势
未来发展的方向与挑战
高精度算法
研究和发展更高精度的算法，以提高数值解的准确性和稳定性。
并行计算
利用并行计算技术，提高计算效率，处理大规模问题。
自适应方法
研究自适应算法，根据问题特性自动调整计算精度和步长。
计算机技术的发展对数值解法的影响
1 2
硬件升级
计算机硬件的升级为数值解法提供了更强大的计 算能力。
它首先使用预估方法（如欧拉方法）得到一个 初步解，然后使用校正方法（如龙格-库塔方法） 对初步解进行修正，以提高精度。
预估校正方法的优点是精度较高，且计算量相 对较小，适用于各种复杂问题。
步长与误差控制
01
在离散化过程中，步长是一个重要的参数，它决定 了离散化的精度和计算量。
02
误差控制是数值逼近的一个重要环节，它通过设定 误差阈值来控制计算的精度和稳定性。
能够给出近似解的近似值，方便快捷，适用范围广。
数值解法的历史与发展
早期发展
早在17世纪，科学家就开始尝 试用数值方法求解常微分方程。
重要进展
随着计算机技术的发展，数值 解法在20世纪取得了重要进展， 如欧拉法、龙格-库塔法等。
当前研究热点
目前，常微分方程初值问题的 数值解法仍有许多研究热点和 挑战，如高精度算法、并行计
软件优化
软件技术的发展为数值解法提供了更多的优化手 段和工具。
微分方程初值问题的数值解法
积分法:
yk 1 yk h f ( xk , yk ) y ( x0 ) y0
积分项利用矩形公式计算
(1) y( xk 1 ) y( xk )
xk 1
xk
f (t , y(t ))dt
(★)
xk 1
xk
f (t , y(t ))dt h f ( xk , yk ) y( xk 1 ) y( xk ) h f ( xk , yk )
引言
初值问题的数值解法:求初值问题的解在一系列节点的值 y ( xn )的近似值 yn 的方法.本章数值解法的特点:都是采用“步进 式”,即求解过程顺着节点排列的次序一步步向前推进. 常微分方程初值问题: dy f ( x, y ), x [a, b] dx y ( x0 ) y0
替 f (x , y)关于 y 满足Lipschitz条件. 除了要保证(1)有唯一解外,还需保证微分方程本身是稳定的,即 (1)的解连续依赖于初始值和函数 f (x , y). 也就是说, 当初始值 y0 及函数 f (x , y)有微小变化时, 只能引起解的微小变化.
注: 如无特别说明,总假设(1)的解存在唯一且足够光滑. 在 f 连续有界, 则 f (x , y)对变量 y 可微的情形下, 若偏导数 y 可取L为
也称折线法 x
2. 梯形法
若采用梯形公式计算(★)中的积分项,则有 h y ( xk 1 ) y ( xk ) [ f ( xk , y ( xk )) f ( xk 1 , y ( xk 1 ))] 2 h yk 1 yk [ f ( xk , yk ) f ( xk 1 , yk 1 )] 2 称之为梯形公式.这是一个隐式公式,通常用迭代法求解.具体做 法: (0) (0) 先用Euler法求出初值 yk ,1 即 ,将其代入梯形公式 yk 1 yk h f ( xk , yk ) 的右端,使之转化为显式公式,即 h ( l 1) (l ) yk 1 yk [ f ( xk , yk ) f ( xk 1 , yk (☆ ) 1 )] 2
第9章 常微分方程初值问题数值解法
数值分析
第9章 常微分方程初值问题数值解法
《常微分方程》中介绍的微分方程主要有:
(1)变量可分离的方程 (2)一阶线性微分方程(贝努利方程) (3)可降阶的一类高阶方程 (4)二阶常系数齐次微分方程 (5)二阶常系数非齐次微分方程 (6)全微分方程 本章主要介绍一阶常微分方程初值问题的数值解法。
进一步: 令
y n1 y n
xn 1 xn
y n 1 y( x n 1 ) , y n y( x n )
f ( x , y( x ))dx h f ( x n , y n )
宽
9
高
实际上是矩形法
数值分析
第9章 常微分方程初值问题数值解法
(3)
用Taylor多项式近似并可估计误差
解决方法：有的可化为显格式，但有的不行 18
数值分析
第9章 常微分方程初值问题数值解法
与Euler法结合，形成迭代算法 ，对n 0，2， 1，
( yn0 )1 yn hf x n , yn ( k 1) h ( yn1 yn f x n , yn f x n1 , ynk )1 2
7
数值分析
第9章 常微分方程初值问题数值解法
建立数值解法的常用方法
建立微分方程数值解法,首先要将微分方程离散 化. 一般采用以下几种方法: (1) 用差商近似导数
dy yx yx x x dx x y
n 1 n n 1 n
n
,
n
进一步: 令
yn1 y( xn1 ) , yn y( xn )
由 x0 , y0 出发取解曲线 y y x 的切线（存在！），则斜率
常微分方程初值问题数值解法
0.4 1.3582 1.3416 0.9 1.7178 1.6733
0.5 1.4351 1.4142 1.0 1.7848 1.7321
7
初值问题（2.2）有解 y ，1按2这x 个解析式子
算出的准确值 y(x同n )近似值 一y起n 列在表9-1中，两者 相比较可以看出欧拉方法的精度很差.
17
所以，局部截断误差可理解为用方法（2.10）计算一步的 误差，也即公式（2.10）中用准确解y(x代) 替数值解产生
的公式误差.
根据定义，显然欧拉法的局部截断误差
Tn1 y( xn1) y( xn ) hf ( xn , y( xn ))
y(xn h) y(xn ) hy(xn )
y(2) n1
yn
hf
( xn1,
y (1) n1
).
11
如此反复进行，得
y (k 1) n1
yn
hf
( xn1,
y(k) n1
),
(k 0,1, ).
（2.6）
由于 f (x,对y) 满足y 利普希茨条件（1.3）. 由（2.6）减 （2.5）得
y (k 1) n 1
yn1
h
f
( xn1,
y(k) n 1
积分曲线上一点 (x的, y切)线斜率等于函数 值.
的f (x, y)
如果按函数 f (在x, y) 平x面y上建立一个方向场，那 么，积分曲线上每一点的切线方向均与方向场在该点的方 向相一致.
基于上述几何解释，从初始点 P0 (x出0 ,发y0，) 先依 方向场在该点的方向推进到 x 上x1一点 ，P然1 后再从 P1 依方向场的方向推进到 x 上x2一点 ，循P2此前进做出
常微分方程初值问题的数值解法
常微分方程初值问题数值解法初值问题：即满足初值条件的常微分方程的解y′=f(x,y),x∈[x0,b]y(x0)=y0.定理1（利普希茨条件）若存在正数L，使得对任意，y1，y2，有|f(x,y1)−f(x,y2)|≤L|(y1−y2)|定理2（解存在性）①若函数f在方区域x∈[a,b],y∈R连续，②函数f关于y 满足利普希茨条件，则对任意x∈[a,b]，常微分方程存在唯一的连续可微数值解.两类问题：①单步法---计算下一个点的值yn+1只需要用到前面一个点的值yn②多步法---计算下一个点的值yn+1需要用到前面l个点的值yl1、欧拉法---下一个点的计算值等于前一个点的计算值加上步长乘以前一个点的函数值•具体过程一些批注：显式欧拉方程指下一步要计算的值，不在迭代方程中；隐式欧拉方程指下一步要计算的值，在迭代方程中。
怎么计算隐式欧拉方程----要借助显示欧拉迭代计算---一般用迭代法-----迭代---将微分方程在区间[xn,xn+1]进行积分，然后函数f进行近似，即可得到迭代方程-----迭代方程收敛性？由函数关于y满足利普希茨条件，可以推出迭代公式收敛。
•局部截断误差：假设前n步误差为0，我们计算第n+1步的误差，将次误差称为局部截断误差，且局部误差为O(hp+1)•p阶精度：由理论证明：若局部误差阶的时间复杂度为O(hp+1)，则整体误差阶为O(hp)我们称公式精度为p。
•显示欧拉法与隐式欧拉法•梯形方法----将显式欧拉迭代方程与隐式欧拉迭代方程做一下加权平均，构造的计算公式.•改进的欧拉方法---思想：因为梯形公式是隐式公式，将显式欧拉公式对下一步的计算值进行预估，用梯形公式对下一步的计算值进行校正.2、龙格-库塔方法思想：根据Lagrange中值定理，下一次的计算值可以用前一次的计算值加上h乘以前一个点的斜率；而这个斜率用该区间上的多个点的斜率的算数平均来逼近。
注意：怎么计算任意斜率Ki？第i个点的斜率Ki有微分方程可以算出f′=f(xn,yn)所以要算的f(xn,yn)值，由欧拉法即可算出, yn+1=yn+hf′•2阶-龙格-库塔方法----类似改进的欧拉法根据Lagrange中值定理，下一次的计算值可以用前一次的计算值加上h乘以斜率；而这个斜率用区间上的端点和中点的斜率的算数平均来逼近。
一阶常微分方程初值问题的数值解fortron
一阶常微分方程初值问题的数值解是微分方程数值求解中的基础问题，对于工程、物理、生物等领域的科学计算和数值模拟具有重要意义。
本文将从常微分方程初值问题的数值解的基本原理和数值方法入手，详细介绍使用Fortran语言进行一阶常微分方程初值问题的数值解的实现过程和注意事项。
1. 常微分方程初值问题的数值解基本原理常微分方程初值问题的数值解是通过数值方法来逼近微分方程的解。
对于一阶常微分方程初值问题：dy/dx = f(x, y)y(x0) = y0其中，f(x, y)是给定的函数，y(x0) = y0是给定的初值条件。
求解该初值问题即是要找到一个函数y(x)近似地满足该微分方程，并且在点x = x0处与给定的初值条件相符。
2. 常微分方程初值问题的数值解的数值方法常见的数值方法包括欧拉方法、改进的欧拉方法、四阶龙格-库塔方法等。
其中，四阶龙格-库塔方法是最常用和最经典的数值方法之一。
该方法通过取若干个函数值点上的斜率的加权平均值来逼近微分方程的解，具有较高的数值精度和稳定性。
3. 使用Fortran语言实现一阶常微分方程初值问题的数值解Fortran是一种古老但经典的科学计算语言，以其高效的数值计算和科学工程计算而闻名。
下面将结合Fortran语言的特点，介绍如何使用Fortran语言实现一阶常微分方程初值问题的数值解。
（1）定义常微分方程的函数f(x, y)在Fortran程序中，首先需要定义常微分方程的函数f(x, y)，并将其定义为一个子程序或函数。
这里以一个简单的一阶线性常微分方程为例：f(x, y) = x + y则在Fortran程序中可以这样定义：```function f(x, y)real :: x, y, ff = x + yend function f```（2）实现四阶龙格-库塔方法在Fortran程序中，可以实现四阶龙格-库塔方法来数值解常微分方程初值问题。
具体做法是按照龙格-库塔方法的算法，在程序中编写相应的代码实现。
数值分析李庆扬第9章常微分方程初值问题数值解法讲义.
② 由 x0 , y0 f x0 , y0 切线 P0P1 ，
切线与 x x1 交点 P1 ： y1 的近似值 ；
③ 再由 x1 , y1 向前推进到 P2 ， 得到折线 P0P1 Pn ，近似 y yx 。
7
2021年5月4日
《数值分析》 黄龙主讲
h
yxn
yxn1
yn1 yn f
h
xn1 , yn1
yn1 yn h f xn1 , yn1
——后退的欧拉公式（隐式）
注意：① 显式计算方便，隐式稳定性较好；
② 上式隐含 yn1 ? ，采用迭代法求解。
12
2021年5月4日
《数值分析》 黄龙主讲来自欧拉公式的另一种理解：
将常微分方程 y f x, y 改写 dy f t , ytdt
“步进式”：顺着节点排列顺序，一步一步地向前推进。
步长：常用等步长 hn xn1 xn ，节点为 xn x0 nh 单步法：计算 yn1 时，只用到前一点的值 yn k 步法：计算 yn1 时，用到前面 k 点的值 yn , yn1 , , ynk1
5
2021年5月4日
《数值分析》 黄龙主讲
对微分方程从 xn 到 xn1 积分
y xn1 yxn
xn1 f t , yt dt
xn
由积分左矩形公式得
xn1 xn
f
t ,
yt dt
hf
xn ,
yxn
例如：
lim
h0
yxn1
h
yxn
yxn
yxn1
h
yxn
yxn
f xn , yxn
计算方法课件第八章常微分方程初值问题的数值解法
整体截断误差与局部截断误差的关系
定理：如果f(x，y)满足李普希兹（Lipschitz）条件
f(x ,y 1 )f(x ,y 2) L y 1y 2
且局部截断误差有界：
|R n|1 2h2M 2
(n1,2, )
则Euler法的整体截断误差n满足估计式：
ne(ba)L 0h 2L M 2(e(ba)L1)
分光滑。初值问题的解析解(理论解）用 y(x表n ) 示, 数值解法的精确解用 y表n 示。
常微分方程数值解法一般分为：
（1）一步法：在计算y n 1 时，只用到x n 1 ,x n和 y，n 即前一步的值。
（2）多步法：计算 y n 1 时，除用到 x n 1 ,x n 和 y n 以外，还要用 x n p 和 y n p (p1 ,2 k;k0) ,即前
其中L为李普希兹常数，b-a为求解区间长度，
M2 mayx(x) 。 axb
证明参见教材。
Remark：该定理表明，整体截断误差比局部截 断误差低一阶。对其它方法，也有类似的结论。
收敛性与稳定性
收敛性定义：如果某一数值方法对于任意固定的
xn=x0+nh，当h0(同时n )时有yn y(xn)，
则称该方法收敛。 稳定性定义 定义 用一个数值方法，求解微分方程初值问 题时，对给定的步长h>0，若在计算 y n 时引入 误差 （n 也称扰动），但由此引起计算后面的 ynk(k1,2, )时的误差按绝对值均不增加，则 称这个数值方法是稳定的。
一般的显式rk方法可以写成型钢截面只需少量加工即可用作构件省工省时成本低但型钢截面受型钢种类及型钢号限制难于完全与受力所需的面积相对应用料较多其中为常数选取这些常数的原则是要求第一式的右端在处泰勒展开后按h型钢截面只需少量加工即可用作构件省工省时成本低但型钢截面受型钢种类及型钢号限制难于完全与受力所需的面积相对应用料较多上述公式叫做n级的rungekutta方法其局部截断误差为显然euler法是一级一阶rk方法
第一章常微分方程初值问题数值解法
（1.2.3）
其中rn,k(t)为插值余项。 代到（1.2.2）式中得
u ( tn +1 ) = u ( tn ) +
舍去余项 并用uj代替u(tj)即得
∫
tn+1 tn
Ln , k ( t ) dt + ∫ t rn , k ( t ) dt
tn+1
n
(1.2.4 (1.2.5)
Rn , k = ∫
⎡ ∑ ⎣α u
j =0
j =0
αk ≠ 0
（1.2.1）
j n+ j
⎤ − hβ j f n + j ⎦ = 0（数值解满足的差分方程）
因此称（1.2.1）为多步法 或 k-步法。 又因为(1.2.1)关于 u n + j , f n + j 是线性的，所以称为线性多步法。 为使多步法的计算能够进行，除给定的初值u0 外，还要 知道附加初值u1,u2,…,uk-1 ，这可用其它方法计算。 若 β k = 0 则称(1.2.1)是显式的； 若 β k ≠ 0 则方法（1.2.1）是隐式的。 例如，一般线性二步法可写成：
f ( t , u ( t ) ) = Ln , k +1 ( t ) + rn , k +1 ( t )
其中rn,k+1(t)为插值余项。 同理即
un +1 = un + h ∑ bk +1i f ( tn −i +1 , un −i +1 )
i =0
k +1
其中
bk +1i
=∫ ∏
−1
j =0 j ≠i
0
k +1
9、常微分方程初值问题数值解法
( k +1) yn +1
xn +1 ∫xn
− yn +1 |≤
hL 2
|
(k ) yn +1
− yn +1 |,
( 只要 hL < 1，则( 2.8)的ynk +1)收敛到(2.7)的yn +1. +1 2
三、单步法的局部截断误差与阶
一阶常微分方程初值问题(1.1)(1.2)的单步法的一般形式 yn +1 = yn + hϕ ( xn , yn , yn +1, h).
clear x=0,yn=1 %初始化 for n=1:10 yp=yn+0.1*(yn-2*x/yn); %预测 x=x+0.1; yc=yn+0.1*(yp-2*x/yp) ; yn=(yp+yc)/2 %校正 end
( 2.2)
作业： 作业：P381, 1, 2(1).
龙格—库塔 库塔(Runge-Kutta)法 §3 龙格 库塔 法
进一步 y ( xn +1 ) = y ( xn ) + ∫
xn +1 xn
f ( x, y ( x))dx,
(Байду номын сангаас.3)
∫
⇒ 其中
xn +1 xn
f ( x, y ( x))dx ≈ h ∑ ci f ( xn + λi h, y ( xn + λi h)).
yn +1 = yn + hϕ ( xn , yn , h),
i =1
r
(3.4) (3.5) 欧拉法r = 1, p = 1.改进
计算方法 常微分方程初值问题数值解法 Euler公式 龙格 库塔法
第12次 常微分方程初值问题数值解法
内容
1. 常微分方程初值问题解的存在性定理 2. Euler公式 3. 梯形公式 4. 两步Euler公式 5. 欧拉法的局部截断误差 6. 改进型Euler公式 7. 龙格-库塔法 8. 算法实现
常微分方程初值问题 解的存在性定理
欧拉法的局部截断误差
9.2.4． 欧拉法的局部截断误差
衡量求解公式好坏的一个主要标准是求解公式的精 度, 因此引入局部截断误差和阶数的概念。
定义9.1 在yi准确的前提下, 即 yi y(xi)时, 用数值
方法计算yi+1的误差：
R i yi( 1 ) x y i 1
称为该数值方法计算时yi+1的局部截断误差。
通常表示成下列平均化形式
yi1的近似 y p y i hf(x i , y i )
yi1的近似 y c y i hf(x i 1 , y p )
yi1的更好 的近似
y i1
1 2
(y p
yc)
i 0, 1, 2 … , n 1
(9.12)
例9.2 用改进欧拉法解初值问题
y
y
2x y
f( y 1 ) x f,( y 2 ) x L y , 1 y 2 , y 1 ,y 2 R
则方程( 9.1 ) 在a, b上存在唯一的连续可微分的 解的解 y=y(x) 。
推论：如果函数f(x,y)对y的偏导数
f(x, y
y)
在带形
区域 R { a x b ,y - }
内有界。
两步欧拉公式的局部截断误差为：
y(i1 x )yi1O3 ()h
从而两步欧拉公式的阶数是2.推导过程省略。
常微分方程初值问题的数值解法
常微分方程初值问题的数值解法在实际应用中，对于某些微分方程，我们并不能直接给出其解析解，需要通过数值方法来求得其近似解，以便更好地理解和掌握现象的本质。
常微分方程初值问题（IVP）即为一种最常见的微分方程求解问题，其求解方法有多种，本文将对常微分方程初值问题的数值解法进行较为详细的介绍。
一、欧拉法欧拉法是最基本的一种数值解法，它采用泰勒级数展开并截断低阶项，从而获得一个差分方程近似求解。
具体来讲，设 t 为独立变量，y(t) 为函数 y 关于 t 的函数，方程为：$$y'(t) = f(t, y(t)), \qquad y(t_0) = y_0$$其中 f(t,y(t)) 为已知的函数，y(t_0) 为已知的初值。
将函数 y(t) 进行泰勒级数展开：$$y(t+h) = y(t) + hf(t, y(t)) + O(h^2)$$其中 h 表示步长，O(h^2) 表示其他高阶项。
为了使误差较小，一般取步长 h 尽可能小，于是我们可以用欧拉公式表示数值解：$$y_{n+1} = y_n + hf(t_n, y_n), \qquad y_0 = y(t_0)$$欧拉法的优点是容易理解和实现，但是由于截取低阶项且使用的单步法，所以误差较大，精度较低，在具体应用时需要慎重考虑。
二、龙格-库塔法龙格-库塔法（Runge-Kutta method）是一种多步法，比欧拉法更加精确。
龙格-库塔法的主要思想是使用不同的插值多项式来计算近似解，并且将时间步长分解，每次计算需要多次求解。
以下简要介绍二阶和四阶龙格-库塔法。
二阶龙格-库塔法将时间步长 h 分解成两步 h/2，得到近似解表达式：$$\begin{aligned} k_1 &= hf(t_n, y_n)\\ k_2 &= hf(t_n+h/2,y_n+k_1/2)\\ y_{n+1} &= y_n+k_2+O(h^3)\\ \end{aligned}$$四阶龙格-库塔法四阶龙格-库塔法是龙格-库塔法中应用最为广泛的一种方法，其需要计算的中间值较多，但是具有更高的精度。
【论文】常微分方程初值问题数值解法的讨论与比较
【关键字】论文《微分方程数值解法》论文常微分方程初值问题数值解法的讨论与比较一、一阶常微分方程的初值问题科学技术中常常要求解常微分方程的定解问题，这类问题最简单的形式是一阶方程的初值问题我们知道，只要函数适当光滑，譬如关于满足Lipschitz条件（1.3）理论上就可以保证初值问题(1.1),(1.2)的解存在并且唯一。
虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特殊类型的方程，实际问题中归结出来的微分方程主要靠数值解法。
2、欧拉法我们知道，在平面上，微分方程(1.1)的解称作它的积分曲线。
积分曲线上一点的切线斜率等于函数的值，如果按函数在平面上建立一个方向场，那么，积分曲线上每一点的切线方向均与方向场在该点的方向相一致，基于上述几何解释，从初始点出发，先依方向场在该点的方向推进到上一点，然后再从依方向场的方向推进到上一点，循此前进推出一条折线，一般地，设已做出该折线的顶点，过依方向场的方向再推进到，显然两个顶点，的坐标有关系即(2.1)这就是著名的欧拉(Euler)公式。
若初值已知，则依公式(2.1)可逐步算出例1 求解初值问题解欧拉公式的具体形式为取步长，计算结果如下表：表1 计算结果对比初值问题(2.2)有解，按这个解析式子算出的准确值同近似值一起列在表1，两者相比较可以看出欧拉方法的精度很差。
三、改进欧拉法为得到比欧拉法精度高的计算公式，在等式如果对方程(1.1)从到积分,得(3.1)右端积分中若用梯形求积公式近似，并用代替，代替，则得(3.2)称为改进欧拉法．改进欧拉方法是隐式单步法，可用迭代法求解．用欧拉方法提供迭代初值，则改进欧拉法的迭代公式为(3.3)为了分析迭代过程的收敛性，将(3.1)式与(3.2)相减，得，于是有，式中为对满足Lipschitz常数，如果选取充分小，使得，则当时有，这说明迭代过程(3.3)是收敛的．例2用改进的欧拉方法求解初值问题(1.1)．解改进的欧拉公式为仍取，计算结果见下表．同例１中欧拉法的计算结果比较，改进欧拉法明显改善了精度．表2 计算结果对比四、线性多步法在逐步推进的求解过程中，计算之前事实上已经求出了一系列的近似值，，如果充分利用前面多步的信息来预测，则可以期望会获得较高的精度．这就是构造所谓线性多步法的基本思想．构造多步法的主要途径是基于数值积分方法和基于泰勒展开方法，前者可直接由方程(1.1)两端积分后利用插值求积公式得到．一般的线性多步法公式可表示为，(4.1)其中为的近似，，，，为常数，及不全为零，则称(4.1)为线性步法，计算时需先给出前面个近似值，再由（4.1)逐次求出．如果，称(4.1)为显式步尖，这时可直接由(4.1)算出；如果，则(4.1)称为隐式步法，求解时与改进欧拉法相同，要用迭代法方可算出，(4.1)中系数及可根据方法的局部截断误差及阶确定，其定义为：设是初值问题(1.1),(1.2)的准确解。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
yn+1 = yn + hf ( xn , yn ) n+ y1 = y0 + hf ( x0 , y0 ), y2 = y1 + hf ( x1, y1 ),
(2 .1)
这就是著名的(显式)欧拉(Euler)公式. 若初值y0已 这就是著名的(显式)欧拉( )公式. 若初值 则依公式(2.1)可逐次逐步算出各点数值解. 可逐次逐步算出各点数值解. 知，则依公式 可逐次逐步算出各点数值解
与准确解 y = 1+ 2x相比，可看出欧拉公式的计算结 + 相比， 果精度很差. 果精度很差
上页
下页
欧拉公式具有明显的几何意义, 就是用折线近似 欧拉公式具有明显的几何意义, 就是用折线近似 代替方程的解曲线，因而常称公式(2.1)为欧拉折线法. 代替方程的解曲线，因而常称公式 为欧拉折线法. 还可以通过几何直观来考察欧拉方法的精度.假 还可以通过几何直观来考察欧拉方法的精度. 设yn=y(xn),即顶点 n落在积分曲线 ,即顶点P 落在积分曲线y=y(x)上，那么， 上 那么， 按欧拉方法做出的折线 PnPn+1便是 便是y=y(x)过点 n 过点P 过点 的切线.从图形上看, 的切线.从图形上看,这 样定出的顶点P 样定出的顶点 n+1显著 地偏离了原来的积分曲 可见欧拉方法是相 线，可见欧拉方法是相 当粗糙的 当粗糙的.
可得欧拉法(2.1)的公式误差为 的公式误差为 可得欧拉法
为了分析计算公式的精度， 为了分析计算公式的精度，通常可用泰勒展开 在 处展开， 将y(xn+1)在xn处展开，则有
2
h h y( xn+1 ) − yn+1 = y′′(ξn ) ≈ y′′( xn ), 2 2
称为此方法的局部截断误差. 称为此方法的局部截断误差. 局部截断误差
上页 下页
设用欧拉公式
( yn0)1 = yn + hf ( xn , yn ) + (0) 用它代入(2.5)式的右端，使之转 式的右端 式的右端， n+1 ，用它代入
给出迭代初值 y
化为显式， 化为显式，直接计算得
( ( yn1+)1 = yn + hf ( xn+1, yn0)1 ), +
然后再用 y
2
2
(2.3)
上页
下页
如果对方程(1.1)从xn到xn+1积分，得 从 积分， 如果对方程
y( xn+1 ) = y( xn ) + ∫
xn+1
xn
f (t, y(t ))dt .
(2.4)
右端积分用左矩形公式 近似， 右端积分用左矩形公式hf(xn,y(xn))近似，再以 n代替 左矩形公式 近似 再以y y(xn)，yn+1代替 n+1)也得到欧拉公式 代替y(x 也得到欧拉公式 也得到欧拉公式(2.1)，局部截 ， ， 断误差也是(2.3). 断误差也是 如果右端积分用右矩形公式 右矩形公式hf(xn+1,y(xn+1))近似， 近似， 如果右端积分用右矩形公式 近似 则得到另一个公式
(1) 代入(2.5)式，又有 式 n+1代入 (2) (1) n+1 n n+1 n+1
y = y + hf ( x , y ).
如此反复进行， 如此反复进行，得
( ( ynk+1) = yn + hf ( xn+1 , ynk)1 ) (k = 0,1,L). (2.6) +1 +
上页 下页
由于f(x, y)对y满足 满足Lipschitz条件 条件(1.3). 由(2.6)减(2.5)得 由于 对 满足 条件 减 得
x1 < x2 < L< xn < xn+1 < L
上的近似值 y1,y2,L,yn,yn+1,L. 相邻两个节点的间距 L L hn=xn+1-xn称为步长 今后如不特别说明，总是假定 称为步长 今后如不特别说明， 步长. hi=h(i=1,2,L)为定数 这时节点为 n=x0+nh(i=0,1,2,L) L 为定数, 这时节点为x L (等距节点 等距节点). 等距节点
y′ = f ( x, y), y( x0 ) = y0 .
(1.1) (1.2)
我们知道，只有f(x, y)适当光滑 譬如关于 满足 适当光滑—譬如关于 我们知道，只有 适当光滑 譬如关于y满足 利普希茨(Lipschitz)条件 利普希茨 条件
f ( x, y) − f ( x, y) ≤ L y − y .
理论上就可以保证初值问题的解y＝ 存在并且唯一. 理论上就可以保证初值问题的解 ＝f(x)存在并且唯一 存在并且唯一
上页 下页
虽然求解常微分方程有各种各样的解析方法， 虽然求解常微分方程有各种各样的解析方法，但 解析方法只能用来求解一些特殊类型的方程， 解析方法只能用来求解一些特殊类型的方程，实际问 题中归结出来的微分方程主要靠数值解法. 题中归结出来的微分方程主要靠数值解法 所谓数值解法 就是寻求解y(x)在一系列离散节点 所谓数值解法, 就是寻求解 数值解法 在一系列离散节点
y
pn
pn +1
y = y ( x)
p ( x n +1 )
xn
xn +1
上页
x
下页
h y( xn+1 ) = y( xn + h) = y( xn ) + hy′( xn ) + y′′(ξn ) 2 h2 = y( xn ) + hf ( xn , yn ) + y′′(ξn ) ξn ∈( xn , xn+1 ). 2 的前提下， 在yn=y(xn)的前提下，f(xn,yn )=f(xn,y(xn))=y′(xn).于是 的前提下
上页 下页
9.2 简单的数值方法与基本概念
9.2.1 欧拉法与后退欧拉法 我们知道， 平面上， 我们知道，在xy平面上，微分方程 平面上 微分方程(1.1)式的解 式 y=f(x)称作它的积分曲线，积分曲线上一点 y)的切 称作它的积分曲线 上一点(x, 的切 称作它的积分曲线，积分曲线上一点 线斜率等于函数f(x, y)的值 如果按 的值. 线斜率等于函数 的值 如果按f(x, y)在xy平面上 在 平面上 建立一个方向场，那么，积分曲线上每一点的切线 建立一个方向场，那么，积分曲线上每一点的切线 方向均与方向场在该点的方向相一致. 方向均与方向场在该点的方向相一致 基于上述几何解释，我们从初始点 出发, 基于上述几何解释，我们从初始点P0(x0, y0)出发 出发 上一点P 先依方向场在该点的方向推进到x=x1上一点 1，然后 先依方向场在该点的方向推进到 再从P 再从 1点依方向场在该点的方向推进到 x=x2 上一点 P2 , 循环前进做出一条折线 0 P1 P2L. 循环前进做出一条折线 折线P
2x1 0.2 y2 = y1 + h( y1 − ) = 1.1+ 0.1(1.1− ) = 1.191818 y1 1.1 L L
上页 下页
依次计算下去，部分计算结果见下表. 依次计算下去，部分计算结果见下表 见下表 xn 欧拉公式数值解 n 准确解 n) 欧拉公式数值解y 准确解y(x 0.2 0.4 0.6 0.8 1.0 1.191818 1.358213 1.508966 1.649783 1.784770 1.183216 1.341641 1.483240 1.612452 1.732051 误差 0.008602 0.016572 0.025726 0.037331 0.052719
上页
下页
9.2.2 梯形方法 为得到比欧拉法精度高的计算公式 在等式 欧拉法精度高的计算公式， 为得到比欧拉法精度高的计算公式，在等式(2.4) 右端积分用梯形求积公式近似, 并用y 代替y(x 右端积分用梯形求积公式近似 并用 n代替 n), yn+1 代替y(x ， 代替 n+1)，则得
h yn+1 = yn + [ f ( xn , yn ) + f ( xn+1, yn+1 )] , 2
上页 下页
初值问题的数值解法有个基本特点， 初值问题的数值解法有个基本特点，他们都采 数值解法有个基本特点 取“步进式”，即求解过程顺着节点排列的次序一步 步进式” 一步地向前推进. 描述这类算法， 一步地向前推进 描述这类算法，只要给出用已知信 息yn,yn-1,yn-2,L计算 n+1的递推公式. 的递推公式 - L计算y 首先，要对微分方程离散化， 首先，要对微分方程离散化，建立求解数值解的 递推公式. 一类是计算 n+1时只用到前一点的值 n，称 递推公式 一类是计算y 时只用到前一点的值y 单步法. 另一类是用到y 点的值y 为单步法 另一类是用到 n+1前面 k 点的值 n,yn-1,L, L yn-k+1，称为 步法 其次，要研究公式的局部截断误差 称为k步法 其次，要研究公式的局部截断误差 步法. 数值解y 与精确解y(x 的误差估计及收敛性， 和阶，数值解 n与精确解 n)的误差估计及收敛性， 还有递推公式的计算稳定性等问题. 还有递推公式的计算稳定性等问题 计算稳定性等问题
称为矩形方法 称为矩形方法. 矩形方法
(2.7)
矩形方法是隐式单步法，用迭代法求解， 矩形方法是隐式单步法，用迭代法求解，同后 退的欧拉方法一样，仍用欧拉法提供迭代初值， 退的欧拉方法一样，仍用欧拉法提供迭代初值，则 矩形迭代公式为 矩形迭代公式为
LL
上页
下页
例1 用欧拉公式求解初值问题 2x ′ = y− (0 < x < 1), y y y(0) = 1.