# Lecture 9 – ODE Toolbox, FFT；實驗七 水火箭

## 課程目標
- 學會使用 MATLAB 的 **ODE toolbox** 解常微分方程式
- 了解 **快速傅立葉轉換 (FFT)** 的基本概念與應用
- 進行 **水火箭實驗**，比較實驗與數值模擬結果

---

## 一、ODE Toolbox in MATLAB

### 1. 代數解法 (Symbolic Toolbox)
使用 `dsolve` 求解可解析的一階 ODE，例如：

syms v(t) b v0
equ = diff(v,t) == -b*v;
cond = [v(0) == v0];
v_sol(t) = dsolve(equ, cond);

% 代入數值
vi = 10; beta = 0.5;
tt = linspace(0, 3*vi/beta, 200);
v_t = subs(v_sol, {v0, b, t}, {vi, beta, tt});

plot(tt, v_t), xlabel('Time (s)'), ylabel('Velocity (m/s)')

👉 但多數實際問題 沒有解析解，需要數值方法。

---
### 2. 常用數值 Solver

MATLAB 提供多種 ODE solver，不同方法適合不同問題：

| Solver | 問題類型 | 精度 | 使用時機 |
|---------|----------|------|----------|
| **ode45** | 非剛性 | 中等 | 預設選擇，最常用 |
| **ode23** | 非剛性 | 低 | 容忍粗略誤差時 |
| **ode113** | 非剛性 | 低–高 | 高精度需求時 |
| **ode15s** | 剛性 | 低–中 | 當 ode45 無法收斂時 |
| **ode23s** | 剛性 | 低 | 適合簡單剛性問題 |
| **ode23t** | 中度剛性 | 低 | 無數值阻尼需求時 |
| **ode15i** | 全隱式 | 低 | 用於隱式方程或 DAE |

👉 一般情況先用 **ode45**，若不收斂再嘗試其他方法。

---

### 3. 範例：雨滴下墜
常微分方程式：
$$
\frac{dv}{dt} = g - b v
$$

MATLAB 程式：

tspan = [0, 10];
v0 = 0; g = 9.8; b = 0.5;

fun = @(t,y) (g - b*y);
[t,y] = ode45(fun, tspan, v0);

plot(t, y)
xlabel('Time (s)')
ylabel('Velocity (m/s)')
title('Falling Rain Drop with Air Resistance')

### 4. 二階 ODE：簡諧振子 (SHO)

方程式：
$$
\frac{d^2x}{dt^2} + \omega^2 x = 0
$$

---

#### 降階成一階系統
令速度 $v =\frac{dx}{dt}$，可寫成：
$$
\frac{dx}{dt} = v, \quad \frac{dv}{dt} = -\omega^2 x
$$

---

#### MATLAB 程式範例

Tp = 0.5;              % 振動週期
w = 2*pi/Tp;           % 角頻率
tspan = [0, 10*Tp];    % 模擬時間

% 定義方程式
fun = @(t,y) [y(2); -w^2*y(1)];

% 使用 ode45 解
[t,y] = ode45(fun, tspan, [2,0]);  % 初始條件 x=2, v=0

% 繪圖
plot(t, y(:,1),'r-', t, y(:,2)/w,'b--')
legend('x(t)','v(t)/w')
xlabel('Time (s)')
ylabel('Amplitude')
title('Simple Harmonic Oscillator')

---
### 5. Event 功能
MATLAB 的 ODE solver 提供 **Event 機制**，可以偵測特殊條件並控制模擬進行。
常見應用：
- 模擬小球落地 (高度 = 0)
- 模擬彈跳：落地後修正速度並繼續模擬
- 設定停止條件 (例如位移達到某值)

👉 Event 可以讓 ODE 模擬更接近真實物理過程。

---

## 二、傅立葉分析與 FFT

### 1. 基本概念
- 時域訊號可以分解為許多不同頻率的正弦波。
- **傅立葉轉換**：將訊號從 **時域 (time domain)** 轉換到 **頻域 (frequency domain)**。
- 頻譜 (spectrum)：顯示不同頻率成分的強度。

---

### 2. 快速傅立葉轉換 (FFT)
- FFT 是一種高效能的離散傅立葉轉換演算法。
- 頻率解析度公式：
$$
\Delta f = \frac{1}{T_f - T_i}
$$
取樣時間越長，頻率解析度越高。

MATLAB 範例：

Fs = 1000;              % 取樣頻率 (Hz)
t = 0:1/Fs:1-1/Fs;      % 時間序列
y = cos(2*pi*50*t) + cos(2*pi*120*t);  % 測試訊號

Y = fft(y);             % 做 FFT
L = length(y);
P2 = abs(Y/L);          % 雙邊頻譜
P1 = P2(1:L/2+1);       % 單邊頻譜
P1(2:end-1) = 2*P1(2:end-1);

f = Fs*(0:(L/2))/L;
plot(f,P1)
xlabel('Frequency (Hz)')
ylabel('|P1(f)|')
title('FFT Spectrum')

## 三、實驗七：水火箭 (Real Projectile)

### 實驗目的
- 觀察水火箭發射過程
- 理解「質量隨時間改變」的動力學特性
- 比較實驗結果與 MATLAB 模擬

---

### 實驗原理
水火箭的質量會隨著水噴出而改變，推進力公式為：
$$
F = m(t)\frac{dv}{dt} + v_e \frac{dm}{dt}
$$

其中：
- \( m(t) \)：火箭的瞬時質量
- \( v_e \)：噴射速度

👉 因為質量隨時間變化，沒有解析解，只能用 **數值方法 (如 ode45)** 模擬。
👉 真實情況還需要考慮空氣阻力與風阻力。

---

### 實驗步驟
1. 製作水火箭，設計尾翼以增加穩定性。
2. 改變不同條件：
- 水量 (多少比例的瓶子裝水)
- 打氣壓力 (初始壓縮空氣壓力)
- 發射角度
3. 發射並記錄：
- 最大高度
- 飛行時間
- 射程距離
4. 在 MATLAB 建立數值模型，模擬相同條件下的飛行軌跡。
5. 比較實驗與模擬的差異，並討論可能原因 (如空氣阻力、摩擦、能量損耗)。

---

### 參考資料

- R. Barrio-Perotti et al., *Theoretical and experimental analysis of the physics of water rockets*, Eur. J. Phys. **31**, 1131 (2010).
- [Wikipedia – Water rocket](https://en.wikipedia.org/wiki/Water_rocket)
- [NASA – Bottle Rocket Simulation](https://web.archive.org/web/20060207020049/http://exploration.grc.nasa.gov/education/rocket/BottleRocket/sim.htm)
- [NPL – Water Rocket Booklet (PDF)](https://www.npl.co.uk/getmedia/9562a6ed-13e3-45e6-8234-205f6a2538d6/wr_booklet_print.pdf)