# 數值分析 7.1 - 7.3

## 7.1

### Q1: 判別差分方程的性質 (Classification)
**題目**：判斷下列方程式是 **單步 (One-step) 還是 多步 (Multistep)**，以及是 **顯式 (Explicit) 還是 隱式 (Implicit)**。

**💡 觀念引導 (How to identify?)**
1. **Step (步數)**：看等號右邊用到了哪些時間點的 $w$？
* 只出現 $w_i$ (當下) $\rightarrow$ **One-step (單步)**。
* 出現了 $w_{i-1}, w_{i-2}$ (過去) $\rightarrow$ **Multistep (多步)**。
2. **Implicit (顯/隱式)**：看等號右邊有沒有出現「未來」的 $w_{i+1}$？
* 沒有 $w_{i+1}$ $\rightarrow$ **Explicit (顯式)** (可以直接算出答案)。
* 有 $w_{i+1}$ $\rightarrow$ **Implicit (隱式)** (未知數卡在函數裡，通常需用牛頓法求解)。

**詳細解答：**

* **(a)** $\frac{w_{i+1} - w_i}{h} = \frac{3}{2}f(t_i, w_i) - \frac{1}{2}f(t_{i-1}, w_{i-1})$
* **分析**：出現了 $i-1$ (多步)；右邊沒有 $w_{i+1}$ (顯式)。
* **答案**：**Multistep, Explicit**

* **(b)** $\frac{w_{i+1} - w_i}{h} = f(t_i + \frac{h}{2}, w_i + \frac{h}{2}f(t_i, w_i))$
* **分析**：雖然函數包得很複雜，但仔細看下標，只依賴 $w_i$ (單步)；右邊沒有 $w_{i+1}$ (顯式)。
* **答案**：**One-step, Explicit**

* **\(c)** $\frac{w_{i+1} - w_i}{h} = \frac{5}{12}f(t_{i+1}, w_{i+1}) + \dots - \frac{1}{12}f(t_{i-1}, w_{i-1})$
* **分析**：出現了 $i-1$ (多步)；右邊函數裡有 $w_{i+1}$ (隱式)。
* **答案**：**Multistep, Implicit**

* **(d)** $\frac{w_{i+1} - w_i}{h} = f(t_{i+1}, w_{i+1})$
* **分析**：只看前後兩步 (單步)；右邊直接就是未知數 $w_{i+1}$ (隱式)。*(註：這是 Backward Euler 法)*
* **答案**：**One-step, Implicit**

* **(e)** $\frac{w_{i+1} - 4w_i + 3w_{i-1}}{h} = -2f(t_{i-1}, w_{i-1})$
* **分析**：出現了 $i-1$ (多步)；右邊只有 $t_{i-1}$，無 $w_{i+1}$ (顯式)。
* **答案**：**Multistep, Explicit**

---

### Q5: 證明 Lipschitz 條件 (Existence and Uniqueness)
**題目**：證明 $f(t, y) = t(y^2 - y)$ 在有界區域 $D$ 上滿足 Lipschitz 條件，但在無界區域 ($y \in \mathbb{R}$) 不滿足。

**💡 觀念引導 (Lipschitz Condition)**
Lipschitz 條件是用來保證微分方程解的**唯一性**與**穩定性**。檢查方法是證明偏導數 $\frac{\partial f}{\partial y}$ 在區域內是**有界 (Bounded)** 的（即最大值不是無限大）。

**解答步驟：**

**1. 計算對 $y$ 的偏導數**
$$
\frac{\partial f}{\partial y} = \frac{\partial}{\partial y}(ty^2 - ty) = t(2y - 1)
$$

**2. 檢查 Case 1: 有界區域 $D = \{(t, y) \mid |t| \le 1, |y| \le M\}$**
我們想知道 $|\frac{\partial f}{\partial y}|$ 在這個範圍內最大可以是多少？
$$
\left| \frac{\partial f}{\partial y} \right| = |t| \cdot |2y - 1|
$$
* 因為 $|t| \le 1$。
* 因為 $|y| \le M$，根據三角不等式 $|2y - 1| \le 2|y| + 1 \le 2M + 1$。
* **Lipschitz 常數 $L$**：$L = 1 \cdot (2M + 1) = 2M + 1$。
* **結論**：因為我們找到了一個固定的數字 $L$，所以**滿足** Lipschitz 條件。

**3. 檢查 Case 2: 無界區域 $D = \{(t, y) \mid |t| \le 1, y \in \mathbb{R}\}$**
* 如果我們固定 $t=1$，導數變成 $|2y - 1|$。
* 當 $y$ 跑到無限大 ($\infty$) 時，導數 $|2y - 1|$ 也會衝向無限大。
* **結論**：找不到天花板 (Unbounded)，所以**不滿足** Lipschitz 條件。

---

## Chapter 7.2

這一節的重點是**尤拉法 (Euler's Method)** 的實作、誤差界限的計算，以及穩定性的分析。

核心公式：
$$w_{i+1} = w_i + h \cdot f(t_i, w_i)$$

---

### Q7: 實作 Euler 法與理論誤差界限
**題目**：$x' = e^t/x, \quad x(0)=1, \quad 0 \le t \le 2, \quad N=4$。
**精確解**：$x(t) = \sqrt{2e^t - 1}$。
**區域限制**：$D = \{(t, x) \mid 0 \le t \le 2, x \ge 1\}$。

#### Step 0: 識別函數與參數
* **方程式**：$x' = \frac{e^t}{x}$。
* **識別 $f$**：$f(t, x) = \frac{e^t}{x}$ (斜率函數)。
* **步長 $h$**：$\frac{2 - 0}{4} = \mathbf{0.5}$。

#### Step 1: 數值計算 (Numerical Calculation)

利用 Euler 公式 $w_{i+1} = w_i + h \cdot f(t_i, w_i)$ 依序計算，結果彙整表：

| 時間 $t$ | Euler 近似值 $w$ | 精確值 $x$ | 實際誤差 $\lvert x-w \rvert$ |
| :--- | :--- | :--- | :--- |
| 0.0 | 1.00000 | 1.00000 | 0.00000 |
| 0.5 | 1.50000 | 1.51573 | 0.01573 |
| 1.0 | 2.04957 | 2.10632 | 0.05674 |
| 1.5 | 2.71271 | 2.82196 | 0.10924 |
| 2.0 | 3.53876 | 3.71189 | 0.17313 |

#### Step 2: 理論誤差界限 (Theoretical Error Bound)
公式：
$$|y_i - w_i| \le \frac{hM}{2L} [e^{L(t_i - a)} - 1]$$

我們計算終點 ($t=2$) 的誤差上界（最壞情況）：

1. **Lipschitz 常數 $L$**：
* 微分：$\left| \frac{\partial f}{\partial x} \right| = \left| -\frac{e^t}{x^2} \right|$。
* 在區域 $D$ ($0 \le t \le 2, x \ge 1$) 內，取分子最大 ($t=2$)、分母最小 ($x=1$)。
* **$L = e^2 \approx 7.38906$**。

2. **二階導數上界 $M$**：
* 經推導 $|x''(t)|$ 為遞增函數，故最大值發生在 $t=2$。
* **$M \approx |x''(2)| \approx 0.92326$**。

3. **計算界限**：
$$\text{Bound} = \frac{0.5 \times 0.92326}{2 \times 7.38906} [e^{7.38906 \times 2} - 1] \approx \mathbf{81,779.5}$$

#### Step 3: 結論與圖形驗證
* **數值比較**：實際誤差 (**0.17313**) 遠小於理論界限 (**81,779.5**)。
* **分析**：理論公式雖然成立，但在 Lipschitz 常數 $L$ 較大時，指數項 $e^{L(b-a)}$ 會導致界限估計變得非常寬鬆 (Loose)，這在數值分析中是常見現象。

#### MATLAB 驗證程式碼 :
:::info

matlab
clc; clear; close all;

% --------------------------------------
% 1. 定義參數
% --------------------------------------
h = 0.5;
t_start = 0;
t_end = 2.0;
N = 4;

% 函數定義
f = @(t, x) exp(t) ./ x;             % x' = e^t / x
exact_sol = @(t) sqrt(2*exp(t) - 1); % 精確解

% 理論界限常數
L = exp(2); % Lipschitz Constant approx 7.389
% M calculation: max |x''(t)| at t=2
% x''(t) = (e^(2t) - e^t) / (2e^t - 1)^(1.5)
M_val = (exp(4) - exp(2)) / (2*exp(2) - 1)^1.5; % approx 0.923

% --------------------------------------
% 2. 計算 Euler 解與誤差
% --------------------------------------
t_vals = linspace(t_start, t_end, N+1);
w_vals = zeros(1, N+1);
w_vals(1) = 1; % Initial value x(0)=1
actual_error = zeros(1, N+1);
theo_bound = zeros(1, N+1);

for i = 1:N
    % Euler Step
    w_vals(i+1) = w_vals(i) + h * f(t_vals(i), w_vals(i));
end

% 計算每個點的誤差與理論界限
for i = 1:N+1
    t = t_vals(i);
    x_true = exact_sol(t);
    
    % 實際誤差
    actual_error(i) = abs(x_true - w_vals(i));
    
    % 理論誤差界限公式: (hM / 2L) * (e^(Lt) - 1)
    theo_bound(i) = (h * M_val / (2 * L)) * (exp(L * t) - 1);
end

% --------------------------------------
% 3. 輸出表格
% --------------------------------------
fprintf('---------------------------------------------------------------------\n');
fprintf('t      | Approx (w) | Exact (x)  | Actual Error | Theo. Bound\n');
fprintf('---------------------------------------------------------------------\n');
for i = 1:N+1
    fprintf('%.1f    | %.5f    | %.5f    | %.5f      | %.5f\n', ...
        t_vals(i), w_vals(i), exact_sol(t_vals(i)), actual_error(i), theo_bound(i));
end
fprintf('---------------------------------------------------------------------\n');

% --------------------------------------
% 4. 繪圖
% --------------------------------------
figure;

% 子圖 1: 解的比較
subplot(2, 1, 1);
t_smooth = linspace(0, 2, 100);
plot(t_smooth, exact_sol(t_smooth), 'b-', 'LineWidth', 2); hold on;
plot(t_vals, w_vals, 'ro--', 'LineWidth', 2, 'MarkerSize', 6);
legend('Exact Solution', 'Euler Approx', 'Location', 'best');
title('7.2 Q7: Solution Comparison');
xlabel('t'); ylabel('x');
grid on;

% 子圖 2: 誤差比較 (使用 Log Scale)
subplot(2, 1, 2);
semilogy(t_vals, theo_bound, 'k--^', 'LineWidth', 2, 'MarkerSize', 6); hold on;
semilogy(t_vals, actual_error, 'r-o', 'LineWidth', 2, 'MarkerSize', 6);
legend('Theoretical Bound (Upper Limit)', 'Actual Error', 'Location', 'best');
title('Error Comparison (Log Scale)');
xlabel('t'); ylabel('Error Magnitude (log)');
grid on;

:::
matlab驗證程式碼執行結果 :

### Q13: 初始條件敏感度 (Stability Analysis)
**題目**：$x' = 2e^t - x, \quad N=5$ ($h=1$)。比較 $x(0)=1$ 與 $x(0)=0.9$ 的結果。

#### Step 0: 識別函數 $f$
* **方程式**：$x' = 2e^t - x$。
* **識別 $f$**：$f(t, x) = 2e^t - x$。

#### Step 1: 推導通式
直接將 $h=1$ 代入 Euler 公式：
$$w_{i+1} = w_i + h \cdot f(t_i, w_i)$$
$$w_{i+1} = w_i + 1 \cdot (2e^{t_i} - w_i)$$
$$w_{i+1} = w_i + 2e^{t_i} - w_i$$
$$w_{i+1} = \mathbf{2e^{t_i}}$$

**關鍵發現**：化簡後，$w_i$ 被消掉了！這表示**下一步的預測值完全不依賴當前的 $w_i$，只取決於時間 $t_i$**。

#### Step 2: 計算比較

我們分別計算兩個 Case 的第一步 $w_1$ ($t=1$)：

* **Case 1 (正確初始值 $x(0)=1$)**:
$$w_1 = 2e^{t_0} = 2e^0 = \mathbf{2}$$

* **Case 2 (錯誤初始值 $x(0)=0.9$)**:
$$w_1 = 2e^{t_0} = 2e^0 = \mathbf{2}$$

#### 結論
儘管初始值不同 ($1$ vs $0.9$)，但在計算第一步之後，兩者的結果完全相同 ($w_1=2$)。這顯示在此特定步長 ($h=1$) 下，Euler 法對這個方程式具有**絕對穩定性**，初始誤差會瞬間消失。

---

### Q20: 驗證全域誤差階數 $O(h)$ (含 MATLAB 深度驗證)
**題目**：驗證 $x' = 1 - x + e^{2t}x^2, \quad x(0)=0, \quad 0 \le t \le 0.9$ 的 Euler 法全域誤差為 $O(h)$。

#### Step 0: 識別函數與參數
* **方程式**：$x' = 1 - x + e^{2t}x^2$。
* **識別 $f$**：$f(t, x) = 1 - x + e^{2t}x^2$。
* **目標**：驗證當步長 $h$ 減半時，全域誤差 $E$ 是否也減半 (即 $\text{Ratio} \approx 2$)。

#### 觀念引導 (Theoretical Basis)
Euler 法是 **一階方法 (First-order method)**，理論上誤差 $E \approx C \cdot h^1$。
這意味著：
> $$\text{Ratio} = \frac{E_{old}}{E_{new}} \approx \frac{C \cdot h}{C \cdot (h/2)} = 2$$

然而，由於本題精確解在 $t \approx 0.94$ 處有奇異點（函數值趨向無限大），在 $t=0.9$ 處導數極大，因此步長 $h$ 必須**足夠小**才能觀察到正確的收斂行為。

### Step 1: 數值實驗結果 (Numerical Results)

我們使用 MATLAB 進行深度驗證，將步長 $h$ 從 $0.05$ 一路減半至 $0.00002$。

**計算結果彙整表：**

| 步長 $h$ | 近似值 $w_{end}$ | 全域誤差 $E$ | 誤差比率 (Ratio) | 收斂階數 $p$ |
| :--- | :--- | :--- | :--- | :--- |
| 0.0500000 | 1.7902182 | 1.8511e+00 | - | - |
| 0.0250000 | 2.2264537 | 1.4149e+00 | 1.3083 | 0.3877 |
| 0.0125000 | 2.6608491 | 9.8049e-01 | 1.4430 | 0.5291 |
| 0.0062500 | 3.0262306 | 6.1511e-01 | 1.5940 | 0.6727 |
| 0.0031250 | 3.2867171 | 3.5463e-01 | 1.7345 | 0.7946 |
| 0.0015625 | 3.4488706 | 1.9247e-01 | 1.8425 | 0.8816 |
| 0.0007813 | 3.5407291 | 1.0061e-01 | 1.9130 | 0.9358 |
| 0.0003906 | 3.5898534 | 5.1490e-02 | 1.9541 | 0.9665 |
| 0.0001953 | 3.6152905 | 2.6053e-02 | 1.9764 | 0.9828 |
| 0.0000977 | 3.6282383 | 1.3105e-02 | 1.9880 | 0.9913 |
| 0.0000488 | 3.6347710 | 6.5724e-03 | 1.9940 | 0.9956 |
| 0.0000244 | 3.6380522 | 3.2912e-03 | 1.9970 | 0.9978 |
-------------------------------------------------------------------------

(註：比率 Ratio = $E_{old} / E_{new}$，收斂階數 $p = \log_2(\text{Ratio})$)

#### Step 2: 結果分析 (Conclusion)
1. **初期階段 (Pre-asymptotic)**：
當 $h$ 較大 ($0.05 \sim 0.006$) 時，收斂階數 $p$ 僅約 $0.3 \sim 0.6$，遠小於 1。這是因為在 $t=0.9$ 處函數變化極劇烈，Euler 法的線性近似產生了較大誤差。
2. **漸近階段 (Asymptotic)**：
隨著步長 $h$ 進一步縮小 ($h < 0.0001$)，**誤差比率 (Ratio)** 穩定趨近於 **2.0** (表中為 1.9970)，且 **收斂階數 (Order)** 趨近於 **1.0** (表中為 0.9978)。
3. **結論**：
實驗數據強烈證實：雖然函數本身具有剛性 (Stiffness) 特質，但只要步長足夠小，Euler's Method 仍然遵循 **$O(h)$** 的一階收斂特性。

#### MATLAB 驗證程式碼 :
:::info

clc; clear;

% 1. 定義函數與精確解
f = @(t, x) 1 - x + exp(2*t) * x^2;
exact_sol = @(t) exp(-t) * tan(exp(t) - 1);

t_start = 0; t_end = 0.9; x0 = 0;

% 2. 設定步長 (從 0.05 開始每次減半)
k_values = 1:12; 
h_values = 0.05 ./ (2 .^ (k_values - 1)); 
num_tests = length(h_values);

true_val = exact_sol(t_end);
results = zeros(num_tests, 5); 

fprintf('\n=========================================================================\n');
fprintf('   7.2 Q20: Euler Method (Deep Verification)\n');
fprintf('   Exact Value at t=0.9 : %.8f\n', true_val);
fprintf('=========================================================================\n');
fprintf('|  Step (h)   | Approx (w)   |  Error (E)   | Ratio (E_old/E) | Order (p) |\n');
fprintf('|-------------|--------------|--------------|-----------------|-----------|\n');

for i = 1:num_tests
    h = h_values(i);
    
    % Euler Iteration
    N = round((t_end - t_start) / h);
    t = t_start;
    w = x0;
    for k = 1:N
        w = w + h * f(t, w);
        t = t + h;
    end
    
    error = abs(true_val - w);
    
    ratio = 0; p = 0;
    if i > 1
        ratio = results(i-1, 3) / error;
        p = log2(ratio);
    end
    
    results(i, :) = [h, w, error, ratio, p];
    
    if i == 1
        fprintf('| %9.7f   | %10.7f   | %10.4e   |      -          |     -     |\n', h, w, error);
    else
        fprintf('| %9.7f   | %10.7f   | %10.4e   |    %8.4f     |  %6.4f   |\n', h, w, error, ratio, p);
    end
end
fprintf('-------------------------------------------------------------------------\n');

:::

## 7.3

这一節的核心觀念是：**尤拉法 (Euler's Method) 其實就是一階的泰勒方法**。為了提高準確度，我們利用更高階的導數 ($f', f'', \dots$) 來修正預測值。

**通用公式 ($n$ 階泰勒方法)**：
$$w_{i+1} = w_i + h \left[ f(t_i, w_i) + \frac{h}{2}f'(t_i, w_i) + \frac{h^2}{6}f''(t_i, w_i) + \dots + \frac{h^{n-1}}{n!}f^{(n-1)}(t_i, w_i) \right]$$

---

### Q4b: 三階泰勒法實作 (3rd-Order Taylor Method)

**題目**：$x' = \sin t - x, \quad \pi \le t \le 2\pi, \quad x(\pi)=1, \quad N=2$。

#### Step 0: 參數與導數推導 (The "Hard" Part)
這一步最容易出錯。對 $f(t, x)$ 微分時，**$x$ 也是 $t$ 的函數**，必須使用連鎖律：
$$\frac{d}{dt} g(t, x) = \frac{\partial g}{\partial t} + \frac{\partial g}{\partial x} \cdot x'$$

1. **步長 $h$**：
$$h = \frac{2\pi - \pi}{2} = \frac{\pi}{2} \approx 1.5708$$
2. **一階導數 ($f$)**：
$$x' = f = \sin t - x$$
3. **二階導數 ($f'$)**：
$$x'' = f' = \frac{d}{dt}(\sin t - x) = \cos t - x'$$
將 $x'$ 代換進來：
$$f' = \cos t - (\sin t - x) = {\cos t - \sin t + x}$$
4. **三階導數 ($f''$)**：
$$x''' = f'' = \frac{d}{dt}(\cos t - \sin t + x) = -\sin t - \cos t + x'$$
再次代換 $x'$：
$$f'' = -\sin t - \cos t + (\sin t - x) = {-\cos t - x}$$

#### Step 1: 第一步計算 ($t_0 = \pi \to t_1 = 3\pi/2$)

我們在點 $(\pi, 1)$ 處計算各階導數的值：
* **$f$** $= \sin(\pi) - 1 = 0 - 1 = \mathbf{-1}$
* **$f'$** $= \cos(\pi) - \sin(\pi) + 1 = -1 - 0 + 1 = \mathbf{0}$
* **$f''$** $= -\cos(\pi) - 1 = -(-1) - 1 = \mathbf{0}$

代入三階泰勒公式：
$$w_1 = w_0 + h \left[ f + \frac{h}{2}f' + \frac{h^2}{6}f'' \right]$$
$$= 1 + \frac{\pi}{2} \left[ -1 + 0 + 0 \right] = 1 - \frac{\pi}{2} \approx \mathbf{-0.57080}$$

#### Step 2: 第二步計算 ($t_1 = 3\pi/2 \to t_2 = 2\pi$)

我們在點 $(3\pi/2, \ 1-\pi/2)$ 處計算各階導數：
* **$f$** $= \sin(3\pi/2) - (1-\pi/2) = -1 - 1 + \pi/2 = \mathbf{\frac{\pi}{2} - 2} \approx -0.42920$
* **$f'$** $= \cos(3\pi/2) - \sin(3\pi/2) + (1-\pi/2) = 0 - (-1) + 1 - \pi/2 = \mathbf{2 - \frac{\pi}{2}} \approx 0.42920$
* **$f''$** $= -\cos(3\pi/2) - (1-\pi/2) = 0 - 1 + \pi/2 = \mathbf{\frac{\pi}{2} - 1} \approx 0.57080$

代入公式計算 $w_2$：
$$
\begin{aligned}
w_2 &= w_1 + h \left[ f + \frac{h}{2}f' + \frac{h^2}{6}f'' \right] \\
&\approx -0.57080 + 1.5708 \times \left[ -0.42920 + \frac{1.5708}{2}(0.42920) + \frac{1.5708^2}{6}(0.57080) \right] \\
&\approx -0.57080 + 1.5708 \times [ -0.42920 + 0.33710 + 0.23471 ] \\
&\approx -0.57080 + 1.5708 \times [ 0.14261 ] \\
&\approx \mathbf{-0.34676}
\end{aligned}
$$

---

### Q8b: 驗證三階誤差 $O(h^3)$ (含 MATLAB 驗證)

**題目**：使用數值實驗證明三階泰勒法的全域誤差是 $O(h^3)$。

#### 觀念引導
* **局部誤差 (Local Error)**：三階方法截斷了 $h^4$ 項，所以局部誤差是 $O(h^4)$。
* **全域誤差 (Global Error)**：因為總共走了 $1/h$ 步，累積誤差會降一階，變成 **$O(h^3)$**。
* **驗證方法 (Ratio Test)**：
如果步長減半 ($h \to h/2$)，誤差應該縮小為原來的 $(1/2)^3 = \mathbf{1/8}$。
也就是說，誤差比率 **Ratio $\approx 8$**。

#### 數值實驗結果 (MATLAB Code & Result)
:::info

clc; clear; close all;

% 1. 定義題目參數 
f = @(t, x) sin(t) - x;
% f' = cos(t) - x' = cos(t) - (sin(t) - x)
df1 = @(t, x) cos(t) - sin(t) + x;
% f'' = -sin(t) - cos(t) + x' = -sin(t) - cos(t) + (sin(t) - x) = -x - cos(t)
df2 = @(t, x) -x -cos(t);

% 精確解 x(t) (用來計算誤差)
% x(t) = 0.5 * (sin(t) - cos(t) + exp(pi - t))
exact_sol = @(t) 0.5 * (sin(t) - cos(t) + exp(pi - t));

t_start = pi;
t_end = 2*pi;
x0 = 1;

% ==========================================
% 2. 執行實驗 (步長逐次減半)
% ==========================================
% 測試 10 組步長：h, h/2, h/4 ...
N_start = 2; % 起始步數
num_tests = 10;
results = zeros(num_tests, 5); % [h, w_end, Error, Ratio, Order]

true_val = exact_sol(t_end);

fprintf('\n==========================================================================\n');
fprintf('   7.3 Q8b: 3rd-Order Taylor Method Convergence Verification\n');
fprintf('   Exact Value at t=2pi : %.8f\n', true_val);
fprintf('   Target Ratio: ~8.0 (2^3) | Target Order: ~3.0\n');
fprintf('==========================================================================\n');
fprintf('|  Step (h)   | Approx (w)   |  Error (E)   | Ratio (E_old/E) | Order (p) |\n');
fprintf('|-------------|--------------|--------------|-----------------|-----------|\n');

for k = 1:num_tests
    N = N_start * 2^(k-1);
    h = (t_end - t_start) / N;
    
    % --- 三階泰勒法迭代 (3rd-Order Taylor Iteration) ---
    t = t_start;
    w = x0;
    
    for i = 1:N
        % 計算各階導數值
        val_f = f(t, w);
        val_df1 = df1(t, w);
        val_df2 = df2(t, w);
        
        % Taylor 3rd Order Formula:
        % w_{i+1} = w_i + h * (f + h/2*f' + h^2/6*f'')
        w = w + h * (val_f + (h/2)*val_df1 + (h^2/6)*val_df2);
        
        t = t + h;
    end
    % -----------------------------------------------------
    
    error = abs(true_val - w);
    
    % 計算比率與階數
    ratio = 0; p = 0;
    if k > 1
        ratio = results(k-1, 3) / error; % Error_old / Error_new
        p = log2(ratio);
    end
    
    results(k, :) = [h, w, error, ratio, p];
    
    if k == 1
        fprintf('| %9.6f   | %10.7f   | %10.4e   |      -          |     -     |\n', h, w, error);
    else
        fprintf('| %9.6f   | %10.7f   | %10.4e   |    %8.4f     |  %6.4f   |\n', h, w, error, ratio, p);
    end
end
fprintf('--------------------------------------------------------------------------\n');

% ==========================================
% 3. 繪製 Log-Log 收斂圖
% ==========================================
figure;
h_vals = results(:, 1);
err_vals = results(:, 3);

loglog(h_vals, err_vals, 'bo-', 'LineWidth', 2, 'MarkerSize', 8, 'MarkerFaceColor', 'b');
hold on;

% 繪製參考線 Slope = 3 (代表 O(h^3))
ref_line = err_vals(end) * (h_vals / h_vals(end)).^3;
loglog(h_vals, ref_line, 'k--', 'LineWidth', 1.5);

grid on;
xlabel('Step Size (h)');
ylabel('Global Error (E)');
title('7.3 Q8b: Convergence of 3rd-Order Taylor Method');
legend('Actual Error', 'Reference Slope = 3 (O(h^3))', 'Location', 'best');
set(gca, 'XDir', 'reverse'); % X軸反轉，讓h由大變小

:::
我們使用與 Q4b 相同的方程式進行測試，計算終點 $t=2\pi$ 的誤差。


**結論**：
觀察上表，隨著步長 $h$ 逐漸變小，誤差比率 $\text{Ratio}$ 逐漸收斂並逼近 **8.0**。
這證實了當 $h \to 0$ 時，誤差縮小的比例符合 $2^3$，故三階泰勒法的全域誤差確實為 **$O(h^3)$**。

---

## Q10
比較 Euler Method、二階 Taylor Method 與四階 Taylor Method 在求解以下初始值問題 (IVP) 的表現：

$$
\frac{dx}{dt} = \frac{t}{x}, \quad (0 \le t \le 5), \quad x(0) = 1
$$

**已知精確解 (Exact Solution)：**
$$x(t) = \sqrt{t^2 + 1}$$

---

### 2. 數學推導 (Mathematical Derivation)

為了使用 Taylor Method，我們必須推導 $x(t)$ 的高階導數。
已知 $x' = f(t, x) = t x^{-1}$。

#### (1) 二階導數 $x''(t)$
對 $x'$ 關於 $t$ 微分（使用連鎖律）：
$$
\begin{aligned}
x'' &= \frac{d}{dt}(t x^{-1}) \\
&= 1 \cdot x^{-1} + t \cdot (-x^{-2}) \cdot x' \\
&= \frac{1}{x} - \frac{t}{x^2} \left( \frac{t}{x} \right) \\
&= \frac{1}{x} - \frac{t^2}{x^3} \\
&= \frac{x^2 - t^2}{x^3}
\end{aligned}
$$

#### (2) 三階導數 $x'''(t)$
對 $x'' = x^{-1} - t^2 x^{-3}$ 再次微分：
$$
\begin{aligned}
x''' &= \frac{d}{dt}(x^{-1} - t^2 x^{-3}) \\
&= (-x^{-2} \cdot x') - [2t \cdot x^{-3} + t^2 \cdot (-3x^{-4}) \cdot x'] \\
&= -x^{-2}\left(\frac{t}{x}\right) - 2tx^{-3} + 3t^2 x^{-4}\left(\frac{t}{x}\right) \\
&= -tx^{-3} - 2tx^{-3} + 3t^3 x^{-5} \\
&= \frac{3t^3 - 3tx^2}{x^5}
\end{aligned}
$$

#### (3) 四階導數 $x^{(4)}(t)$
對 $x''' = -3tx^{-3} + 3t^3 x^{-5}$ 再次微分：
$$
\begin{aligned}
x^{(4)} &= \frac{d}{dt}(-3tx^{-3} + 3t^3 x^{-5}) \\
&= [-3x^{-3} - 3t(-3x^{-4})x'] + [9t^2 x^{-5} + 3t^3(-5x^{-6})x'] \\
&= -3x^{-3} + 9tx^{-4}\left(\frac{t}{x}\right) + 9t^2 x^{-5} - 15t^3 x^{-6}\left(\frac{t}{x}\right) \\
&= -3x^{-3} + 9t^2 x^{-5} + 9t^2 x^{-5} - 15t^4 x^{-7} \\
&= -3x^{-3} + 18t^2 x^{-5} - 15t^4 x^{-7} \\
&= \frac{-3x^4 + 18t^2 x^2 - 15t^4}{x^7}
\end{aligned}
$$

---

### 3. 數值方法公式 (Numerical Schemes)

令 $h$ 為步長，$t_i = t_0 + i h$，$w_i \approx x(t_i)$。

#### Method 1: Euler's Method (一階)
截斷誤差為 $O(h)$。
$$w_{i+1} = w_i + h f(t_i, w_i) = w_i + h \left( \frac{t_i}{w_i} \right)$$

#### Method 2: Taylor Method of Order 2
截斷誤差為 $O(h^2)$。
$$w_{i+1} = w_i + h \left[ x'(t_i) + \frac{h}{2} x''(t_i) \right]$$
代入導數公式：
$$w_{i+1} = w_i + h \left[ \frac{t_i}{w_i} + \frac{h}{2} \left( \frac{w_i^2 - t_i^2}{w_i^3} \right) \right]$$

#### Method 3: Taylor Method of Order 4
截斷誤差為 $O(h^4)$。
$$w_{i+1} = w_i + h \left[ x' + \frac{h}{2} x'' + \frac{h^2}{6} x''' + \frac{h^3}{24} x^{(4)} \right]$$
其中各項導數值需代入 $(t_i, w_i)$ 計算。

---

### 4. 誤差比較與分析 (Error Analysis)

設定步長 $h=0.1$，計算區間為 $t \in [0, 5]$。我們關注終點 $t=5$ 時的絕對誤差 $|x(5) - w_N|$。

*(註：以下數據基於 MATLAB 計算結果)*

#### 數值結果表

| 方法 (Method) | 近似值 $w(5)$ | 精確值 $x(5) = \sqrt{26}$ | 絕對誤差 (Error) | 誤差階數 (Order) |
| :--- | :--- | :--- | :--- | :--- |
| **Euler** | 5.08564818 | 5.09901951 | **0.01337134** | $O(h)$ (低) |
| **Taylor 2nd** | 5.09950574 | 5.09901951 | **0.00048623** | $O(h^2)$ (中) |
| **Taylor 4th** | 5.09901909 | 5.09901951 | **0.00000042** | $O(h^4)$ (高) |

#### 結論 (Conclusion)
1. **Euler Method** 的誤差最大，僅準確到小數點後第 1 位。這是因為它僅利用切線斜率進行逼近，隨著 $t$ 增大，累積誤差迅速擴大。
2. **Taylor 2nd Order** 顯著改進了精度，誤差降低了約 100 倍（與 $h=0.1$ 到 $h^2=0.01$ 的比例相符）。
3. **Taylor 4th Order** 表現最佳，誤差極小 ($10^{-6}$ 級別)。雖然其單步計算量較大（需要計算到四階導數），但在追求高精度時，它是最有效率的選擇。

---
::: info

clc; clear; close all;

% === 1. 參數設定 ===
f = @(t, x) t./x;             % ODE: x' = t/x
exact = @(t) sqrt(t.^2 + 1);  % 精確解: x(t) = sqrt(t^2 + 1)

t_start = 0;
t_end = 5.0;     % 題目要求區間 [0, 5]
x0 = 1;          % 初始條件 x(0) = 1
h = 0.1;         % 步長 (可根據需要調整，如 0.05 或 0.1)
N = (t_end - t_start) / h;

% 初始化陣列以儲存結果 (用於畫圖)
T = t_start:h:t_end;
X_exact = exact(T);

% === 2. Euler Method (一階) ===
w_eu = zeros(1, N+1);
w_eu(1) = x0;
t = t_start;

for i = 1:N
    w = w_eu(i);
    w_eu(i+1) = w + h * (t/w);
    t = t + h;
end
err_eu = abs(X_exact(end) - w_eu(end));

% === 3. Taylor Method (二階) ===
w_t2 = zeros(1, N+1);
w_t2(1) = x0;
t = t_start;

for i = 1:N
    w = w_t2(i);
    % 計算導數值
    d1 = t/w;                % x'
    d2 = (w^2 - t^2)/w^3;    % x''
    
    % 更新公式: w + h(x' + h/2 * x'')
    w_t2(i+1) = w + h * (d1 + (h/2)*d2);
    t = t + h;
end
err_t2 = abs(X_exact(end) - w_t2(end));

% === 4. Taylor Method (四階) ===
w_t4 = zeros(1, N+1);
w_t4(1) = x0;
t = t_start;

for i = 1:N
    w = w_t4(i);
    
    % 計算各階導數值
    d1 = t/w;                                     % x'
    d2 = (w^2 - t^2)/w^3;                         % x''
    d3 = (3*t^3 - 3*t*w^2)/w^5;                   % x'''
    d4 = (-3*w^4 + 18*t^2*w^2 - 15*t^4)/w^7;      % x''''
    
    % 更新公式: w + h(x' + h/2*x'' + h^2/6*x''' + h^3/24*x'''')
    term1 = d1;
    term2 = (h/2) * d2;
    term3 = (h^2/6) * d3;
    term4 = (h^3/24) * d4;
    
    w_t4(i+1) = w + h * (term1 + term2 + term3 + term4);
    t = t + h;
end
err_t4 = abs(X_exact(end) - w_t4(end));

% === 5. 結果輸出與比較 ===
fprintf('Comparison at t = %.1f with h = %.2f\n', t_end, h);
fprintf('--------------------------------------------------\n');
fprintf('Method          | Result      | Absolute Error\n');
fprintf('--------------------------------------------------\n');
fprintf('Exact Solution  | %.8f  | 0\n', X_exact(end));
fprintf('Euler Method    | %.8f  | %.8f\n', w_eu(end), err_eu);
fprintf('Taylor Order 2  | %.8f  | %.8f\n', w_t2(end), err_t2);
fprintf('Taylor Order 4  | %.8f  | %.8f\n', w_t4(end), err_t4);

% === 6. 繪圖 (Optional) ===
figure;
plot(T, abs(X_exact - w_eu), 'r-o', 'LineWidth', 1.5, 'DisplayName', 'Euler Error');
hold on;
plot(T, abs(X_exact - w_t2), 'g-^', 'LineWidth', 1.5, 'DisplayName', 'Taylor 2 Error');
plot(T, abs(X_exact - w_t4), 'b-s', 'LineWidth', 1.5, 'DisplayName', 'Taylor 4 Error');
xlabel('t'); ylabel('Absolute Error');
title('Error Comparison of Numerical Methods');
legend('Location', 'NorthWest');
grid on;
set(gca, 'YScale', 'log'); % 使用對數座標更能看出數量級差異

:::