[隨筆] 關於朋友

今天走在路上特別想到的一個問題：人與人之間如何才會變成朋友？或許彼此之間有相同興趣，相同喜好，或者有某些共同人格特點。我認為過程大概是這樣的：在某個當下，有那麼一個人出現在某個特定情境，你們彼此共同經歷了一些事情，也許過程有苦有樂，然後最後很幸運的，他成為了你的朋友。

但是，有沒有可能 同樣 的一個人，我們讓他在不同時間與不同場景出現，也許 你 與 這個同一個人 就不再是朋友，甚至你們根本不可能是朋友。對嗎？ 因為人會變，一個人所經歷的人事物 會 漸漸改變 他看事情的角度，喜好，思想以及習慣。

如果在不同時間點相遇的你們已經沒有了共同興趣，或者共同喜好，那麼要變成彼此的朋友可能就變得困難。當我越這樣想的時候，越發現人與人之間能成為彼此朋友，或者能成為一生的朋友實在是一件非常奇妙的事情，因為這大概只能發生在對的時間與對的人才有可能。

願我們都能 好好珍惜彼此身邊的朋友。

[數值方法] Forward Euler's Method 與 其改良型 求解 ODE

給定 初始值問題(Initial Value Problem, IVP)
\[
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 來獲取更加精確或者穩定的數值解，但未免離題，在此不做贅述。