常微分方程初值问题数值解法
问题
一阶常微分方程的初值问题:
y
′
=
f
(
x
,
y
)
y^{'}=f(x,y)
y′=f(x,y)
y
(
x
0
)
=
y
0
y(x_0)=y_0
y(x0)=y0
定理一
设f在区域
D
=
(
x
,
y
)
∣
a
≤
x
≤
b
,
y
∈
R
D={(x,y)|a\leq x\leq b,y\in R}
D=(x,y)∣a≤x≤b,y∈R上连续,关于y满足
∣
f
(
x
,
y
1
)
−
f
(
x
,
y
2
)
∣
≤
L
∣
y
1
−
y
2
∣
,
∀
x
0
∈
[
a
,
b
]
,
y
0
∈
R
|f(x,y_1)-f(x,y_2)|\leq L|y_1-y_2|,\forall x_0\in[a,b],y_0\in R
∣f(x,y1)−f(x,y2)∣≤L∣y1−y2∣,∀x0∈[a,b],y0∈R(这条件叫做Lipschitz条件)常微分方程的初值问题存在唯一的连续可微解
y
(
x
)
y(x)
y(x)
一、一阶常微分方程初值问题的有限差分方法与误差分析
只有一些特殊的常微分方程能求解析解,实际问题中主要靠数值解,所谓数值解法,就是寻求解y(x)在一系列离散节点
x
0
<
x
1
<
x
2
.
.
.
.
.
.
x_0<x_1<x_2......
x0<x1<x2......上的近似值
y
0
,
y
1
,
y
2
.
.
.
.
.
.
y_0,y_1,y_2......
y0,y1,y2......
整体截断误差
g
i
=
∣
y
i
−
y
(
x
i
)
∣
g_i=|y_i-y(x_i)|
gi=∣yi−y(xi)∣
其中
y
(
x
i
)
y(x_i)
y(xi)就是我们的精确解,
y
i
y_i
yi就是给定计算格式(迭代公式,初值是精确值)之后所计算出来的数值解。
局部截断误差
e
i
=
∣
y
(
x
i
+
1
)
−
w
i
+
1
∣
e_i=|y(x_{i+1})-w_{i+1}|
ei=∣y(xi+1)−wi+1∣
其中
w
i
+
1
w_{i+1}
wi+1是 初值为
y
(
x
i
)
y(x_i)
y(xi) 时迭代出的数值解
对于某种方法,局部截断误差满足
e
i
=
O
(
h
p
+
1
)
e_i=O(h^{p+1})
ei=O(hp+1),则称该方法具有p阶精度。
二、向前Euler法及误差分析
1.向前Euler法
先设步长相等,大小为h,对常微分方程
y
′
=
f
(
x
,
y
)
y^{'}=f(x,y)
y′=f(x,y)中的倒数用均差近似,即
y
(
x
i
+
1
)
−
y
(
x
i
)
h
≈
y
′
=
f
(
x
i
,
y
(
x
i
)
)
\frac{y(x_{i+1})-y(x_i)}{h}\approx y^{'}=f(x_i,y(x_i))
hy(xi+1)−y(xi)≈y′=f(xi,y(xi))
若初值已知,利用迭代格式
{
y
0
=
y
(
x
0
)
y
i
+
1
=
y
i
+
h
f
(
x
i
,
y
i
)
\begin{cases} y_0=y(x_0)\\ y_{i+1}=y_i+hf(x_i,y_i) \end{cases}
{y0=y(x0)yi+1=yi+hf(xi,yi)
可逐次算出
y
1
,
y
2
.
.
.
.
.
y_1,y_2.....
y1,y2.....
2.误差分析
对于精确值做泰勒展开式
y
(
x
i
+
1
)
=
y
(
x
i
)
+
h
y
′
(
x
i
)
+
h
2
2
y
′
′
(
c
)
y(x_{i+1})=y(x_i)+hy^{'}(x_i)+\frac{h^2}{2}y^{''}(c)
y(xi+1)=y(xi)+hy′(xi)+2h2y′′(c)其中
c
∈
(
x
i
,
x
i
+
1
)
c \in(x_i,x_{i+1})
c∈(xi,xi+1),其中
y
′
(
x
i
)
=
f
(
x
i
,
y
(
x
i
)
)
y^{'}(x_i)=f(x_i,y(x_i))
y′(xi)=f(xi,y(xi))
对于欧拉法的迭代格式,由局部截断误差定义可得
y
i
+
1
=
y
(
x
i
)
+
h
y
′
(
x
i
)
y_{i+1}=y(x_i)+hy^{'}(x_i)
yi+1=y(xi)+hy′(xi)
局部截断误差
e
i
=
h
2
2
y
′
′
(
c
)
e_i=\frac{h^2}{2}y^{''}(c)
ei=2h2y′′(c)
3.后退Euler法
通过利用均差近似倒数 y ′ ( x i + 1 ) y^{'}(x_{i+1}) y′(xi+1)即
y
(
x
i
+
1
)
−
y
(
x
i
)
h
≈
y
′
(
x
i
+
1
)
=
f
(
x
i
+
1
,
y
(
x
i
+
1
)
)
\frac{y(x_{i+1})-y(x_i)}{h}\approx y^{'}(x_{i+1})=f(x_{i+1},y(x_{i+1}))
hy(xi+1)−y(xi)≈y′(xi+1)=f(xi+1,y(xi+1))
因此后退欧拉法的递推公式为
y
i
+
1
=
y
i
+
h
f
(
x
i
+
1
,
y
i
+
1
)
y_{i+1}=y_i+hf(x_{i+1},y_{i+1})
yi+1=yi+hf(xi+1,yi+1)(1式)
相比欧拉法,后退欧拉法的递推公式并不能直接进行计算,而需要解方程,即此公式是隐式的,而欧拉法的递推公式是可以直接计算的,即为显式的。隐式方程通常用迭代法求解,而迭代过程的实质是逐步隐式化。
{
y
i
+
1
0
=
y
i
+
h
f
(
x
i
,
y
i
)
y
i
+
1
1
=
y
i
+
h
f
(
x
i
,
y
i
+
1
0
)
.
.
.
.
.
.
y
i
+
1
k
+
1
=
y
i
+
h
f
(
x
i
,
y
i
+
1
k
)
(
2
式
)
\begin{cases} y_{i+1}^{0}=y_i+hf(x_i,y_i)\\ y_{i+1}^{1}=y_i+hf(x_i,y_{i+1}^{0}) \\ ......\\ y_{i+1}^{k+1}=y_i+hf(x_i,y_{i+1}^{k}) (2式) \end{cases}
⎩⎪⎪⎪⎨⎪⎪⎪⎧yi+10=yi+hf(xi,yi)yi+11=yi+hf(xi,yi+10)......yi+1k+1=yi+hf(xi,yi+1k)(2式)
由于f(x,y)满足Lipschitz条件,利用(1)式减去(2)式得到
∣
y
i
+
1
k
+
1
−
y
i
+
1
∣
=
h
(
f
(
x
i
,
y
i
+
1
k
)
−
h
f
(
x
i
+
1
,
y
i
+
1
)
)
≤
h
L
∣
y
i
+
1
k
−
y
i
+
1
∣
|y_{i+1}^{k+1}-y_{i+1}|=h(f(x_i,y_{i+1}^{k})-hf(x_{i+1},y_{i+1}))\leq hL|y_{i+1}^{k}-y_{i+1}|
∣yi+1k+1−yi+1∣=h(f(xi,yi+1k)−hf(xi+1,yi+1))≤hL∣yi+1k−yi+1∣
若
h
L
<
1
hL<1
hL<1则迭代法收敛到解
y
i
+
1
y_{i+1}
yi+1。
三、改进欧拉公式
补充:梯形方法
y
i
+
1
=
y
i
+
h
2
(
f
(
x
i
+
1
,
y
i
+
1
)
+
f
(
x
i
,
y
i
)
)
y_{i+1}=y_i+\frac{h}{2}(f(x_{i+1},y_{i+1})+f(x_{i},y_{i}))
yi+1=yi+2h(f(xi+1,yi+1)+f(xi,yi))
先用欧拉法求得一个初步的近似值
y
‾
i
+
1
\overline{y}_{i+1}
yi+1,称为预报值,然后用它替代梯形法右端的
y
i
+
1
y_{i+1}
yi+1,得到左端校正值
y
i
+
1
{y}_{i+1}
yi+1,可以校正n次,但是改进欧拉法中,我们只校正一次,这样建立的预报-校正系统称为改进的欧拉格式:
{
y
‾
i
+
1
=
y
i
+
h
f
(
x
i
,
y
i
)
y
i
+
1
=
y
i
+
h
2
(
f
(
x
i
+
1
,
y
‾
i
+
1
)
+
f
(
x
i
,
y
i
)
)
\begin{cases} \overline{y}_{i+1}=y_i+hf(x_i,y_i)\\ y_{i+1}=y_i+\frac{h}{2}(f(x_{i+1},\overline{y}_{i+1})+f(x_{i},y_{i})) \end{cases}
{yi+1=yi+hf(xi,yi)yi+1=yi+2h(f(xi+1,yi+1)+f(xi,yi))
也可以写成下面一个形式:
{
y
p
=
y
n
+
h
f
(
x
n
,
y
n
)
y
c
=
y
n
+
h
f
(
x
n
,
y
p
)
y
n
+
1
=
1
2
(
y
p
+
y
c
)
\begin{cases} {y}_{p}=y_n+hf(x_n,y_n)\\ {y}_{c}=y_n+hf(x_n,y_p)\\ y_{n+1}=\frac{1}{2}(y_p+y_c) \end{cases}
⎩⎪⎨⎪⎧yp=yn+hf(xn,yn)yc=yn+hf(xn,yp)yn+1=21(yp+yc)
四、单步法局部截断误差与阶
初值问题
y
′
=
f
(
x
,
y
)
y^{'}=f(x,y)
y′=f(x,y),
y
(
x
0
)
=
y
0
y(x_0)=y_0
y(x0)=y0的单步法可用一般形式表示为
y
n
+
1
=
y
n
+
h
φ
(
x
n
,
y
n
,
y
n
+
1
,
h
)
y_{n+1}=y_{n}+h\varphi(x_n,y_n,y_{n+1},h)
yn+1=yn+hφ(xn,yn,yn+1,h)
其中
φ
(
x
n
,
y
n
,
y
n
+
1
,
h
)
\varphi(x_n,y_n,y_{n+1},h)
φ(xn,yn,yn+1,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
)
=
f
(
x
,
y
)
\varphi(x_n,y_n,h)=f(x,y)
φ(xn,yn,h)=f(x,y)
对于显式单步法的局部截断误差
T
n
+
1
=
y
(
x
n
+
1
)
−
y
(
x
n
)
−
h
φ
(
x
n
,
y
(
x
n
)
,
h
)
T_{n+1}=y(x_{n+1})-y(x_{n})-h\varphi(x_n,y(x_n),h)
Tn+1=y(xn+1)−y(xn)−hφ(xn,y(xn),h)
五、龙格—库塔方法
1.定义
对于方程 y ′ = f ( x , y ) y^{'}=f(x,y) y′=f(x,y)的等价积分形式
y ( x n + 1 ) − y ( x n ) = ∫ x n x n + 1 f ( x , y ( x ) ) d x y(x_{n+1})-y(x_{n})=\int_{x_{n}}^{x_{n+1}} {f(x,y(x))} \ dx y(xn+1)−y(xn)=∫xnxn+1f(x,y(x)) dx
方程右端可用求积公式表示为
∫
x
n
x
n
+
1
f
(
x
,
y
(
x
)
)
d
x
≈
h
∑
i
=
1
r
c
i
f
(
x
n
+
λ
i
h
,
y
(
x
n
+
λ
i
h
)
)
\int_{x_{n}}^{x_{n+1}} {f(x,y(x))} \ dx \approx h\sum_{i=1}^rc_if(x_n+\lambda_ih,y(x_n+\lambda_ih))
∫xnxn+1f(x,y(x)) dx≈h∑i=1rcif(xn+λih,y(xn+λih))
一般来说点数r越多精度越高,上式右端相当于增量函数
φ
(
x
n
,
y
n
,
h
)
\varphi(x_n,y_n,h)
φ(xn,yn,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
n
c
i
K
i
\varphi(x_n,y_n,h)=\sum_{i=1}^n c_iK_i
φ(xn,yn,h)=∑i=1nciKi
K
1
=
f
(
x
n
,
y
n
)
K_1=f(x_n,y_n)
K1=f(xn,yn)
K
2
=
f
(
x
n
+
λ
2
h
,
y
n
+
h
u
21
K
1
)
,
λ
2
=
h
K_2=f(x_n+\lambda_2h,y_n+hu_{21}K_1),\lambda_2=h
K2=f(xn+λ2h,yn+hu21K1),λ2=h
…
K
i
=
f
(
x
n
+
λ
i
h
,
y
n
+
h
×
∑
j
=
1
i
−
1
u
i
j
K
j
)
,
λ
i
=
∑
j
=
1
i
−
1
u
i
j
K_i=f(x_n+\lambda_ih,y_n+h\times \sum_{j=1}^{i-1} u_{ij}K_j),\lambda_i=\sum_{j=1}^{i-1} u_{ij}
Ki=f(xn+λih,yn+h×∑j=1i−1uijKj),λi=∑j=1i−1uij
将此方法称为龙格库r级塔显示方法
2.常用的龙格库塔方法
常用二阶显式龙格库塔方法
y
n
+
1
=
y
n
+
h
K
2
y_{n+1}=y_n+hK_2
yn+1=yn+hK2
K
1
=
f
(
x
n
,
y
n
)
K_1=f(x_n,y_n)
K1=f(xn,yn)
K
2
=
f
(
x
n
+
h
2
,
y
n
+
h
2
K
1
)
K_2=f(x_n+\frac{h}{2},y_n+\frac{h}{2}K_1)
K2=f(xn+2h,yn+2hK1)
常用三阶
y
n
+
1
=
y
n
+
h
6
(
K
1
+
c
2
K
2
)
y_{n+1}=y_n+\frac{h}{6}(K_1+c2K_2)
yn+1=yn+6h(K1+c2K2)
K
1
=
f
(
x
n
,
y
n
)
K_1=f(x_n,y_n)
K1=f(xn,yn)
K
2
=
f
(
x
n
+
h
2
,
y
n
+
h
2
K
1
)
K_2=f(x_n+\frac{h}{2},y_n+\frac{h}{2}K_1)
K2=f(xn+2h,yn+2hK1)
K
3
=
f
(
x
n
+
h
,
y
n
−
h
K
1
+
2
h
K
2
)
K_3=f(x_n+h,y_n-hK_1+2hK_2)
K3=f(xn+h,yn−hK1+2hK2)
常用四阶(记下来,要考,哈哈哈,matlab计算内核所在)
y
n
+
1
=
y
n
+
h
6
(
K
1
+
2
K
2
+
2
K
3
+
K
4
)
y_{n+1}=y_n+\frac{h}{6}(K_1+2K_2+2K_3+K_4)
yn+1=yn+6h(K1+2K2+2K3+K4)
K
1
=
f
(
x
n
,
y
n
)
K_1=f(x_n,y_n)
K1=f(xn,yn)
K
2
=
f
(
x
n
+
h
2
,
y
n
+
h
2
K
1
)
K_2=f(x_n+\frac{h}{2},y_n+\frac{h}{2}K_1)
K2=f(xn+2h,yn+2hK1)
K
3
=
f
(
x
n
+
h
2
,
y
n
+
h
2
K
2
)
K_3=f(x_n+\frac{h}{2},y_n+\frac{h}{2}K_2)
K3=f(xn+2h,yn+2hK2)
K
4
=
f
(
x
n
+
h
,
y
n
+
h
K
3
)
K_4=f(x_n+h,y_n+hK_3)
K4=f(xn+h,yn+hK3)
六、单步法的收敛性与稳定性
1.收敛性与相容性
定义:若一种数值方法(如单步法)对于固定的
x
n
=
x
0
+
n
h
x_n=x_0+nh
xn=x0+nh,当
h
→
0
h\rightarrow 0
h→0 ,有
y
n
=
y
(
x
n
)
y_n=y(x_n)
yn=y(xn),则称该方法收敛。
定理:假设单步法具有p阶精度,且增量函数
φ
(
x
,
y
,
h
)
\varphi(x,y,h)
φ(x,y,h)满足利普希茨条件
∣
φ
(
x
,
y
,
h
)
−
φ
(
x
,
y
‾
,
h
)
∣
≤
L
φ
∣
y
−
y
‾
∣
|\varphi(x,y,h)-\varphi(x,\overline{y},h)|\leq L_{\varphi}|y-\overline{y}|
∣φ(x,y,h)−φ(x,y,h)∣≤Lφ∣y−y∣
又设初值准确,则其整体截断误差
y
(
x
n
)
−
y
n
=
O
(
h
p
)
y(x_n)-y_n=O(h^p)
y(xn)−yn=O(hp)
定义:若单步法
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
,
y
,
0
)
=
f
(
x
,
y
)
\varphi(x,y,0)=f(x,y)
φ(x,y,0)=f(x,y),则称单步法与初值问题相容(前述的单步法是用了泰勒展开式略去高阶项近似,有用求积公式离散得到的数值方法,相容性指数值方法逼近微分方程,在
h
→
0
h\rightarrow0
h→0时,可得到
y
′
(
x
)
=
f
(
x
,
y
)
y^{'}(x)=f(x,y)
y′(x)=f(x,y))
2.绝对稳定性与绝对稳定域
定义:若一种数值方法在节点值
y
n
y_n
yn上大小扰动为
δ
\delta
δ,于以后的节点值
y
m
(
m
>
n
)
y_m(m>n)
ym(m>n)上产生的偏差均不超过
δ
\delta
δ,则称该方法是稳定的
定义:对于模型方程
y
′
=
λ
y
y^{'}=\lambda y
y′=λy,使用单步法解模型方程
y
n
+
1
=
E
(
λ
h
)
y
n
y_{n+1}=E(\lambda h)y_n
yn+1=E(λh)yn,若得到
∣
E
(
h
λ
)
∣
<
1
|E(h\lambda )|<1
∣E(hλ)∣<1,则称该方法是绝对稳定的(以此可以来确定步长)。
常用:对于四阶龙格库塔方法:
−
2.78
<
h
λ
<
0
-2.78<h\lambda<0
−2.78<hλ<0,即
0
<
h
<
−
2.78
/
λ
0<h<-2.78/\lambda
0<h<−2.78/λ