考慮常微分方程
d
y
d
t
=
−
y
2
,
t
∈
[
0
,
a
]
(
2
)
{\displaystyle {\frac {dy}{dt}}=-y^{2},\ t\in [0,a]\quad \quad (2)}
初始條件是
y
(
0
)
=
1
{\displaystyle y(0)=1}
。考慮格點
t
k
=
a
k
n
{\displaystyle t_{k}=a{\frac {k}{n}}}
,0 ≤ k ≤ n,意思是,時間間隔是
Δ
t
=
a
/
n
,
{\displaystyle \Delta t=a/n,}
,且
y
k
=
y
(
t
k
)
{\displaystyle y_{k}=y(t_{k})}
。用最簡單的顯式和隱式方法將此方程式离散化,分別是「前向歐拉方法」及「後向歐拉方法」,並且比較其差異。
前向歐拉方法
用不同的積分法所得的結果
y
′
=
−
y
2
,
t
∈
[
0
,
5
]
,
y
0
=
1
{\displaystyle y'=-y^{2},\;t\in [0,5],\;y_{0}=1}
,
Δ
t
=
5
/
10
{\displaystyle \Delta t=5/10}
.
前向欧拉方法
(
d
y
d
t
)
k
≈
y
k
+
1
−
y
k
Δ
t
=
−
y
k
2
{\displaystyle \left({\frac {dy}{dt}}\right)_{k}\approx {\frac {y_{k+1}-y_{k}}{\Delta t}}=-y_{k}^{2}}
可得
y
k
+
1
=
y
k
−
Δ
t
y
k
2
(
3
)
{\displaystyle y_{k+1}=y_{k}-\Delta ty_{k}^{2}\quad \quad \quad (3)\,}
對所有
k
=
0
,
1
,
…
,
n
.
{\displaystyle k=0,1,\dots ,n.}
,這是
y
k
+
1
{\displaystyle y_{k+1}}
的顯式公式。
後向歐拉方法
用後向歐拉方法(英语:backward Euler method)
y
k
+
1
−
y
k
Δ
t
=
−
y
k
+
1
2
{\displaystyle {\frac {y_{k+1}-y_{k}}{\Delta t}}=-y_{k+1}^{2}}
可以得到
y
k
+
1
{\displaystyle y_{k+1}}
的隱式方程
y
k
+
1
+
Δ
t
y
k
+
1
2
=
y
k
{\displaystyle y_{k+1}+\Delta ty_{k+1}^{2}=y_{k}}
比較上式和公式(3),公式(3)的
y
k
+
1
{\displaystyle y_{k+1}}
可以直接求得,而此處是方程式中的未知數,需要求解。
這是一元二次方程,有一個正根和一個負根,因為其初值為正,選擇其正根,則下一步的
y
{\displaystyle y}
為
y
k
+
1
=
−
1
+
1
+
4
Δ
t
y
k
2
Δ
t
.
(
4
)
{\displaystyle y_{k+1}={\frac {-1+{\sqrt {1+4\Delta ty_{k}}}}{2\Delta t}}.\quad \quad (4)}
大部份的隱式方程中,要求解的方程會比一元二次方程複雜的多,也有可能不存在解析解,因此需要用其他求根算法(例如牛顿法)來求得數值解。
克兰克-尼科尔森方法
用克兰克-尼科尔森方法
y
k
+
1
−
y
k
Δ
t
=
−
1
2
y
k
+
1
2
−
1
2
y
k
2
{\displaystyle {\frac {y_{k+1}-y_{k}}{\Delta t}}=-{\frac {1}{2}}y_{k+1}^{2}-{\frac {1}{2}}y_{k}^{2}}
可以求得
y
k
+
1
{\displaystyle y_{k+1}}
的隱式方程
y
k
+
1
+
1
2
Δ
t
y
k
+
1
2
=
y
k
−
1
2
Δ
t
y
k
2
{\displaystyle y_{k+1}+{\frac {1}{2}}\Delta ty_{k+1}^{2}=y_{k}-{\frac {1}{2}}\Delta ty_{k}^{2}}
這可以用求根算法(例如牛顿法)來求得
y
k
+
1
{\displaystyle y_{k+1}}
的數值解。
克兰克-尼科尔森方法可以視為是通用的IMEX(Implicit-Explicit,隱式-顯式)架構。
前向-後向歐拉方法
用前向歐拉方法以及前向-後向歐拉方法,在
a
=
5
{\displaystyle a=5}
和
n
=
30
{\displaystyle n=30}
下的結果
為了應用IMEX架構,考慮另一個微分方程:
d
y
d
t
=
y
−
y
2
,
t
∈
[
0
,
a
]
(
5
)
{\displaystyle {\frac {dy}{dt}}=y-y^{2},\ t\in [0,a]\quad \quad (5)}
可以得到
(
d
y
d
t
)
k
≈
y
k
+
1
−
y
k
2
,
t
∈
[
0
,
a
]
{\displaystyle \left({\frac {dy}{dt}}\right)_{k}\approx y_{k+1}-y_{k}^{2},\ t\in [0,a]}
因此
y
k
+
1
=
y
k
(
1
−
y
k
Δ
t
)
1
−
Δ
t
(
6
)
{\displaystyle y_{k+1}={\frac {y_{k}(1-y_{k}\Delta t)}{1-\Delta t}}\quad \quad (6)}
針對
k
=
0
,
1
,
…
,
n
{\displaystyle k=0,1,\dots ,n}