\[
y'(t) = f(t,y(t)),\;\;\;\; y(0) = y_0
\] 一般而言，上述 一般的 IVP 問題並沒有辦法寫下解析解，但我們可以退而求其次詢問是否有合適的數值方法求解上述 初始值問題，最為基本的想法是利用所謂的 Forward Euler method:
Forward Euler Method:
\[
y_{n+1} = y_n + h \cdot f(t_n, y_n)
\]其中 $t_n = n h$。令終止時間 $T$ 在 $n$ 步之後到達，則 $T = n h$。
Comment:
1. 上式中 $h$ 稱為 迭代步長 (step size)。
2. $t_{n+1} = t_n + h$
3. 儘管 Forward Euler Method 非常簡便，但其數值誤差相當大，以下我們看個實際例子：
Example:
考慮 IVP $ y'=t-y $ 且 $ y(0) = 0 $ ，
(a) 試求解上述 IVP
(b) 現在令 $h := 0.1$ ，試求用 Forward Euler Method 求解 $y_2$
(Note: $y_2 = y(t_2) = y(2h)$)
(c) 比較 (a) 與 (b)
Solution (a):
首先改寫原式：
\[y'\left( t \right) = - y\left( t \right) + t
\]由上式可知此為線性常係數 ODE 其解可立即求得
\[\begin{array}{l}
y\left( t \right) = {e^{ - 1}}y\left( 0 \right) + \int_0^t {{e^{ - 1\left( {t - s} \right)}}sds} \\
\Rightarrow y\left( t \right) = 0 + \underbrace {{e^{ - t}}\int_0^t {{e^s}sds} }_{ = {e^{ - t}}\left( {{e^t}(t - 1) + 1} \right)}\\
\Rightarrow y\left( t \right) = (t - 1) + {e^{ - t}}
\end{array}\]
Solution (b):
由 $y_0 := y(0)=0$ 出發，則我們可以計算
\[\begin{array}{*{20}{l}}
{{y_1} = {y_0} + h \cdot f({t_0},{y_0})}\\
{ \Rightarrow {y_1} = 0 + \left( {0.1} \right) \cdot \left( {0 \cdot h - {y_0}} \right)}\\
{ \Rightarrow {y_1} = \left( {0.1} \right) \cdot \left( {0 - 0} \right) = 0}
\end{array}\]同理，我們接著計算 $y_2$
\[\begin{array}{l}
{y_2} = {y_1} + h \cdot f({t_1},{y_1})\\
\Rightarrow {y_2} = {y_1} + h \cdot \left( {1 \cdot h - {y_1}} \right)\\
\Rightarrow {y_2} = 0 + \left( {0.1} \right) \cdot \left( {0.1 - 0} \right)\\
\Rightarrow {y_2} = 0.01
\end{array}
\]
Solution (c)
由 (a) 可知
\[\begin{array}{l}
y\left( {2h} \right) = (2h - 1) + {e^{ - 2h}}\\
\Rightarrow y\left( {0.2} \right) =(0.2 - 1) + {e^{ - 0.2}} \approx {\rm{0.01873}}
\end{array}
\]但由 (b) 可知我們得到的 $y_2 = y(2h) =0.01$ 亦即具有誤差
\[
|y(2h) - y_{2}| = |0.01873 - 0.01| = 0.0873
\]
以下為簡單的 MATLAB code 實現 Forward Euler Method:
t0 = 0; %Initial time
y0 = 0; %Initial value
h = 0.1; %Step size
tfinal = 1; %Final time
F = @(t,y) t - y;
y = y0;
yout = y;
for t = 0: h: tfinal-h
y = y + h * F(t,y); %Euler Method
yout =[yout; y];
end
t_interval = 0: h: tfinal
plot(t_interval, yout, 'o-')
Improved Euler Method
為了改良上述誤差過大的問題，我們可採用以下 Predictor-Corrector 架構來改善 Forward Euler Method，亦即我們首先透過 Forward Euler Method 建構
Predictor Equation
\[
\tilde{y}_{n+1} := y_n + h \cdot f(t_n,y_n)
\]接著再修正上述的預測結果:
Corrector Equation:
\[
y_{n+1} = y_n + \frac{h}{2}(f(t_n,y_n) + f(t_{n+1}, \tilde{y}_{n+1}))
\]
我們以下再次使用前述例子來看看改良型 Euler 是否有得到更加精準的結果：
Example (Revisit):
考慮 IVP $y'=t-y$ 且 $y(0)=0$ ，
(a) 令 $h := 0.1$ ，試求用 Improved Euler Method 求解 $y_2$
Solution:
由 $y_0 = 0$ 開始，為了計算 $y_1$ 我們先透過 predictor 計算 $\tilde{y}_{1} $:
\[
\tilde{y}_{n+1} := y_n + h \cdot f(t_n,y_n) = 0
\]接著用Corrector Equation:
\[\begin{array}{l}
{y_1} = {y_0} + \frac{h}{2}\left( {f({t_0},{y_0}) + f({t_1},{{\tilde y}_1})} \right)\\
\Rightarrow {y_1} = 0 + \frac{h}{2}\left( {\left( {{t_0} - {y_0}} \right) + \left( {{t_1} - {{\tilde y}_1}} \right)} \right)\\
\Rightarrow {y_1} = \frac{h}{2}\left( {\left( {0 \cdot h - 0} \right) + \left( {1 \cdot h - 0} \right)} \right)\\
\Rightarrow {y_1} = \frac{h}{2}h = \frac{1}{2}{\left( {0.1} \right)^2} = 0.005
\end{array}
\]有了 $y_1$ 我們再透過 Predictor 計算 $\tilde{y}_{2} $:
\[\begin{array}{l}
{{\tilde y}_2} = {y_1} + h \cdot f({t_1},{y_1})\\
\Rightarrow {{\tilde y}_2} = 0.005 + h \cdot \left( {{t_1} - {y_1}} \right)\\
\Rightarrow {{\tilde y}_2} = 0.005 + h \cdot \left( {1 \cdot h - 0.005} \right)\\
\Rightarrow {{\tilde y}_2} \approx {\rm{0.01450}}
\end{array}\]再用 Corrector 計算 $y_2$
\[\begin{array}{l}
{y_2} = {y_1} + \frac{h}{2}(f({t_1},{y_1}) + f({t_2},{{\tilde y}_2}))\\
\Rightarrow {y_2} = 0.005 + \frac{h}{2}(\left( {{t_1} - {y_1}} \right) + \left( {{t_2} - {{\tilde y}_2}} \right))\\
\Rightarrow {y_2} = 0.005 + \frac{{0.1}}{2}(\left( {1 \cdot h - 0.005} \right) + \left( {2 \cdot h - {\rm{0.01450}}} \right))\\
\Rightarrow {{y}_2} \approx {\rm{0.0190}}
\end{array}
\]現在回憶前述例子中，我們知道 $0.01873$ 此表明此例中，improved Euler 確實大幅降低計算誤差。
Comments:
1. 實際數值應用中，大多數數值方法採用所謂 Runge-Kutta Method 來獲取更加精確或者穩定的數值解，但未免離題，在此不做贅述。
沒有留言:
張貼留言