研究生数值分析-第6章 微分方程数值解法

为方程的解。 一般称为方程的通解。 如果
y(0) 1 则有
y x 1
2
为方程满足定解条件的解。
解的图示
y' 2 x y (0) 1
方程的通解 微分关系（方程）

y x2 C 1 C

y x2 1
满足定解条件的解
本教材重点讨论定解问题(初值问题）
y( x1 ) y( x 2 ) y( x0 ) 2h
y( x2 ) y( x0 ) 2h f ( x1 , y( x1 ))
x0
x1
x2
yi 1 yi 1 2h f ( xi , yi ) i 1, ... , n 1
假设 yi 1 y( xi 1 ), yi y( xi ) , 则可以导出 Ri y( xi 1 ) yi 1 O(h3 ) 即中点公式具有 2 阶精度。
y' f ( x , y ) y ( x 0 ) y0
定解条件（初始条件）
是否能够找到定解问题的解取决于 f ( x , y )
仅有极少数的方程可以通过“常数变易法”、“可分
离变量法”等特殊方法求得初等函数形式的解，绝大
部分方程至今无法理论求解。
如
y' sin(xy ),
2
y' 1 x
第6章
§6.1 §6.2
常微分方程数值解法
引 言 欧拉方法
§6.3
龙格—库塔方法
§6.1 引Leabharlann 言微分方程数值解一般可分为：常微分方程数值解和偏微分 方程数值解。自然界与工程技术中的许多现象，其数学表达式 可归结为常微分方程（组）的定解问题。一些偏微分方程问题 也可以转化为常微分方程问题来（近似）求解。Newton最早采 用数学方法研究二体问题，其中需要求解的运动方程就是常微 分方程。许多著名的数学家，如 Bernoulli（家族），Euler、 Gauss、Lagrange和Laplace等，都遵循历史传统，研究重要 的力学问题的数学模型，在这些问题中，许多是常微分方程的 求解。作为科学史上的一段佳话，海王星的发现就是通过对常 微分方程的近似计算得到的。本章主要介绍常微分方程数值解 的若干方法。
算法有p 阶精度。
欧拉法的局部截断误差：
yi y( xi )
y '( xi ) f ( xi , y( xi ))
Ri y( xi 1 ) yi 1
[ y( xi ) hy( xi )
h2 2
y( xi ) O(h3 )] [ yi hf ( xi , yi )]
, yi r 1
§6.2 Euler方法 1、Euler格式
第一步：连续变量离散化
y' f ( x , y ) y ( x 0 ) y0
x x0 , x1 ,, xk , xn ,
p1
第二步：用直线步进
y1 y0 f ( x0 , y0 ) x1 x0
h2 3 y ( xi ) O(h ) 2
O(h )
2
Ri 的主项
/* leading term */
欧拉法具有 1 阶精 度。
例1: 用欧拉公式求解初值问题
2 y 2 xy y 0 1
（0 x 1.2）
取步长 h 0.1。 解: 应用Euler公式于题给初值问题的具体形式为：
方 法 显式欧拉 隐式欧拉 梯形公式 中点公式
简单 稳定性最好
精度提高 精度提高, 显式

精度低 精度低, 计算量大
计算量大 多一个初值, 可能影响精度

改进欧拉法 /* modified Euler’s method */
Step 1: 先用显式欧拉公式作预测，算出
yn1 yn hf xn , yn
y xi yi
0.009901 0.018462 0.024153 0.026320 0.025250 0.021852 0.017213 0.012262 0.007626 0.003642 0.000422 0.002053
2 这个初值问题的准确解为 y x 1 1 x ，
Step 2: 再将 y n 1
代入隐式梯形公式的右边作校正，得到
h yn 1 yn f xn , yn f xn 1 , yn 1 2


上式还常写成 1 yn 1 yn 2 (k1 k2 ) k1 hf ( xn , yn ) k2 hf ( xn h, yn k1 ) (n 0,1, 2,...) 该式称为改进Euler方法, 亦可写成 h yn 1 yn [ f ( xn , yn ) f ( xn h, yn hf ( xn , yn ))] 2
果74卷。
2、欧拉法的局部截断误差
定义2.1 在假设 yi = y(xi)，即第 i 步计算是精确的前提
下，考虑的截断误差 Ri = y(xi+1) yi+1 称为局部截断 误差 /* local truncation error */。
定义2.2 若某算法的局部截断误差为O(hp+1)，则称该
5、微分方程的数值解法需要解决的主要问题
(1) 如何将微分方程离散化，并建立求其数值解的迭代公式？
(2) 如何估计迭代公式的局部截断误差与整体误差？ (3) 如何保证迭代公式的稳定性与收敛性?
二、初值问题解的存在唯一性
考虑一阶常微分方程的初值问题 /* Initial-Value Problem */:
yi
1.000000 0.980000 0.941584 0.888389 0.825250 0.757147 0.688354 0.622018 0.560113 0.503642 0.452911 0.407783
y xi
0.990099 0.961538 0.917431 0.863069 0.800000 0.735294 0.671141 0.609756 0.552486 0.500000 0.452489 0.409836
L 0
s.t.
f ( x, y1 ) f ( x, y2 ) L y1 y2 , x [a, b], y1 , y2 [ y ( x) , y ( x) ]
4、 迭代格式的构造
(1) 构造思想：将连续的微分方程及初值条件离散为线性方程 组加以求解。由于离散化的出发点不同，产生出各种不同的数 值方法。基本方法有：有限差分法（数值微分）、有限体积法 （数值积分）、有限元法（函数插值）等等。 (2) 一般构造方法： 离散点函数值集合 + 线性组合结构 → 近似公式
可用来检验近似解的准确程度。 从上表最后一列，我们看到取步长 h 0.1
进行计算，数值解已达到了一定的精度。
3、 欧拉公式的改进： 隐式欧拉法 /* implicit Euler method */
向后差商近似导数
y( x1 ) y( x1 ) y( x0 ) h
y ( x1 ) y 0 h f ( x1 , y ( x1 ))
y2 y1 f ( x1 , y1 ) x2 x1
·····
pk
pn pn1
p0
Euler格式
yn1 yn f ( xn , yn ) xn1 xn
yn1 yn hf ( xn , yn )
x0 x1 xk xn xn 1
18世纪最杰出的数学家之一，13岁 时入读巴塞尔大学，15岁大学毕业， 16岁获得硕士学位。 1727年-1741年（20岁-34岁）在彼
dy f ( x, y) dx y ( a ) y0 x [a , b ]
只要 f ( x, y) 在a, b R1上连续, 且关于 y 满足 Lipschitz 条件， 即存在与 x, y 无关的常数 L 使
| f ( x, y1 ) f ( x, y2 ) | L | y1 y2 |
的方法称为微分方程的数值解法。
y1 ,
, yn 称为微分方程的数值解。
称节点间距 hi xi 1 xi (i 0, ... , n 1) 为步长，
通常采用等距节点，即取 hi = h (常数)。
3、相关定义 记 D {( x, y) a x b, y( x) y y( x) } 称 f ( x, y) 在区域D上对 y 满足Lipschitz条件是指:
得堡科学院从事研究工作，在分析学、
数论、力学方面均有出色成就，并应 俄国政府要求，解决了不少地图学、 造船业等实际问题。 24岁晋升物理学教授。
1735年（28岁）右眼失明。
1741年 - 1766（34岁-59岁）任德国科学院物理数学所所
长，任职25年。在行星运动、刚体运动、热力学、弹道学、人 口学、微分方程、曲面微分几何等研究领域均有开创性的工作。 1766年应沙皇礼聘重回彼得堡，在1771年（64岁）左眼失 明。 Euler是数学史上最多产的数学家，平均以每年800页的速 度写出创造性论文。他去世后，人们用35年整理出他的研究成
x0
x1
yi1 yi hf xi1, yi1 i 0, n 1
由于未知数 yi+1 同时出现在等式的两边，不能直接 得到，故称为隐式 /* implicit */ 欧拉公式，而前者 称为显式 /* explicit */ 欧拉公式。
一般先用显式计算一个初值，再迭代求解。
yk y( xk ) k 1,2,, n 工程师关注
解的光滑性 yk y( xk ) 解的振动性 * 解的周期性