第五章 常微分方程初值问题数值解法
常微分方程初值问题的数值解法
1 1 2 1 , 2 p 2
这里有 3 个未知 数， 2 个方程。
存在无穷多个解。所有满足上式的格式统称为2阶龙格 - 库 塔格式。注意到，p 1, 1 2 1 就是改进的欧拉法。
Step 1: 将 K2 在 ( xi , yi ) 点作 Taylor 展开
K 2 f ( xi ph, yi phK1 ) f ( xi , yi ) phf x ( xi , yi ) phK1 f y ( xi , yi ) O( h2 )
y( xi ) phy( xi ) O(h2 )
d f ( x, y) dx 首先希望能确定系数 1、2、p，使得到的算法格式有 2阶 dy 精度，即在 yi y( xi ) 的前提假设下，使得 f x ( x, y) f y ( x, y) dx Ri y( xi 1 ) yi 1 O(h3 ) f x ( x, y) f y ( x, y) f ( x, y) y( x )
y( x0 ) y0 yk 1 yk h f ( xk , yk 1 )
, k 0,1,...
隐式欧拉法的求解: 利用迭代的思路进行.
yi 1 yi hf ( xi , yi 1 )
变换为
y
( k 1) i 1
yi hf ( xi , y )
y i 1 K1 K2
1 1 y i h K 1 K 2 2 2 f ( xi , yi ) f ( xi h, yi hK 1 )
步长一定是一个h 吗？
§2 Runge-Kutta Method
常微分方程初值问题数值解法
常微分方程初值问题的数值解法在自然科学、工程技术、经济和医学等领域中，常常会遇到一阶常微分方程初值问题：(,),,(),y f x y a x b y a y '=≤≤⎧⎨=⎩ (1) 此处f 为,x y 的已知函数，0y 是给定的初始值。
本章讨论该问题的数值解法，要求f 在区域{(,)|,}G x y a x b y =≤≤<∞内连续，并对y 满足Lipschitz 条件，从而初值问题(1)有唯一的连续可微解()y y x =，且它是适定的。
1 几个简单的数值积分法1.1 Euler 方法（1）向前Euler 公式（显式Euler 公式）10(,),0,1,2,,(),n n n n y y hf x y n y y a +=+=⎧⎨=⎩(2) 其中h 为步长。
由此便可由初值0y 逐步算出一阶常微分方程初值问题(1)的解()y y x =在节点12,,x x 处的近似值12,,y y 。
该公式的局部截断误差为2()O h ，是一阶方法。
（2）向后Euler 公式（隐式Euler 公式）1110(,),0,1,2,,(),n n n n y y hf x y n y y a +++=+=⎧⎨=⎩(3) 这是一个隐格式，也是一阶方法。
这类隐格式的计算比显格式困难，一般采用迭代法求解。
首先用向前Euler 公式提供迭代初值，然后迭代计算：(0)1(1)()111(,),(,),0,1,2,n n n n k k n n n n y y hf x y y y hf x y k +++++⎧=+⎨=+=⎩ (4)1.2 梯形方法1110[(,)(,)],2(),(0,1,2,)n n n n n n h y y f x y f x y y y a n +++⎧=++⎪⎨⎪=⎩= (5) 这也是一个隐格式，是二阶方法。
一般也采用迭代法求解。
迭代公式如下：(0)1(1)()111(,),[(,)(,)],0,1,2,2n n n n k k n n n n n n y y hf x y h y y f x y f x y k +++++⎧=+⎪⎨=++=⎪⎩ (6)1.3 改进的Euler 方法11110(,),[(,)(,)],0,1,2,,2(),n n n n n n n n n n y y hf x y h y y f x y f x y n y y a ++++⎧=+⎪⎪=++=⎨⎪=⎪⎩(7) 为了便于上机编程计算，(7)可改写为110(,),(,),0,1,2,,1(),2(),p n n n cn n p n p c y y hf x y y y hf x y n y y y y y a ++=+⎧⎪=+⎪⎪=⎨=+⎪⎪=⎪⎩(8) 该格式是显式，也是二阶方法。
第5章 - 常微分方程初值问题
xn b x
的近似曲线 ， 折线 p1 p2 pn , 称为欧拉折线，所以欧拉方法又称为
折线法。
二、隐式Euler方法和梯形方法 (Euler方法的改进) 1. 隐式 Euler方法 将 y( xk ) 在 x xk 1 点进行Taylor展开 y(k ) 2 2 y( xk ) y( xk 1 ) hk f ( xk 1 , y( xk 1 )) hk ,k [ xk , xk 1 ]，忽略 hk 项, 用 yk , yk 1 , fk 1 f ( xk 1 , yk 1 ) 分别近似 y( xk ), y( xk 1 ), f ( xk 1 , y( xk 1 )), 可得隐式欧拉方法: (13) yk 1 yk hk f ( xk 1 , yk 1 ), k 0,1,, n 1 或 yk 1 yk hk fk 1 , k 0,1,, n 1 说明: (1)隐式 Euler方法也可用向后差商,即用 dy y ( xk ) y ( xk 1 ) ， 近似微分 dx x x hk y y( x ) 或者用右矩数值求积公式来建立. (2)隐式 Euler方法(13)是关于 yk 1 的方程，若求 yk 1 需要 解方程(13).
y(x)存在且适定.
二、 初值问题数值解的基本概念 因为初值问题的数值解法是通过微分方程离散化而给出解在某些 为了讨论问题方便,引入以下概念。 离散点上(节点上)的近似值， 在 [a , b]上引入节点{ x k }n 0 , a x0 x1 x2 xn b, hk xk xk 1 k ba h , 则有 xk a kh( k 0,, n). ( k 1,2,, n) 称为步长。 常用等步长: n (1)，(2)的准确解记为y(x), y( xk ) 的近似解记为 yk ,记 f ( xk , yk ) fk .
常微分方程初值问题数值解法
数值解法的必要性
实际应用需求
许多实际问题需要求解常微分方程初值问题，如物理、 化学、生物、工程等领域。
解析解的局限性
对于复杂问题，解析解难以求得或不存在，因此需要 采用数值方法近似求解。
数值解法的优势
未来发展的方向与挑战
高精度算法
研究和发展更高精度的算法，以提高数值解的准确性和稳定性。
并行计算
利用并行计算技术，提高计算效率，处理大规模问题。
自适应方法
研究自适应算法，根据问题特性自动调整计算精度和步长。
计算机技术的发展对数值解法的影响
1 2
硬件升级
计算机硬件的升级为数值解法提供了更强大的计 算能力。
它首先使用预估方法（如欧拉方法）得到一个 初步解，然后使用校正方法（如龙格-库塔方法） 对初步解进行修正，以提高精度。
预估校正方法的优点是精度较高，且计算量相 对较小，适用于各种复杂问题。
步长与误差控制
01
在离散化过程中，步长是一个重要的参数，它决定 了离散化的精度和计算量。
02
误差控制是数值逼近的一个重要环节，它通过设定 误差阈值来控制计算的精度和稳定性。
能够给出近似解的近似值，方便快捷，适用范围广。
数值解法的历史与发展
早期发展
早在17世纪，科学家就开始尝 试用数值方法求解常微分方程。
重要进展
随着计算机技术的发展，数值 解法在20世纪取得了重要进展， 如欧拉法、龙格-库塔法等。
当前研究热点
目前，常微分方程初值问题的 数值解法仍有许多研究热点和 挑战，如高精度算法、并行计
软件优化
软件技术的发展为数值解法提供了更多的优化手 段和工具。
常微分方程数值解法欧拉法
)
f ( xn1, yn1)
hL
y(k ) n 1
yn1
L
hL
k 1
y(0) n 1
yn1
Q
hL 1,
y (k 1) n 1
yn1 (k
)
在迭代公式中取极限，有
yn1 yn h f ( xn1, yn1 ) 因此yn(k1)的极限就是隐式方程的解
几何意义
y
设已知曲线上一点 Pn (xn , yn ),过该 点作弦线，斜率为(xn+1 , yn +1 ) 点的 方向场f(x,y)方向,若步长h充分小， 可用弦线和垂线x=xn+1的交点近似 曲线与垂线的交点。
式。隐式公式不能直接求解，一般需要用Euler显式公式
得到初值，然后用Euler隐式公式迭代求解。因此隐式公
式较显式公式计算复杂，但稳定性好
y0 n1
yn
h
y(k 1) n1
yn
h
f (xn , yn )
f
( xn1 ,
y(k) n1
)
收敛性
y (k 1) n 1
yn1
h
f
( xn1,
y(k ) n 1
如何求解
解析解法：（常微分方程理论）
只能求解极少一类常微分方程；实际中给定的问题不一 定是解析表达式，而是函数表，无法用解析解法。
数值解法: 求解所有的常微分方程
计算解函数 y(x) 在一系列节点 a = x0< x1<…< xn= b
处的近似值 yi y( xi ) (i 1, ... , n)
y(xn1) y(xn ) hy(xn ) y(xn ) yn
y(xn1) yn1 yn h f (xn , yn )
第五章 常微分方程初值问题数值解法
则有
yn 1 yn hf ( xn , yn )
( 5.2 ) Euler格式
例5.1 用Euler格式解初值问题
2x y y y y (0) 1
取步长h=0.1.
(0 x 1)
Euler格式的具体形式为
y n 1 y n hf ( x n , y n ) 2 xn yn 0.1( yn ) yn 0.2 xn 1.1 yn yn
计算公式的精度 常以Taylor展开为工具来分析计算公式的精度. 为简化分析，假定yn是准确的，即在 yn y ( xn ) 的前提下估计误差 y ( xn 1 ) yn 1 Euler格式的局部截断误差 由 从而 局部截断误差
f ( xn , yn ) f ( xn , y ( xn )) y '( xn ) y ( xn 1 ) yn 1 y ( xn 1 ) ( yn hf ( xn , yn )) y ( xn 1 ) y ( xn ) hy '( xn )
y ( xn ), y ( xn 1 ), 的近似值 y1 , y2 , , yn , yn 1 ,
相邻两个节点的间距 h xi 1 xi 称为步长，步 长可以相等，也可以不等.本章总是假定h为定数， 称为定步长，这时节点可表示为
xn x0 nh , n 0,1, 2,
由f ( xn 1 , yn 1 ) f ( xn 1 , y ( xn 1 )) f y ( xn 1 , )( yn 1 y ( xn 1 )) f ( xn 1 , y ( xn 1 )) y '( xn 1 )（在xn点Taylor展开） h2 y '( xn ) hy ''( xn ) y '''( xn ) ... 2 3 2 h h 因此yn 1 y ( xn ) hy '( xn ) y ''( xn ) y '''( xn ) 2 4 h f y ( xn 1 , )( yn 1 y ( xn 1 )) 2 h2 h3 y ( xn 1 ) y ( xn ) hy '( xn ) y ''( xn ) y '''( xn ) 2 3!
第5章_常微分方程数值解法
（5.2.6）
由于方程关于 uk +1 是隐式形式，所以式（5.2.6）称为隐式 Euler 公式。前面显式和隐式 Euler 公式在计
u '(tk ) ≈
得到的递推公式：
u (tk +1 ) − u (tk −1 ) 2h
（5.2.7）
uk +1 ≈ uk −1 + 2hf (tk , uk )
在计算 uk +1 时，需要用到前两步结果 uk −1 , uk ，称为两步法公式。 （2）积分近似方法 将（5.2.1）式的微分方程写成 du = f (t , u )dt ，在区间 [tk , tk +1 ] 上积分，有：
5.2.2 Runge-Kutta 方法 Euler 方法比较简单，但它的收敛阶数低。可以利用 Taylor 展开式构造高阶的单步方法。Euler 公式 可以看成是由一阶 Taylor 展开式得到的，所以应用高阶 Taylor 展开就可以得到高阶单步法。例如：将 u (tk +1 ) 在 tk 处作 q 阶 Taylor 展开：
dy = a − by (t ) dt
是一阶常微分方程，而
2 ∂ 2 u ( x, t ) 2 ∂ u ( x, t ) a = ∂t 2 ∂x 2
（5.1.1）
（5.1.2）
是二阶偏微分方程。 所有使微分方程成为等式的函数，都是微分方程的解；在 n 阶微分方程中，将微分方程的含有 n 个任 意常数的解称为该微分方程的通解。为确定微分方程通解中的任意常数而需要的条件称为定解条件；定解 条件可以分为初始条件和边界条件两类。由微分方程和定解条件一起构成的问题称为微分方程定解问题。 根据定解条件的不同，常微分方程分为初值问题和边值问题；若定解条件是描述函数在一点（或初始 点）处状态的，则称为初值问题，一阶常微分方程初值问题的一般形式为：
常微分方程初值问题解法
为了克服欧拉方法精度不足的问题，可以对方法进行改进。一种常见的方法是使用更高阶的离散近似，例如使用 二阶或更高阶的离散化公式。这些改进可以减小数值误差，提高解的精度。
龙格-库塔方法
总结词
龙格-库塔方法是求解常微分方程初值问题 的一种高精度和高稳定性的数值方法。
详细描述
龙格-库塔方法是一种迭代方法，通过构造 一系列近似解来逼近微分方程的精确解。该 方法采用多步策略，每一步使用微分方程的 离散近似来更新未知数的值，同时考虑了更 多的信息，从而提高了数值解的精度和稳定 性。龙格-库塔方法在许多领域都有广泛的 应用，如物理、工程和科学计算等。
初值问题的定义
定义
常微分方程的初值问题由一个微分方程 和一个初始条件组成。给定一个初始状 态，我们需要找出该状态随时间变化的 规律。
VS
形式
dy/dt = f(t, y) with y(t0) = y0，其中f是 关于时间t和状态y的函数，t0是初始时间， y0是初始状态。
02
初值问题的解法
欧拉方法
05
结论与展望
研究成果总结
数值解法
常微分方程初值问题数值解法是当前研究的热点，包括欧拉法 、龙格-库塔法等多种方法，这些方法在精度和稳定性方面取
得了显著进展。
稳定性分析
对于数值解法的稳定性分析，研究者们通过分析数值解法 的收敛性和误差估计，为算法的改进提供了理论支持。
实际应用
常微分方程初值问题在物理、工程、生物等领域有广泛的应用 ，研究成果在实际问题中得到了验证，为解决实际问题提供了
04
实际应用与案例分析
物理问题中的应用
1 2 3
自由落体运动
描述物体在重力作用下的运动轨迹，可以通过常 微分方程求解物体在不同时刻的速度和位置。
常微分方程初值问题的数值解法60608
5.2Adams方法
显式Adams方法 隐式Adams方法
5.3预测-校正方法
y(0)n+1 = yn + h f(tn,yn)
(5.20)
y n+1 = (1/2)[y(0)n+1 + yn + hf(tn+1, y(0)n+1 )] (5.21)
n = 0,1,2,…N-1 ; y0= η
定理4 若Ф(t,y,h) 对于 a ≤t ≤ b, 0 ≤ h≤ h0以及一 切实数y，关于t,y,h满足Lipschitz条件，则单步 法(4.1)是稳定的。
定义4 对给定的微分方程和给定的步长h，若有
单步法（显式或隐式）计算yn时有大小为δ的误
差，即计算得yn` = yn + δ,而引起其后值ym(m>n)
(5.20)起预测y n+1 的作用， (5.21)起校正作用。 记f(i)n=f(tn,y(i)n).用P表预测过程，C表校正过程，E表计 算f的过程
P: y(0)n+1 = yn + hfn, E: f(0)n+1 = f(t n+1 , y(0)n+1 ) C: y n+1 = yn + (h/2)(fn+f(0)n+1)
5.4Hamming方法
Milne方法 建立线性多步法的待定系数法 Hamming方法
7线性多步法的相容性、收敛性和稳定 性
定义1 若求解初值问题(1.1)的线性k步法(5.1)至少是一阶 方法，则称他们是相容的。 记 ρ(λ)= αk λk+ αk -1 λk-1+…+ α1 λ + α0， σ(λ)= βk λk+ βk-1 λk-1+…+ β1 λ + β0。 他们由线性k步法(5.1)完全确定。反之，若给定了ρ(λ)和 σ(λ)，则他们唯一确定一个线性k步法。我们称ρ(λ)为线 性k步法(5.1)的特征多项式。
常微分方程初值问题的数值解法
常微分方程初值问题数值解法初值问题：即满足初值条件的常微分方程的解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乘以斜率；而这个斜率用区间上的端点和中点的斜率的算数平均来逼近。
常微分方程初值问题数值解法综述
常微分方程初值问题数值解法朱欲辉(浙江海洋学院数理信息学院, 浙江舟山316004)[摘要]：在常微分方程的课程中讨论的都是对一些典型方程求解析解的方法. 然而在生产实际和科学研究中所遇到的问题往往很复杂, 在很多情况下都不可能给出解的解析表达式. 本篇文章详细介绍了常微分方程初值问题的一些数值方法, 导出了若干种数值方法, 如Euler法、改进的Euler法、Runge—Kutta法以及线性多步法中的Adams显隐式公式和预测校正公式, 并且对其稳定性及收敛性作了理论分析. 最后给出了数值例子, 分别用不同的方法计算出近似解, 从得出的结果对比各种方法的优缺点.[关键词]：常微分方程; 初值问题; 数值方法; 收敛性; 稳定性; 误差估计Numerical Method for Initial-Value ProblemsZhu Yuhui(School of Mathematics, Physics, and Information Science,Zhejiang Ocean University, Zhoushan, Zhejiang 316004)[Abstract]: In the course about ordinary differential equations, the methods for analytic solutions of some typical equations are often discussed. However, in scientific research, the problems are very complex and the analytic solutions about these problems can 'bte expressed explicitly. In this paper, some numerical methods for the initial-value problems are introduced. these methods include Euler method, improved Euler method, Runge-Kutta method and some linear multistep method (e.g. Adams formula and predicted-corrected formula). The stability and convergence about the methods are presented. Some numerical examples are give to demonstrate the effectiveness and accuracy of theoretical analysis. [Keywords]: Ordinary differential equation; Initial-value problem; Numerical method;Convergence; Stability ; Error estimate1刖言自然界和工程技术中的很多现象，例如自动控制系统的运行、电力系统的运行、飞行器的运动、化学反应的过程、生态平衡的某些问题等，都可以抽象成为一个常微分方程初值问题•其真解通常难以通过解析的方法来获得，至今有许多类型的微分方程还不能给出解的解析表达式，一般只能用数值的方法进行计算•有关这一问题的研究早在十八世纪就已经开始了，特别是计算机的普遍应用，许多微分方程问题都获得了数值解，从而能使人们认识解的种种性质及其数值特征，为工程技术等实际问题提供了定量的依据•关于常微分方程初值问题的数值计算方法，许多学者己经做了大量的工作• 1768年, Euler提出了关于常微分方程初值问题的方法，1840年，Cauchy第一次对初值问题进行了仔细的分析，早期的常微分方程数值解的问题来源于天体力学•在1846年，当Adams还是一个学生的时候，和Le Verrier —起根据天王星轨道中出现的己知位置，预测了它下一次出现的位置.1883年，Adams 提出了Adams —Bashforth 和Adams —Moulton 方法.Rull （1895 年）、Heun （1900年）和Kutta （1901 年）提出Runge.Kutta 方法.二十世纪五十年代，Dahlquist建立了常微分方程数值解法的稳定性理论，线性多步法是常微分方程初值问题的一种数值方法.由于通常的数值方法，其绝对稳定区域是有限的，不适用于求解刚性常微分的初值问题.刚性微分方程常常出现于航空、航天、热核反应、自动控制、电子网络及化学动力学等一系列与国防和现代化建设密切相关的高科技领域，具有无容置疑的重要性.因此，刚性微分方程的研究工作早在二十世纪五十年代就开始了，1965年，在爱丁堡举行的IFIP会议后，更进一步地认识刚性方程的普遍性和重要性.自从六十年代初，许多数值分析家致力于探讨刚性问题的数值方法及其理论，注意到刚性问题对传统数值积分方法所带来的挑战.这一时期，人们的研究主要集中在算法的线性稳定性上,就是基于试验方程y，- C）数值解的稳定性研究.在此领域发表了大量的论文，取得了许多重要的理论成果.例如，1963年，Dahlquist给出A稳定性理论，1967年, Widlund给出A（〉）一稳定性理论，1969年，Gear将A —稳定性减弱，给出刚性（Stiff）稳定性理论，并找到了当k空6的k步k阶的刚性稳定方法，1969年Dill找到刚性稳定的7阶和8阶以及1970年Jain找到刚性稳定的9阶到11阶，但可用性没有检验.这些稳定性理论和概念都是在线性试验方程的框架下推导出的，从严格的数学意义上来说，这些理论只适用于常系数线性自治系统.但从实用的观点来说，这些理论无疑是合理和必要的，对刚性问题的算法设计具有重要的指导意义•在八十至九十年代，国内也有一些学者研究线性理论，主要有匡蛟勋、陈果良、项家祥、李寿佛、黄乘明、李庆扬和费景高等线性理论虽然对一般问题具有指导作用，但其不能作为非线性刚性问题算法的稳定性理论研究基础•为了将线性理论推广到非线性问题中，人们开始对非线性模型问题进行研究•但是,早期文献主要致力于数值方法基于经典Lipschitz条件下的经典收敛理论,即认为良好的稳定性加上经典相容性和经典相容阶就足以描述方法的整体误差性态•直到1974年, Prothero和Robi nson首先注意到算法的经典误差估计由于受刚性问题巨大参数的影响而严重失真，产生阶降低现象，这时人们认识到经典收敛理论对于非线性刚性问题以及线性模型的不足•于是，1975年，Dahlquist和Butcher分别提出了单支方法和线性多步法的G一急定概念和B稳定概念.这两个概念填补了非线性稳定性分析理论，引起了计算数学家们的极大关注，在上述理论的基础上，1975年至1979年，Burrage和Butcher提出了AN 一稳定性与BN 一稳定性概念，并相应地建立了基本的 B 一稳定及代数稳定理论.1981至1985年，Frank, schneid和ueberhuber建立Runge一kutta方法的旷收敛理论.B稳定与B一收敛理论统称B—理论，它是常微分方程数值解法研究领域的巨大成就之一，是刚性问题算法理论的突破性进展，标志着刚性问题研究从线性向非线性情形深入发展.国内也有众多学者致力于 B 一理论的研究，如李寿佛、曹学年等.1989年，李寿佛将Dahlquist的G—稳定概念推广到更一般的（C,P,Q）代数稳定，克服了G急定的线性多步法不能超过二阶的限制.对于一般线性方法，李寿佛建立了一般线性方法的（K,P,Q）稳定性理论及（K,P,Q）弱代数稳定准则和多步Runge- Kutta法的一系列代数准则.此外，Dahquist, Butcher和Hairer分别深刻地揭示了单支方法、一般线性方法和Runge.Kutta方法线性与非线性稳定性之间的内在关系.为了求解刚性微分方程，不少文献中构造含有稳定参数的线性多步方法，利用适当选择稳定参数来扩大方法的稳定区域.所有改进的思想，都是通过构造一些特殊的显式或隐式线性多步法，使其具有增大的稳定域，或使A（＞）—稳定的〉角增大.八十年代，就成为国内外学者所研究的一个课题，学者主要有Rodabaugh和Thompson、Feinberg、李旺尧、李寿佛、包雪松、徐洪义、刘发旺、匡蛟勋、项家祥、蒋立新、李庆扬、谢敬东和李林忠等.当前国内外研究刚性问题的一个主要趋势就是在E一理论指导下寻找更为有效的新算法.另一个发展趋势就是力图突破单边lipschitz常数和内积范数的局限，建立比B—理论更为普遍的定量分析收敛理论.近年来，刚性延迟系统的算法研究成为刚性问题的另一个热点研究领域，张诚坚将Burrage等人创立的针对刚性常微分系统的B—理论拓展到非线性刚性延迟系统常微分方程的数值算法发展到今天己有了线性多步法、龙格一库塔法和在此基础上发展起来的单支方法、分块方法、循环方法、外推法、混合方法、二阶导数法以及各种常用的估校正算法.其中经常用到的线性多步法公式有Euler公式、Heun公式、中点公式、Milne公式、Adams公式、simpson公式、Hamming公式,Gear方法、Adams预估一校正法和Mile预估一Hamming校正法公式等，此外还包含许多迄今尚末探明的新公式.Burage曾将线性多步法和Runge—Kutta法比作大海中的两座小岛，在浩瀚的汪洋之中，还有许多到现在没有发现的新方法.本篇文章详细介绍了常微分方程初值问题的一些数值方法，导出了若干种数值方法，如Euler法、改进的欧拉法、显式龙格-库塔法、隐式龙格-库塔法以及线性多步法中的Adams显隐式公式和预测校正公式，并且对其稳定性及收敛性作了理论分析.最后给出了数值例子，进行了计算机程序算法的分析与实现，以计算机的速度优势来弥补计算量大的不足. 从得出的结果对比各种方法的优缺点.2常微分方程初值问题的数值解法2.1数值方法的基本思想考虑一阶常微分方程初值问题dxf(x,y),x [a,b],dyy(x c)=y o (2.1) 的数值解法，其中f是x和y的己知函数，y0是给定的初始值.对于常微分方程初值问题(2.1)数值解法，就是要算出精确解y(x)在区间[a,b]上的一系列离散节点x0x^ I : x n4x n= b的函数值y(x)) , y(xj, Hl,y(x n)的近似值y0, %, •…，y n.相邻两个节点的间距x i称为步长，这时节点也可以表示为X ih (i =1,2, III,n).数值解法需要把连续性的问题加以离散化,从而求出离散节点的数值解.通常微分方程初值问题(2.1)的数值方法可以分为两类代入式(2.3),并用y n 近似代替式中y(x n )即可得到梯形公式hy ny n -[ f (X n ,『n ) f (X n 1, Y n 1 )]由于数值积分的梯形公式比矩形公式精度高 ，所以梯形公式(2.4)比欧拉公式(2.2)的精度高一个数值方法.式(2.4)的右端含有未知的 y n1,它是关于y n d 的函数方程，这类方法称为隐式方法(1)单步法计算y(x)在x =X n 1处的值仅取决于x n 处的应变量及其导数值(2)多步法-----计算y(x)在X 二X n 1处的值需要应变量及其导数在 X n 1之前的多个网个节点出的值.2.2 Euler 方法 2.2.1 Euler 公式若将函数y(x)在点x n 处的导数y '(xj 用两点式代替，即y(X n ),竺4血，h再用y n 近似地代替y(x n ),贝U 初值问题(2.1)变为y n 1 =y n hf (X n ,y n )I y^= y(x o),0,1,2，川(2.2)式就是著名的欧拉(Euler)公式.(2.2)2.2.2 梯形公式欧拉法形式简单精度低,为了提高精度，对方程y = f (x, y)的两端在区间[X n ,X n 1]上积分得xny(X n+)=y(X n ) +Jf[x,y(x)]d(2.3)改用梯形方法计算其积分项xn 1x f[x, y(x)]dx xPmXn^X n )) "-,"“1))](2.4)223改进的欧拉公式梯形公式实际计算时要进行多次迭代，因而计算量较大•在实用上，对于梯形公式(2.4)只迭代一次,即先用欧拉公式算出 y n1的预估值了宀，再用梯形公式(2.4)进行一次 迭代得到校正值y n 勺，即+y n 1 =y n hf(X n ,y n )<h(2・5)| yn 1 = yn* ^[ f(Xn , yn)f (xn 1, yn -1)]2.2.4欧拉法的局部截断误差衡量求解公式好坏的一个主要标准是求解公式的精度 ，因此引入局部截断误差和阶数概念.定义2.1 在y n 准确的前提下，即y n 二y(x n )时，用数值方法计算y n1的误差 尺1二y(X n 1)- y n1，称为该数值方法计算 yn 1时的局部截断误差•对于欧拉公式，假定y n 二y(x n )，则有y n 1 =y(X n ) h[f(X n ，y(X n ))Hy(X n ) hy '(X n )而将真解y(X)在X n 处按二阶泰勒展开式有, h 2................Y n ^y(X n ) hy '(X n ) = 2( )，X ，X . 1)2!因此有定义2.2若数值方法的局部截断误差为0(h p1)，则称这种数值方法的阶数是长(h<1)越小，P 越高，则局部截断误差越小，计算精度越高2.3 Runge — Kutta 方法2.3.1 R unge — Kutta 方法的基本思想欧拉公式可改写成y n 1 = y n K^ f (X n , y n )『(X n 1)Tn 1£2!Y ''() P.步则y n 1的表达式与y(x n1)的泰勒展开式的前两项完全相同,即局部截断误差为0(『).改进的欧拉公式又可改写成hy^ = y n十2(心+心)& = f (X n 朴y n +hkj1 L上述两组公式在形式上有一个共同点：都是用f(x,y)在某些点上值的线形组合得出y(X.,)的近似值y n d,而且增加计算f(x,y)的次数，可提高截断误差的阶•如欧拉公式每步计算一次f(x, y)的值，为一阶方法•改进的欧拉公式需计算两次 f (x, y)的值，它是二阶方法•它的局部截断误差为0(h3).于是可考虑用函数f (x, y)在若干点上的函数值的线形组合来构造近似公式，构造时要求近似公式在(x n, y n)处的泰勒展开式与解y( x)在x n处的泰勒展开式的前面几项重合，从而使近似公式达到所需要的阶数•既避免求偏导，又提高了计算方法精度的阶数•或者说，在［x^x^］这一步内多预报几个点的斜率值,然后将其加权平均作为平均斜率•则可构造出更高精度的计算格式，这就是Runge—Kutta方法的基本思想•2•3•2 Runge — Kutta 方法的构造一般地，Runge —Kutta方法设近似公式为(下面的公式修改了)' p%卅= yn+h瓦qK i! 7“K^ f (x n，y n) (1 =2,3，|||，p)i 4K j = f(焉 +ahy n + 运b j K j)I. j^其中a j，q，G都是参数，确定它们的原则是使近似公式在y m处的泰勒展开式与y(x)在x n处的泰勒展开式的前面的项尽可能多的重合，这样就使近似公式有尽可能高的精度•以此我们可以通过一个复杂的计算过程得出常用的的三阶和四阶Runge—Kutta 公式yn^ = y n +：(心+4K2 +K3)6K i=f(X n,y n)1 h hK2 = f(X n +q，y n +君1)、K3 = f(Xn+h,y n—hKi+2hK2)和•h*卑7.+6(心+2心+2心+心)K i = f (X n, y n)h hK2 = f (Xn + yn + K i) (2.6)2 2h hK3 = f (X n y n K2)2 2K^ f (X n h, y n hK3)式(2.6)称为经典Runge —Kutta方法.2.4线性多步法在逐步推进的求解过程中，计算y n1之前事实上已经求出了一系列的近似值y0, %, y n.如果充分利用前面多步的信息来预测ym,则可期望获得较高的精度，这就是构造多步法的基本思想.线性k步方法的一般公式为k 4 k4y n1=»a j y n4 h' b j f(X n4,y n_j) (2.7)j=e j=J其中a j, b j均为与n无关的常数，|a k』+b k j HO.当b_^ =0时为显格式；当时1为隐格式.特别当k = 1, a0 = b0 = 1,b^ = 0时为Euler 公式；当k = 1,a0 =1,b0 =b_i = 2 时为梯形公式.定义2.3称kJ kJR n 1二y(X n 1)—[y n 1八2 b j f (人」小」)]j」j =J为k步公式(2.7)在x n d处的局部截断误差•当R n l =0(h P )时称式(2.7)是p阶的.应用方程y (x)二f (x, y(x))可知局部截断误差也可写成为k A kJI R n+ = y(X n 卅)—[y n 卑二瓦a j y n」+h 瓦b j y (X n」)] j=0 j=J定义2.4如果线性k步方法(2.7)至少是1阶的，则称是相容的；如果线性k步法(2.7)是 p阶的，贝U称是p阶相容的.2.4.1A dams 外插法将微分方程生=f (x, y),x [a,b],dy的两端从x n到x n 1进行积分，得到x n +y(x nJ=y(x) x f(x,y(x))dx (2.8)xn我们用插值多项式代替右端的被积函数.Adams外插法选取k个点x n, x n斗=111, x n乂十作为插值基点构造f (x, y)的k-1阶多项式Adams外插法的计算公式为kAY n ^Y nh a^ j f m7j =0,1,HI根据此递推公式，可逐个的计算q (j =0,1,|||),表2.1给出了 aj 的部分数值2.4.2 Adams 内插法根据插值理论知道，插值节点的选择直接影响着插值公式的精度，同样次数的内插公式的精度要比外插公式的高 .仍假定已按某种方法求得问题 (2.1)的解y(x)在xXi =0,1,l||n)处的数值y i ,并选取插 值节点x n 」4,III, X n iH ,X n+p P 是正整数，用Lagrange 型插值多项式L P iL(x)构造 y= f( x, y)可以导出解初值问题(2.1)的Adams 内插公式：冷1丫风十)=y(X n ) + J L nku (x)dx(2.9)X n当p = 0时上式就退化为内插公式.内插公式(2.9)除了包含f 在X n»1,lH,X n 处的已知值外，还包含了在点X n ,川，X n.p ,处的未知值.因此内插公式(2.9)只给出了未知量y n ,IILy n.p 的关系式，实际计算时仍需要其中a j 满足如下代数递推式aj1 1才』3a j「'解联立方程组• p =1的内插公式是最适用的，L^k j(x)采用Newt on向后插值公式得到Adams内插公式ky n y n a j j f mj=0其中系数a］定义为(-1)j■dij =0,1,HI其中a j满足如下代数递推式：_ 1,H0_ 0,j 0根据此递推公式，可逐个的计算a］(j =0,1,1"),表2.2给出了a］的部分数值2.22.5算法的收敛性和稳定性2.5.1相容性初值问题(2.1):f(x,y),x [a,b], dyy(«) = y。
常微分方程初值问题的数值解法
就得到初值问题(7.1),(7.2)的解y(t )的解析表达式。然而
在实际问题和科学研究中所遇到的微分方程往往很复
杂，很多情况下不可能求出它的解析解。有时侯即使
能求出解析解，也会由于很难从解析解中计算函数y(t )
的值而不实用。
例如，容易求出初值问题
y' 1 y cos t,0 t T
注意：这是“折
线法”而非“切
yN
线法”除第一个
点是曲线切线外，
其他点不是切线
而是折线(如右 图所示)。
y2 yy10
t0 a t1 t2
tN b x
§7.2.1 显式单步法的一般形式 显式单步法的一般形式是
yn1 yn h (tn , yn , h), n 0,1, , M 1 (7.2.4)
普希兹)条件。常数L称为函数f 在D0中的Lipschitz常数。
例1 函数f (t, y) t y 在区域D0 (t, y) | 1 t 2, 3 y 4
关于y满足Lipschitz条件，相应的Lipschitz常数可取为L 2
3 存在性定理 定理1 设函数f (t, y)在凸集D R2中有定义，若存在常数
(7.4)
若k 2,则数值解法(7.3)统称为多步法，或具体称 为k步法。
显示法、隐式法与截断误差 若在差分方程(7.3)中，ynk能表示为tn , yn, yn1, , ynk-1, h 的显函数，即
ynk G(tn , yn , yn1, , ynk-1, h), n 0,1, , M - k (7.5)
y0 yn1
Hale Waihona Puke y(t0 yn)
常微分方程初值问题的数值解法
1
1、基本概念： 基本概念：
•一阶常微分方程初值问题是： 一阶常微分方程初值问题是： 一阶常微分方程初值问题是 y ′ =f(x,y) (1.1) y(x0)=y0 (1.2) 其中f是已知的xy平面上某个区域D上连续函数，式（1.1）是微分 方程，有无穷多解，式（1.2）是确定解的初始条件。 如一元函数y(x)对一切a≤x≤b 满足 （1） （x,y(x))∈D （2） y(x0)=y0 （3） y′存在，且y′(x)=f(x,y(x)) 则称y(x)是初值问题（1.1）、（1.2）在[a,b]上的解。 定义1.1 如果存在正常数L>0,使得对一切x∈[a,b]及y1 、y2 ∈[c,d] 有 |f(x,y1)-f(x,y2)|≤ L|y1- y2| (1.3) 则称f(x,y)满足对y的lipschitz(李普希滋）条件，其中L称为 lipschitz常数。
~ = y + hf ( x , y ) y k +1 k k k y k +1
2011-12-12
h y = y k + ( f ( x k , y k ) + f ( x k +1 , ~k +1 )) 2
12
Euler公式的几何意义 公式的几何意义
Y=y(x)
a
b
2011-12-12
13
例题：用Euler的数值解(取步长h=0.1计算到y3)：
2011-12-12 4
求初值问题，是给出它的解在某些节点数值的近似值，这称为 数值离散方法，若要求出在[a,b} 上的离散解，引入点列{xK} ， 其中 xk= xk-1 + hk-1 ，k=1，2，3….. xk称为节点， hk 称为步长。 通常，步长h不变，取为等距步长 hk =h= （b-a）/N ，N为等份 区间[a,b]分割数,节点变为等距节点，此时有 kh, xk= xk-1 + h= x0 + kh 记式（1.1），（1.2）的准确解y(x)在节点xk 之值为y(xk) ，而 记y(xk)的近似值为yk ，又记fk=f（xk ， yk ），
常微分方程初值问题的数值解法
常微分方程初值问题的数值解法在实际应用中，对于某些微分方程，我们并不能直接给出其解析解，需要通过数值方法来求得其近似解，以便更好地理解和掌握现象的本质。
常微分方程初值问题（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)的准确解。
《数值分析》第五章课件
取 h = 0.2 ，要求保留六位小数.
校正： cn+1 = y n + 2 ( y n' + mn' +1 )
解：Euler 迭代格式为
校正的改进：
1 y n +1 = c n +1 + ( p n +1 − c n+1 ) 5
yk +1 = yk + 0.2(− yk − xk yk2 ) = 0.8 yk − 0.2 xk yk2
差分方程：关于未知序列的方程.
例如： y n +3 = 5 y n + 2 − 3 y n +1 + 4 y n
例如： y ' ' ( x) − a ( x) y '+b( x) y + c( x) = 0
3
4
微分方程的应用情况
实际中，很多问题的数学模型都是微分方程. 常微分方程作为微分方程的基本类型之一，在 理论研究与工程实际上应用很广泛. 很多问题 的数学模型都可以归结为常微分方程. 很多偏 微分方程问题，也可以化为常微分方程问题来 近似求解.
且
可得，
y(xn+1) − yn+1 = hf y (xn+1,η)[ y(xn+1) − yn+1] − h2 '' y (xn ) + O(h3 ) 2
f (xn+1, y(xn+1)) = y' (xn+1) = y' (x n ) + hy'' (xn ) + O(h2 )
19
20
2 考虑到 1 − hf y ( xn+1 ,η ) = 1 + hf y ( xn+1 ,η ) + O(h ) ，则有
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
y=y(x) Pn
x0
x1
xi
xi+1
xn
这样,从x0逐个算出 x1 , x2 , xn 对应的数值解
y1 , y2 , yn
从图形上看,就获得了一条近似于曲线y=y(x) 的折线 PP 1 2P 3
P n.
Euler格式算法实现 （1）计算步骤 ba ① 输入 a, b, h, f ( x, y), x0 , y0 ,并计算出 N ; h ② 使用Euler格式进行计算
向前推进，描述这类算法，要求给出用已知信息 差分格式 yn , yn1 , yn2 , 计算 yn1 的递推公式. 建立这类递推公式的基本方法是在这些节点上用数 值积分、数值微分、泰勒展开等离散化方法，对初 值问题 y f ( x, y ) y ( x0 ) y 0 中的导数 y 进行不同的离散化处理.
输出 x1, y1 n
n=N ? n+1 n
y 结束
x1 x0 y1 y0
例5.1 用Euler格式解初值问题
2x y y y y(0) 1
取步长h=0.1.
(0 x 1)
Euler格式的具体形式为
yn1 yn hf ( xn , yn ) 2 xn yn 0.1( yn ) yn 0.2 xn 1.1yn yn
y x y
2
2
这个一阶微分方程就不能用初等函数及其积分来 表达它的解.
再如，方程
y y y (0) 1
的解 y e x ,虽然有表可查,但对于表上没 有给出 e x 的值,仍需插值方法来计算
从实际问题当中归纳出来的微分方程，通常主要依 靠数值解法来解决.本章主要讨论一阶常微分方程初 值问题
5.2.2 后退的Euler格式
对 y '( x) f ( x, y( x)),在点xn+1处为
y '( xn1 ) f ( xn1, y( xn1 ))
y ( xn 1 ) y ( xn ) y ( xn 1 ) y ( xn ) 取 y '( xn 1 ) xn 1 xn h
则有 y( xn1 ) y( xn ) hf ( xn1, y( xn1 )) 设y(xn)的近似值yn已知，用它带入上式右端进行计
算，并取计算结果yn+1作为y(xn+1)的近似值,则可得
yn1 yn hf ( xn1 , yn1 )
( 5.3 )
后退的Euler格式
隐式方程(5.3) yn1 yn hf ( xn1 , yn1 )常用迭代法求解 (0) (0) 设用Euler格式 yn 给出迭代初值 y hf ( x , y ) y 1 n n n n 1
长可以相等，也可以不等.本章总是假定h为定数，
称为定步长，这时节点可表示为
,
xn x0 nh, n 0,1, 2,
数值解法需要把连续性的问题加以离散化，从而求 出离散节点的数值解.
对常微分方程数值解法的基本出发点就是离散 化.其数值解法有个基本特点，它们都采用“步进 式”，即求解过程顺着节点排列的次序一步一步地
对R内任意两个 y1 , y2 都成立,则方程( 5.1 )的解 y y ( x) 在a, b上存在且唯一.
数值方法的基本思想 对常微分方程初值问题(5.1)式的数值解法，就 是要算出精确解y(x)在区间a,b上的一系列离散节
点 x1 x2
xn xn1 处的函数值 y( x1 ), y( x2 ), y( xn ), y( xn1 ), 的近似值 y1, y2 , , yn , yn1, 相邻两个节点的间距 h xi 1 xi 称为步长，步
从而有
h2 h3 y( xn1 ) yn 1 y ''( xn ) y '''( xn ) 2 3 hf y ( xn 1 , )( y( xn 1 ) yn 1 ) 移项整理得 1 h2 h3 y( xn1 ) yn1 ( y ''( xn ) y '''( xn ) 1 hf y ( xn1 , ) 2 3 1 由 1 x 2 x3 可知 1 x 1 2 1 hf y ( xn1 , ) (hf y ( xn 1 , )) 1 hf y ( xn1 , )
以 f ( xn , yn )为斜率作直线 y yn f ( xn , yn )(x xn )
当 x xn1 时,得 yn1 yn f ( xn , yn )(x x1 xn ) 取 y( xn1 ) yn1
P1 P0 P1
Pi+1 Pn Pi Pi+1 Pi
程的阶数.如果未知函数y及其各阶导数
y, y,, y
( n)
都是一次的,则称它是线性的,否则称为非线性的.
在《常微分方程》中，对于常微分方程的求解，
给出了一些典型方程求解析解的基本方法，如可分
离变量法、常系数齐次线性方程的解法、常系数非
齐次线性方程的解法等.但能求解的常微分方程仍 然是有限的，大多数的常微分方程是不可能给出解 析解. 譬如
当 x x2 时,得 y2 y1 f ( x1 , y1 )(x2 x1 ) 由此获得了P2的坐标.
P1 P0 P1
Pi+1 Pn Pi Pi+1 Pi
y=y(x) Pn
x0
x1
xi
xi+1
xn
重复以上过程,就可获得一系列的点:P1,P1,…,Pn.对 已求得点 Pn ( xn , y n )
y f ( x, y ) y ( x0 ) y 0 的数值解法，首先要解决的问题就是如何对微分方 程进行离散化，建立求数值解的递推公式. 递推公式通常有两类，一类是计算yn+1时只用到xn+1, xn和yn，即前一步的值，因此有了初值以后就可以逐 步往下计算，此类方法称为单步法. 另一类是计算yn+1时，除用到xn+1,xn和yn以外，还要 用到 xn1, yn1, xn2 , yn2 , xnk , ynk , 即前面k步的值，此 类方法称为多步法.
即有 y '( x) f ( x, y( x))
在点xn处为 y '( xn ) f ( xn , y( xn ))
y( xn1 ) y( xn ) y( xn1 ) y( xn ) 取 y '( xn ) xn1 xn h
y( xn1 ) y( xn ) hf ( xn , y( xn )) yn1 yn 设y(xn)的近似值yn已知，用它带入上式右端进行计 算，并取计算结果yn+1作为y(xn+1)的近似值,则可得
计算公式的精度 常以Taylor展开为工具来分析计算公式的精度. 为简化分析，假定yn是准确的，即在 yn y( xn ) 的前提下估计误差 y( xn1 ) yn1 Euler格式的局部截断误差 由 从而
局部截断误差
f ( xn , yn ) f ( xn , y( xn )) y '( xn ) y( xn1 ) yn1 y( xn1 ) ( yn hf ( xn , yn )) y( xn1 ) y( xn ) hy '( xn ) h2 y ''( xn ) 2
(k ) 如果迭代过程收敛，则极限值 yn 1 lim yn 1 k
必满足隐式方程(5.3),从而获得后退的Euler格式的解.
后退的Euler格式的局部截断误差(假设 yn y( xn ) )
yn1 yn hf ( xn1 , yn1 ) y( xn ) h f ( xn1 , yn1 )
用它带入(5.3)的右端，使之转化为显式，直接计算得
(1) ( 0) yn y hf ( x , y 1 n n1 n1 )
然后再用 y
(1) n 1 ( 2) n1
带入(5.3)的右端，又有 (1) y yn hf ( xn1, yn 1 )
如此反复进行迭代，得
( k 1) (k ) yn y hf ( x , y 1 n n1 n1 ), k 0,1,
y1 y0 hf ( x0 , y0 ) x1 x0 h
③ 输出 x1 , y1，并使 x1 x0 , y1 y0 转到 ② 直至n > N 结束.
（ ） 格 式 的 流 程 图
开始 输入 a,b,h,x0, y0,f(x,y) N=(b-a)/h;n=1
2 Euler
y0+hf(x0,y0 ) y1 x0+h x1
y f ( x, y ) ( 5.1 ) y ( x ) y 0 0 在区间a≤x≤b上的数值解法. 可以证明,如果函数在带形区域 R={a≤x≤b,-∞＜y＜∞} 内连续，且关于y满足Lipschitz条件，即存在常数 L(它与x,y无关)使 f ( x, y1 ) f ( x, y2 ) L y1 y2
对于初值问题
§5.2 Euler方法 5.2.1 Euler格式
Euler方法是解初值问题的最简单的数值方法.
初值问题
y f ( x, y ) y ( x0 ) y 0 的解y=y(x)称为它的积分曲线.积分曲线上每一点 ( x, y ) 处的切线的斜率 y ( x) 等于函数 f ( x, y ) 在这 点的值.
由f ( xn1, yn1 ) f ( xn1, y( xn1 )) f y ( xn1,)( yn1 y( xn1 ))
f ( xn1, y( xn1 )) y '( xn1 )（在xn点Taylor展开） h2 y '( xn ) hy ''( xn ) y '''( xn ) ... 2 3 h 2 因此yn 1 y ( xn ) hy '( xn ) h y ''( xn ) y '''( xn ) 2 hf y ( xn 1 , )( yn 1 y ( xn 1 )) 2 3 h h y ( xn1 ) y ( xn ) hy '( xn ) y ''( xn ) y '''( xn ) 2 3!