- 一、显式欧拉法
- 二、显式欧拉法的改进
- 1. 隐式欧拉法
- 2. 梯形法
- 3. 两步欧拉法 (中点法)
- 三、预报-校正法 (改进欧拉法)
- 四、龙格 - 库塔 (Runge-Kutta) 法
- 2 阶龙格-库塔法
- 4 阶经典龙格-库塔法 (最常用)
- 一般的龙格 - 库塔法
考虑一阶常微分方程的初值问题
{
d
y
d
x
=
f
(
x
,
y
)
x
∈
[
a
,
b
]
y
(
a
)
=
y
0
{begin{cases} {frac{{dy}}{{dx}} = f(x,y)quad x in [a,b]} \ {y(a) = {y_0}} end{cases}}
{dxdy=f(x,y)x∈[a,b]y(a)=y0
只要 f ( x , y ) f(x, y) f(x,y) 在 [ a , b ] × R 1 [a, b] times mathbb{R}^{1} [a,b]×R1 上连续，且关于 y y y 满足 Lipschitz 条件，即存在与 x , y x, y x,y 无关的常数 L L L, s . t . s.t. s.t. ∣ f ( x , y 1 ) − f ( x , y 2 ) ∣ ≤ L ∣ y 1 − y 2 ∣ left|fleft(x, y_{1}right)-fleft(x, y_{2}right)right| leq Lleft|y_{1}-y_{2}right| ∣f(x,y1)−f(x,y2)∣≤L∣y1−y2∣, 对任意定义在 [ a , b ] [a, b] [a,b] 上的 y 1 ( x ) y_{1}(x) y1(x) 和 y 2 ( x ) y_{2}(x) y2(x) 成立，则上述初值问题的连续可微解 y ( x ) y(x) y(x) 在 [ a , b ] [a, b] [a,b] 上存在且唯一。
常微分方程数值解要计算出解函数
y
(
x
)
y(x)
y(x) 在一系列节点
a
=
x
0
<
x
1
<
…
<
x
n
=
b
a=x_{0} 向前差商近似导数
{
y
0
=
y
(
x
0
)
,
y
i
+
1
=
y
i
+
h
f
(
x
i
,
y
i
)
(
i
=
0
,
1
,
⋯
,
n
−
1
)
begin{cases} {y_0} = y({x_0}),\ {y_{i + 1}} = {y_i} + {h}f({x_i},{y_i}) quad (i = 0,1, cdots ,n - 1) end{cases}
{y0=y(x0),yi+1=yi+hf(xi,yi)(i=0,1,⋯,n−1) 定义 在假设
y
i
=
y
(
x
i
)
y_{i}=yleft(x_{i}right)
yi=y(xi) ，即第
i
i
i 步计算是精确的前提下，
R
i
=
y
(
x
i
+
1
)
−
y
i
+
1
R_{i}=yleft(x_{i+1}right)-y_{i+1}
Ri=y(xi+1)−yi+1 称为局部截断误差。 定义 若某算法的局赔截断误差为
O
(
h
p
+
1
)
Oleft(h^{p+1}right)
O(hp+1) ，则称该算法有
p
p
p 阶精度。 显式欧拉格式具有 1 阶精度 Python 实现: 向后差商近似导数
y
i
+
1
=
y
i
+
h
f
(
x
i
+
1
,
y
i
+
1
)
(
i
=
0
,
…
,
n
−
1
)
y_{i+1}=y_{i}+h fleft(x_{i+1}, y_{i+1}right) quad(i=0, ldots, n-1)
yi+1=yi+hf(xi+1,yi+1)(i=0,…,n−1) 隐式欧拉格式具有 1 阶精度 即显、隐式两种算法的平均
y
i
+
1
=
y
i
+
h
2
[
f
(
x
i
,
y
i
)
+
f
(
x
i
+
1
,
y
i
+
1
)
]
(
i
=
0
,
.
.
.
,
n
−
1
)
{y_{i + 1}} = {y_i} + frac{h}{2}[f({x_i},{y_i}) + f({x_{i + 1}},{y_{i + 1}})]quad (i = 0,;...;,n - 1)
yi+1=yi+2h[f(xi,yi)+f(xi+1,yi+1)](i=0,...,n−1) 梯形格式具有 2 阶精度 中心差商近似导数
y
i
+
1
=
y
i
−
1
+
2
h
f
(
x
i
,
y
i
)
(
i
=
1
,
.
.
.
,
n
−
1
)
{y_{i + 1}} = {y_{i - 1}} + 2h,f({x_i},{y_i})quad (i = 1,;...;,n - 1)
yi+1=yi−1+2hf(xi,yi)(i=1,...,n−1) 中点格式具有 2 阶精度 先用显式欧拉格式作预报, 算出
y
ˉ
i
+
1
=
y
i
+
h
f
(
x
i
,
y
i
)
bar{y}_{i+1}=y_{i}+h fleft(x_{i}, y_{i}right)
yˉi+1=yi+hf(xi,yi) 再将
y
ˉ
i
+
1
bar{y}_{i+1}
yˉi+1 代入隐式梯形格式的右边作校正，得到 即 预报-校正法具有 2 阶精度
{
y
i
+
1
=
y
i
+
h
[
λ
1
K
1
+
λ
2
K
2
]
K
1
=
f
(
x
i
,
y
i
)
K
2
=
f
(
x
i
+
p
h
,
y
i
+
p
h
K
1
)
left{begin{array}{l} y_{i+1}=y_{i}+hleft[lambda_{1} K_{1}+lambda_{2} K_{2}right] \ K_{1}=fleft(x_{i}, y_{i}right) \ K_{2}=fleft(x_{i}+p h, y_{i}+p h K_{1}right) end{array}right.
⎩⎨⎧yi+1=yi+h[λ1K1+λ2K2]K1=f(xi,yi)K2=f(xi+ph,yi+phK1)
λ
1
+
λ
2
=
1
lambda _1 + {lambda_2} = 1
λ1+λ2=1,
λ
2
p
=
1
2
{lambda_2}p = frac{1}{2}
λ2p=21, 算法格式有2阶精度
{
y
i
+
1
=
y
i
+
h
6
(
K
1
+
2
K
2
+
2
K
3
+
K
4
)
K
1
=
f
(
x
i
,
y
i
)
K
2
=
f
(
x
i
+
h
2
,
y
i
+
h
2
K
1
)
K
3
=
f
(
x
i
+
h
2
,
y
i
+
h
2
K
2
)
K
4
=
f
(
x
i
+
h
,
y
i
+
h
K
3
)
left{begin{array}{l} y_{i+1}=y_{i}+frac{h}{6}left(K_{1}+2 K_{2}+2 K_{3}+K_{4}right) \ K_{1}=fleft(x_{i}, y_{i}right) \ K_{2}=fleft(x_{i}+frac{h}{2}, y_{i}+frac{h}{2} K_{1}right) \ K_{3}=fleft(x_{i}+frac{h}{2}, y_{i}+frac{h}{2} K_{2}right) \ K_{4}=fleft(x_{i}+h, y_{i}+h K_{3}right) end{array}right.
⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧yi+1=yi+6h(K1+2K2+2K3+K4)K1=f(xi,yi)K2=f(xi+2h,yi+2hK1)K3=f(xi+2h,yi+2hK2)K4=f(xi+h,yi+hK3) Python 实现:
{
y
i
+
1
=
y
i
+
h
[
λ
1
K
1
+
λ
2
K
2
+
…
+
λ
m
K
m
]
K
1
=
f
(
x
i
,
y
i
)
K
2
=
f
(
x
i
+
α
2
h
,
y
i
+
β
21
h
K
1
)
K
3
=
f
(
x
i
+
α
3
h
,
y
i
+
β
31
h
K
1
+
β
32
h
K
2
)
⋯
⋯
K
m
=
f
(
x
i
+
α
m
h
,
y
+
β
m
1
h
K
1
+
β
m
2
h
K
2
+
…
+
β
m
m
−
1
h
K
m
−
1
)
left{begin{aligned} y_{i+1} &=y_{i}+hleft[lambda_{1} K_{1}+lambda_{2} K_{2}+ldots+lambda_{m} K_{m}right] \ K_{1} &=fleft(x_{i}, y_{i}right) \ K_{2} &=fleft(x_{i}+alpha_{2} h, y_{i}+beta_{21} h K_{1}right) \ K_{3} &=fleft(x_{i}+alpha_{3} h, y_{i}+beta_{31} h K_{1}+beta_{32} h K_{2}right) \ & cdots cdots \ K_{m} &=fleft(x_{i}+alpha_{m} h, y+beta_{m 1} h K_{1}+beta_{m 2} h K_{2}+ldots+beta_{m m-1} h K_{m-1}right) end{aligned}right.
⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧yi+1K1K2K3Km=yi+h[λ1K1+λ2K2+…+λmKm]=f(xi,yi)=f(xi+α2h,yi+β21hK1)=f(xi+α3h,yi+β31hK1+β32hK2)⋯⋯=f(xi+αmh,y+βm1hK1+βm2hK2+…+βmm−1hKm−1) 由于龙格-库塔法的导出基于泰勒展开, 故精度主要受解函数的光滑性影响. 对于光滑性不太好的解, 最好采用低阶算法而将步长
h
h
h 取小.import numpy as np
#显式欧拉法的Python实现
def Euler(f,x0,y0,h,l):
#f函数, x0,y0初值, h步长, l数量
xn, yn = x0, y0
n = 0
ns, xs, ys = [n], [x0], [y0]
while n <= l:
n += 1
xn += h
yn = yn + f(xn,yn) * h
ns.append(n)
xs.append(xn)
ys.append(yn)
return (ns,xs,ys)
二、显式欧拉法的改进
1. 隐式欧拉法
y
i
+
1
=
y
i
+
h
2
[
f
(
x
i
,
y
i
)
+
f
(
x
i
+
1
,
y
ˉ
i
+
1
)
]
y_{i+1}=y_{i}+frac{h}{2}left[fleft(x_{i}, y_{i}right)+fleft(x_{i+1}, bar{y}_{i+1}right)right]
yi+1=yi+2h[f(xi,yi)+f(xi+1,yˉi+1)]
y
i
+
1
=
y
i
+
h
2
[
f
(
x
i
,
y
i
)
+
f
(
x
i
+
1
,
y
i
+
h
f
(
x
i
,
y
i
)
)
]
(
i
=
0
,
…
,
n
−
1
)
y_{i+1}=y_{i}+frac{h}{2}left[fleft(x_{i}, y_{i}right)+fleft(x_{i+1}, y_{i}+h fleft(x_{i}, y_{i}right)right)right] quad(i=0, ldots, n-1)
yi+1=yi+2h[f(xi,yi)+f(xi+1,yi+hf(xi,yi))](i=0,…,n−1)import numpy as np
#四阶龙格库塔法的Python实现
def RK4(f,x0,y0,h,maxx):
#f函数, x0,y0初值, h步长, maxx: x最大值
xn, yn = x0, y0
n = 0
ns, xs, ys = [n], [xn], [yn]
while xn < maxx:
n += 1
xn += h
k1 = f(xn,yn)
k2 = f(xn + (h / 2),yn + (h * k1) / 2)
k3 = f(xn + (h / 2),yn + (h * k2) / 2)
k4 = f(xn + h,yn + h * k3)
yn = yn + (h / 6) * (k1 + 2 * k2 + 2 * k3 + k4)
ns.append(n)
xs.append(xn)
ys.append(yn)
return (ns,xs,ys)
一般的龙格 - 库塔法
