跳轉到

10.3 專案三：科學計算工作流#

本節目標#

  • 掌握四段式接力：SymPy 推導公式 → lambdify 轉函數 → NumPy 批量計算 → Matplotlib 視覺化
  • 以「含空氣阻力的拋體運動」為完整範例走完全流程
  • 學會 Jupyter Notebook 的組織技巧
  • 建立可重現的專案資料夾（requirements.txt、隨機種子）
  • 實驗：把之前章節的實驗打包成 notebook

四段接力總覽#

階段 工具 產出
推導 SymPy（符號） 公式（運算式物件）
轉換 lambdify 快速 Python/NumPy 函數
計算 NumPy 陣列、大量數值結果
呈現 Matplotlib 圖形、結論

為什麼要接力？ SymPy 擅長「推出公式」但算大量數值很慢； NumPy 擅長「批量計算」但不會微分積分。 lambdify 是兩個世界的橋梁（詳見 4.1）。

完整範例：含空氣阻力的拋體運動#

第一段：SymPy 推導#

一個質量 \(m\) 的拋體，速度 \(v\) ，受線性空氣阻力 \(-bv\) （一維垂直分量先做）：

\[ m\frac{dv}{dt} = -mg - bv \]

用 SymPy 解出 \(v(t)\) 與高度 \(y(t)\) ：

import sympy as sp

t = sp.Symbol('t', nonnegative=True)
v = sp.Function('v')
m, b, g = sp.symbols('m b g', positive=True)

ode = sp.Eq(m*v(t).diff(t), -m*g - b*v(t))
sol_v = sp.dsolve(ode, v(t), ics={v(0): sp.Symbol('v0')})
print(sp.simplify(sol_v.rhs))
# v(t) = (b*v0 + g*m*exp(-b*t/m) - g*m)/b ... 化簡後

# 積分得到高度
y_expr = sp.integrate(sol_v.rhs, t)
y_expr = sp.simplify(y_expr.subs(t, t) )
print(y_expr)

再推「落地時間」與「射程」的符號式（垂直分量）：

# 最大高度：v(t) = 0 的時刻
t_top = sp.solve(sp.Eq(sol_v.rhs, 0), t)[0]
h_max = sp.simplify(sp.integrate(sol_v.rhs, (t, 0, t_top)))
print(h_max)   # 符號表示的最大高度

第二段：lambdify 轉函數#

v0_, b_ = 30.0, 0.5      # 初速 30 m/s，阻力係數 0.5
m_, g_ = 1.0, 9.8

params = {v0: v0_, b: b_, m: m_, g: g_}
v_func = sp.lambdify(t, sol_v.rhs.subs(params), 'numpy')
h_func = sp.lambdify(t, h_max.subs(params), 'numpy')

第三段：NumPy 批量計算#

import numpy as np

tt = np.linspace(0, 6, 300)
vv = v_func(tt)
print(f"最大高度（符號公式）: {float(h_max.subs(params)):.2f} m")

如果我們要比較 10 種阻力係數，lambdify 出的函數可以直接批量餵陣列—— 這正是 NumPy 的強項。

第四段：Matplotlib 視覺化#

import matplotlib.pyplot as plt

plt.figure(figsize=(8, 4))
for b_val in [0, 0.2, 0.5, 1.0]:
    p = {v0: v0_, b: b_val, m: m_, g: g_}
    f = sp.lambdify(t, sp.integrate(sol_v.rhs, t).subs(p), 'numpy')
    plt.plot(tt, f(tt), label=f'b = {b_val}')
plt.xlabel('t (s)'); plt.ylabel('y (m)')
plt.title('不同空氣阻力下的高度曲線')
plt.legend(); plt.grid(True)
plt.show()

解讀： \(b = 0\) 是理想拋物線；阻力越大，最大高度越低、 下落越慢——符號推導一次、數值比較多組，這就是接力的威力。

Jupyter Notebook 組織技巧#

  1. Markdown 混排：每段程式碼前後用 markdown 儲存格說明 「要做什麼 → 程式 → 觀察到什麼」，讓 notebook 自己會說話。
  2. %matplotlib inline：讓圖形直接嵌在輸出區（新版 Jupyter 已內建，但寫著無害）。
  3. 一個 notebook 一個問題：不要把三個專案塞在同一份。
  4. Restart & Run All：交付前從頭執行一遍，確認沒有依賴「亂序執行」的隱藏狀態。
  5. 用 display(sp.Eq(...))：把符號式以 LaTeX 排版呈現，比 print 易讀得多。

建議的 notebook 骨架：

# 標題：含空氣阻力的拋體運動
## 1. 問題與假設        (markdown)
## 2. 符號推導          (code: sympy)
## 3. 轉為數值函數      (code: lambdify)
## 4. 批量計算與比較    (code: numpy)
## 5. 視覺化            (code: matplotlib)
## 6. 結論              (markdown)

資料夾與可重現性#

一個可重現的專案長這樣：

projectile-lab/
├── notebook.ipynb        # 主工作區
├── requirements.txt      # 依賴清單
├── data/                 # 原始資料（若有）
├── output/               # 產出的圖與表
└── README.md             # 怎麼執行、結論摘要

requirements.txt（固定版本，別人才裝得出同樣環境）：

sympy==1.13
numpy==1.26
matplotlib==3.8
jupyter==1.0

隨機種子：只要程式裡有隨機性（蒙地卡羅、隨機取樣），開頭就固定種子：

import numpy as np
rng = np.random.default_rng(seed=42)   # 固定種子，結果可重現

檢查清單：換一台機器 clone 下來，pip install -r requirements.txt 後 Restart & Run All 能跑出同樣結果，才算「可重現」。

實驗：把之前任一章的實驗包成 notebook#

以 10.1 的「旋轉體體積」為例，打包步驟：

  1. 開新 notebook，第一格放說明 markdown（題目、假設、預期答案）。
  2. 第二格放 imports 與參數（集中管理，方便改）：
import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

x = sp.Symbol('x', nonnegative=True)
f = sp.sqrt(x)          # 想換題目只改這裡
a, b_ = 0, 4
  1. 推導格：V = sp.pi * sp.integrate(f**2, (x, a, b_)) 並 display(V)。
  2. 檢查格：V.evalf() 與 NumPy 取樣對照。
  3. 圖格：sp.plot(f, (x, a, b_)) 加上旋轉體的示意（圓盤半徑 \(f(x)\) ）。
  4. 結論格：用 markdown 寫下「答案、檢查方式、學到什麼」。

驗收標準：送給同學，對方能不看你本人就重跑出全部結果與結論。

小結與練習#

本節建立了 SymPy → lambdify → NumPy → Matplotlib 的四段接力工作流， 以空氣阻力拋體為完整範例演練；並整理了 Jupyter 組織技巧與 可重現性要求（requirements.txt、隨機種子、Restart & Run All）。 從下一節開始，我們要走向更廣的科學計算生態。

練習

  1. 把 10.1 的「極值應用題（菜園問題）」打包成 notebook， 要求包含符號推導、數值對照、畫圖三段。
  2. 在拋體範例中改用二次阻力 \(-bv^2\) （只對向下運動合理即可）， 推導並數值比較 \(b = 0.2, 0.5, 1.0\) 的最大高度。
  3. 為你的專案資料夾寫 README.md 與 requirements.txt， 交給同學測試能否在另一台機器重現。
  4. 用固定種子的蒙地卡羅法估計 \(\pi\) （隨機投點法）， notebook 裡畫出「投點數 vs 誤差」的收斂圖，並說明種子如何保證可重現。