\qquad
常微分方程是描述连续变化的数学语言,微分方程的求解就是确定给定方程的可微方程
y
(
t
)
y(t)
y(t)。考虑一阶常微分方程的初值问题
y
′
=
f
(
x
,
y
)
,
x
∈
[
x
0
,
b
]
,
(
1
)
y
(
x
0
)
=
y
0
(
2
)
y^{'}=f(x,y),\quad x\in [x_{0},b],\quad(1)\\ y(x_{0})=y_{0}\quad(2)
y′=f(x,y),x∈[x0,b],(1)y(x0)=y0(2)
\qquad
如果存在实数
L
>
0
L>0
L>0,使得
∣
f
(
x
,
y
1
)
−
f
(
x
,
y
2
)
∣
≤
L
∣
y
1
−
y
2
∣
,
∀
y
1
,
y
2
∈
R
,
|f(x,y_{1})-f(x,y_{2})|\leq L|y_{1}-y_{2}|,\forall y_{1},y_{2}\in\mathbb{R},
∣f(x,y1)−f(x,y2)∣≤L∣y1−y2∣,∀y1,y2∈R,则称
f
f
f关于
y
y
y满足利普西茨(
L
i
p
s
c
h
i
t
e
Lipschite
Lipschite)条件,
L
L
L称为
f
f
f的利普西茨常数(简称
L
i
s
p
.
Lisp.
Lisp.常数)
\qquad
定理1 设
f
f
f在区域
D
=
{
(
x
,
y
)
∣
a
≤
x
≤
b
,
y
∈
R
}
D=\{(x,y)|a\leq x\leq b,y\in\mathbb{R}\}
D={(x,y)∣a≤x≤b,y∈R}上连续,关于
y
y
y满足利普西茨条件,则对任意
x
0
∈
[
a
,
b
]
,
y
0
∈
R
x_{0}\in [a,b],y_{0}\in\mathbb{R}
x0∈[a,b],y0∈R,常微分方程初值问题满足(1)(2)式当
x
∈
[
a
,
b
]
x\in [a,b]
x∈[a,b]时存在唯一的连续可微解
y
(
x
)
y(x)
y(x)。
\qquad
所谓数值解法,就是寻求解
y
(
x
)
y(x)
y(x)在一系列离散节点
x
1
<
x
2
<
⋯
<
x
n
<
⋯
x_{1}<x_{2}<\cdots<x_{n}<\cdots
x1<x2<⋯<xn<⋯
\quad
上的近似值
y
1
,
y
2
,
⋯
 
,
y
n
,
⋯
y_{1},y_{2},\cdots,y_{n},\cdots
y1,y2,⋯,yn,⋯相邻两节点的间距
h
n
=
x
n
+
1
−
x
n
h_{n}=x_{n+1}-x_{n}
hn=xn+1−xn称为步长。求解过程首先对常微分方程离散化,建立求数值解的递推公式。当计算
y
n
+
1
y_{n+1}
yn+1时只用到前一点
y
n
y_{n}
yn,称为单步法。用到前面
k
k
k个点的值
y
n
,
y
n
−
1
,
⋯
y_{n},y_{n-1},\cdots
yn,yn−1,⋯时称为
k
k
k步法。
\qquad
欧拉法 设做出折线(由切线相连(
x
0
,
x
1
,
⋯
x
n
x_{0},x_{1},\cdots x_{n}
x0,x1,⋯xn处作切线))的顶点
P
n
P_{n}
Pn,过
P
n
(
x
n
,
y
n
)
P_{n}(x_{n},y_{n})
Pn(xn,yn)依方向场的方向推进到
P
n
+
1
(
x
n
+
1
,
y
n
+
1
)
P_{n+1}(x_{n+1},y_{n+1})
Pn+1(xn+1,yn+1),显然两个顶点
P
n
,
P
n
+
1
P_{n},P_{n+1}
Pn,Pn+1的坐标有关系
y
n
+
1
−
y
n
x
n
+
1
−
x
n
=
f
(
x
n
,
y
n
)
\frac{y_{n+1}-y_{n}}{x_{n+1}-x_{n}}=f(x_{n},y_{n})
xn+1−xnyn+1−yn=f(xn,yn)
\qquad
对微分方程式(1)中的导数用均差近似,即
y
n
+
1
−
y
n
h
≈
y
′
(
x
n
)
=
f
(
x
n
,
y
(
x
n
)
)
\frac{y_{n+1}-y_{n}}{h}\approx y^{'}(x_{n})=f(x_{n},y(x_{n}))
hyn+1−yn≈y′(xn)=f(xn,y(xn))
\qquad
依均差近似逐次计算。为分析计算公式的精度,使用泰勒展开将
y
n
+
1
y_{n+1}
yn+1在
x
n
x_{n}
xn处展开,则有
y
(
x
n
+
1
)
=
y
(
x
n
+
h
)
=
y
(
x
n
)
+
y
′
(
x
n
)
h
+
h
2
2
y
′
′
(
ξ
n
)
,
ξ
n
∈
(
x
n
,
x
n
+
1
)
y(x_{n+1})=y(x_{n}+h)=y(x_{n})+y^{'}(x_{n})h+\frac{h^{2}}{2}y^{''}(\xi_{n}),\quad \xi_{n}\in(x_{n},x_{n+1})
y(xn+1)=y(xn+h)=y(xn)+y′(xn)h+2h2y′′(ξn),ξn∈(xn,xn+1)
\qquad
在
y
n
=
y
x
n
y_{n}=y_{x_{n}}
yn=yxn的前提下,
f
(
x
n
,
y
n
)
=
f
(
x
n
,
y
(
x
n
)
)
=
y
′
(
x
n
)
f(x_{n},y_{n})=f(x_{n},y(x_{n}))=y^{'}(x_{n})
f(xn,yn)=f(xn,y(xn))=y′(xn),于是欧拉法的误差
y
(
x
n
+
1
)
−
y
n
+
1
=
h
2
2
y
′
′
(
ξ
n
)
≈
h
2
2
y
′
′
(
x
n
)
y(x_{n+1})-y_{n+1}=\frac{h^{2}}{2}y^{''}(\xi_{n})\approx\frac{h^{2}}{2}y^{''}(x_{n})
y(xn+1)−yn+1=2h2y′′(ξn)≈2h2y′′(xn)称为此方法的局部截断误差。显式单步法的局部截断误差一般性的表示方法为:
T
n
+
1
=
y
(
x
n
+
1
)
−
y
(
x
n
)
−
h
φ
(
x
n
,
y
(
x
n
)
,
h
)
=
O
(
h
p
+
1
)
T_{n+1}=y(x_{n+1})-y(x_{n})-h\varphi(x_{n},y(x_{n}),h)=O(h^{p+1})
Tn+1=y(xn+1)−y(xn)−hφ(xn,y(xn),h)=O(hp+1)则称该方法具有
p
p
p阶精度。
\qquad
如果对微分方程从
x
n
x_{n}
xn到
x
n
+
1
x_{n+1}
xn+1积分,得
y
(
x
n
+
1
)
=
y
(
x
n
)
+
∫
x
n
x
n
+
1
f
(
t
,
y
(
t
)
)
d
t
y(x_{n+1})=y(x_{n})+\int_{x_{n}}^{x_{n+1}}f(t,y(t))dt
y(xn+1)=y(xn)+∫xnxn+1f(t,y(t))dt右端积分使用右矩形公式
h
f
(
x
n
+
1
,
y
(
x
n
+
1
)
)
hf(x_{n+1},y(x_{n+1}))
hf(xn+1,y(xn+1))近似得到另一个公式
y
n
+
1
=
y
n
+
h
f
(
x
n
+
1
,
y
n
+
1
)
y_{n+1}=y_{n}+hf(x_{n+1},y_{n+1})
yn+1=yn+hf(xn+1,yn+1)通常使用迭代法求解(实现
x
n
→
x
n
+
1
x_{n}\to x_{n+1}
xn→xn+1的运算),设用欧拉公式
y
n
+
1
(
0
)
=
y
n
+
h
f
(
x
n
,
y
n
)
y_{n+1}^{(0)}=y_{n}+hf(x_{n},y_{n})
yn+1(0)=yn+hf(xn,yn)给出迭代初值
y
n
+
1
(
0
)
y_{n+1}^{(0)}
yn+1(0)带入上式得
y
n
+
1
(
1
)
=
y
n
+
h
f
(
x
n
+
1
,
y
n
+
1
(
0
)
)
y_{n+1}^{(1)}=y_{n}+hf(x_{n+1},y^{(0)}_{n+1})
yn+1(1)=yn+hf(xn+1,yn+1(0))反复迭代得到
y
n
+
1
(
k
+
1
)
=
y
n
+
h
f
(
x
n
+
1
,
y
n
+
1
(
k
)
)
,
k
=
0
,
1
,
⋯
y_{n+1}^{(k+1)}=y_{n}+hf(x_{n+1},y^{(k)}_{n+1}),k=0,1,\cdots
yn+1(k+1)=yn+hf(xn+1,yn+1(k)),k=0,1,⋯由利普西茨条件得(
证
明
迭
代
法
收
敛
\color{#F00}{证明迭代法收敛}
证明迭代法收敛)
∣
y
n
+
1
(
k
+
1
)
−
y
n
+
1
∣
=
h
∣
f
(
x
n
+
1
,
y
n
+
1
(
k
)
)
−
f
(
x
n
+
1
,
y
n
+
1
)
∣
≤
h
L
∣
y
n
+
1
(
k
)
−
y
n
+
1
∣
|y^{(k+1)}_{n+1}-y_{n+1}|=h|f(x_{n+1},y^{(k)}_{n+1})-f(x_{n+1},y^{}_{n+1})|\leq hL|y^{(k)}_{n+1}-y_{n+1}|
∣yn+1(k+1)−yn+1∣=h∣f(xn+1,yn+1(k))−f(xn+1,yn+1)∣≤hL∣yn+1(k)−yn+1∣由此,只要
h
L
<
1
hL<1
hL<1使用迭代法就能收敛到解
y
n
+
1
y_{n+1}
yn+1。
\qquad
梯形方法(增加精度) 公式为
y
n
+
1
=
y
n
+
h
2
[
f
(
x
n
,
y
n
)
+
f
(
x
n
+
1
,
y
n
+
1
)
]
y_{n+1} = y_{n}+\frac{h}{2}[f(x_{n},y_{n})+f(x_{n+1},y_{n+1})]
yn+1=yn+2h[f(xn,yn)+f(xn+1,yn+1)],迭代公式为:
y
n
+
1
(
0
)
=
y
n
+
h
f
(
x
n
,
y
n
)
y
n
+
1
(
k
+
1
)
=
y
n
+
h
2
(
f
(
x
n
,
y
n
)
+
f
(
x
n
+
1
,
y
n
+
1
(
k
)
)
,
k
=
0
,
1
,
2
,
⋯
y_{n+1}^{(0)}=y_{n}+hf(x_{n},y_{n})\\ y_{n+1}^{(k+1)}=y_{n}+\frac{h}{2}(f(x_{n},y_{n})+f(x_{n+1},y_{n+1}^{(k)}),\quad k=0,1,2,\cdots
yn+1(0)=yn+hf(xn,yn)yn+1(k+1)=yn+2h(f(xn,yn)+f(xn+1,yn+1(k)),k=0,1,2,⋯
\qquad
改进欧拉公式(增加了欧拉公式的精度,减小梯形算法迭代带来的计算次数),建立预测-校正系统。
预测 y ˉ n + 1 = y n + h f ( x n , y n ) \bar{y}_{n+1}=y_{n}+hf(x_{n},y_{n}) yˉn+1=yn+hf(xn,yn)
校正
y
n
+
1
=
y
n
+
h
2
[
f
(
x
n
,
y
n
)
+
f
(
x
n
+
1
,
y
ˉ
n
+
1
)
]
y_{n+1}=y_{n}+\frac{h}{2}[f(x_{n},y_{n})+f(x_{n+1},\bar{y}_{n+1})]
yn+1=yn+2h[f(xn,yn)+f(xn+1,yˉn+1)]
\qquad
龙格-库塔法,由公式
y
(
x
n
+
1
)
=
y
(
x
n
)
+
∫
x
n
x
n
+
1
f
(
t
,
y
(
t
)
)
d
t
y(x_{n+1})=y(x_{n})+\int_{x_{n}}^{x_{n+1}}f(t,y(t))dt
y(xn+1)=y(xn)+∫xnxn+1f(t,y(t))dt要使近似精度提高就是提高右式的近似精度,近似公式为:
∫
x
n
x
n
+
1
f
(
t
,
y
(
t
)
)
d
t
≈
h
∑
i
=
1
r
c
i
f
(
x
n
+
λ
i
h
,
y
(
x
n
+
λ
i
h
)
)
\int_{x_{n}}^{x_{n+1}}f(t,y(t))dt\approx h\sum_{i=1}^{r}c_{i}f(x_{n}+\lambda _{i}h,y(x_{n}+\lambda _{i}h))
∫xnxn+1f(t,y(t))dt≈hi=1∑rcif(xn+λih,y(xn+λih))分割的段越多,
r
r
r越大,上式右端(相当于增量函数
φ
(
x
n
,
y
(
x
n
)
h
)
\varphi(x_{n},y(x_{n})h)
φ(xn,y(xn)h))的近似精度越高。由公式
y
n
+
1
=
y
n
+
h
φ
(
x
n
,
y
n
,
h
)
y_{n+1} = y_{n} + h\varphi(x_{n},y_{n},h)
yn+1=yn+hφ(xn,yn,h)其中,
φ
(
x
n
,
y
n
,
h
)
=
∑
i
=
1
r
c
i
K
i
K
1
=
f
(
x
n
,
y
n
)
K
i
=
f
(
x
n
+
λ
i
h
,
y
n
+
∑
j
=
1
i
−
1
μ
i
j
K
j
)
i
=
2
,
3
,
⋯
 
,
r
\varphi(x_{n},y_{n},h)=\sum_{i=1}^{r}c_{i}K_{i}\\ K_{1}=f(x_{n},y_{n})\\ K_{i}=f(x_{n}+\lambda_{i}h,y_{n}+\sum_{j=1}^{i-1}\mu_{ij}K_{j})\quad i=2,3,\cdots,r
φ(xn,yn,h)=i=1∑rciKiK1=f(xn,yn)Ki=f(xn+λih,yn+j=1∑i−1μijKj)i=2,3,⋯,r 通过给定局部截断误差精度求解未知系数。对于步长的选择,由给定的精度
ε
\varepsilon
ε和偏差
Δ
=
∣
y
n
+
1
(
h
2
)
−
y
n
+
1
(
h
)
∣
\Delta=|y_{n+1}^{(\frac{h}{2})}-y_{n+1}^{(h)}|
Δ=∣yn+1(2h)−yn+1(h)∣(
折
半
前
后
两
次
的
偏
差
\color{#F00}{折半前后两次的偏差}
折半前后两次的偏差)决定。当精度大于偏差时,减小步长(对折),反之则反。
\qquad
单步法的收敛性 假设单步法具有
p
p
p阶精度,且增量函数
φ
(
x
,
y
,
h
)
\varphi(x,y,h)
φ(x,y,h)关于
y
y
y满足利普西茨条件
∣
φ
(
x
,
y
,
h
)
−
φ
(
x
,
y
ˉ
,
h
)
∣
≤
L
φ
∣
y
−
y
ˉ
∣
|\varphi(x,y,h)-\varphi(x,\bar{y},h)|\leq L_{\varphi}|y-\bar{y}|
∣φ(x,y,h)−φ(x,yˉ,h)∣≤Lφ∣y−yˉ∣又设初值
y
0
y_{0}
y0是准确的,即
y
0
=
y
(
x
0
)
y_{0}=y(x_{0})
y0=y(x0),其整体的截断误差
y
(
x
n
)
−
y
n
=
O
(
h
p
)
y(x_{n})-y_{n}=O(h^{p})
y(xn)−yn=O(hp)稳定性 由模型方程
y
′
=
λ
y
y^{'}=\lambda y
y′=λy的欧拉公式为:
y
n
+
1
=
(
1
+
h
λ
)
y
n
y_{n+1}=(1+h\lambda)y_{n}
yn+1=(1+hλ)yn设在节点值
y
n
y_{n}
yn上有一个扰动值
ε
n
\varepsilon_{n}
εn,它的传播使节点
y
n
+
1
y_{n+1}
yn+1产生
ε
n
+
1
\varepsilon_{n+1}
εn+1的扰动值,扰动值满足
ε
n
+
1
=
(
1
+
h
λ
)
ε
n
\varepsilon_{n+1}=(1+h\lambda)\varepsilon_{n}
εn+1=(1+hλ)εn为使差分方程的解非增(解稳定),选择
h
h
h充分小,使得
∣
1
+
h
λ
∣
≤
1
|1+h\lambda|\leq1
∣1+hλ∣≤1,即
−
2
<
λ
<
0
-2<\lambda<0
−2<λ<0。