计算方法-常微分方程

第八章 常微分方程数值解
/* Numerical Methods for Ordinary Differential Equations */
常微分方程是常用的数学模型。 比如，描述人口增长的著名
人口模型：
y ry ( x x0 )
y(
x0
)
y0
y 为人口数量，r为人口增长率， y0为时刻x0的人口数量。
§1 欧拉方法 /* Euler’s Method */
➢ 欧拉公式：
对于xi∈[a,b], 有 y( xi ) f ( xi , y( xi ))
xi
xi+
1
向前差商近似导数
y( xi )
y( xi1) h
y( xi )
则有
y( xi1 ) y( xi ) hf ( xi , y( xi ))
Ri 的主项
/* leading term */
➢ 欧拉公式的改进：
隐式欧拉法 /* implicit Euler method */
向后差商近似导数
y( xi1 )
y( xi1) h
y( xi )
xi
y(xi+1)≈yi+hf(xi+1, y(xi+1))
yi1 yi h f ( xi1 , yi1 ) (i 0, ... , n 1)
7.1 Euler’s Method xi+1
由于未知数 yi+1 同时出现在等式的两边，不能直接得到，故 称为隐式 /* implicit */ 欧拉公式，而前者称为显式 /* explicit */ 欧拉公式。
隐式欧拉法的局部截断误差：
Ri
y( xi1 ) yi1
h2 2
y(xi ) O(h3 )
要计算出解函数 y(x) 在一系列节点 a = x0< x1<…< xn= b 处的近似值 yi y( xi ) (i 1, ... , n) 节点间距hi=xi+1-xi称为步长，当 hi = h 为常数时称为等步长。
利用y0求节点x1处的近似值y1, 再从y1来求出y2…, 直至 求出所有的yn. 称之为步进法.
即隐式欧拉公式具有 1 阶精度。
7.1 Euler’s Method
梯形公式 /* trapezoid formula */ — 显、隐式两种算法的平均
yi1 yi h f ( xi , yi ) wenku.baidu.comi 0, ... , n 1)
Ri
h2 2
y(xi ) O(h3 )
yi1 yi h f ( xi1 , yi1 ) (i 0, ... , n 1)
2xi yi
)
取步长分别为h=0.1和h=0.05, 进行计算, 结果有
7.1 Euler’s Method
k
xi
0 0.0
1 0.1
2 0.2
3 0.3
4 0.4
5 0.5
.
.
.
.
y(xi)
1.0 1.0954 1.1832 1.2649 1.3416 1.4142
. .
yi (h=0.1)
1.0 1.1000 1.1918 1.2774 1.3582 1.4351
dy f ( x, y) x [a, b] dx y(a) y0
只要 f (x, y) 在[a, b] R1 上连续，且关于 y 满足 Lipschitz 条 件，即存在与 x, y 无关的常数 L 使 | f (x, y1) f (x, y2) | L | y1 y2 | 对任意定义在 [a, b] 上的 y1(x) 和 y2(x) 都成立，则上述IVP存 在唯一解。
error */。
定义 若某算法的局部截断误差为O(hp+1)，则称该算法有p 阶 精度。
欧拉法的局部截断误差：
Ri
y(xi1)
yi1
[ y(xi ) hy(xi )
h2 2
y(xi ) O(h3)][ yi
hf (xi ,
yi )]
h2 2
y( xi ) O(h3 )
欧拉法具有 1 阶精度。
将上式中的函数值y(xi)都用近似值yi来表示, 则有数值计算 格式
yi1 yi h f ( xi , yi ) (i 0, ... , n 1)
例7.1 求解初值问题
y
y
2x y
,
y( 0 ) 1
0 x1
解: 该方程的解析解是y=(1+2x)1/2. 欧拉格式是
yi1
yi
h( yi
)
f ( xi1,
yi1 )]
yi1
yi
h 2
f ( xi , yi ) f
xi1, yi h f ( xi , yi )
(i 0, ..., n 1)
或者
y
p
yi
hf (xi ,
yi )
yc yi hf ( xi1, y p )
yi
1
( yp
yc ) /
2
7.1 Euler’s Method
例7.2 求解初值问题
y
y
2x y
,
y( 0 ) 1
该方程有解析解y(x)=y0er(x-xo)。 求常微分方程解析解的方法有多种多样，但是利用这些方法， 我们只能对极少数特殊类型的常微分方程求解；科学研究和 工程技术上的大量常微分方程的求解需借助于数值计算方法。
考虑一阶常微分方程的初值问题 /* Initial-Value Problem */:
. .
yi (h=0.05)
1.0 1.0977 1.1876 1.2713 1.3502 1.4174
. .
7.1 Euler’s Method
局部截断误差：
定义 在假设 yi = y(xi)，即第 i 步计算是精确的前提下，考虑 的截断误差 Ri = y(xi+1) yi+1 称为局部截断误差 /* local truncation
改进欧拉法 /* modified Euler’s method */
7.1 Euler’s Method
Step 1: 先用显式欧拉公式作预测，算出 yi1 yi h f ( xi , yi )
Step 2: 再将 yi1 代入隐式梯形公式的右边作校正，得到
yi 1
yi
h 2
[
f
(
x
i
,
yi
Ri
h2 2
y(xi ) O(h3)
yi 1
yi
h[ 2
f (xi ,
yi )
f ( xi1,
yi1 )]
(i 0, ..., n 1)
Ri O(h3 )
注：的确有局部截断误差 Ri y( xi1 ) yi1 O(h3 ) ， 即梯形公式具有2 阶精度，比欧拉方法有了进步。 但注意到该公式是隐式公式，不便于实际计算。