1977 年 SPH 流體粒子法：無網格的天體漂流#

報案人：星體碰撞與自由表面。案情：網格一遇大變形就扭曲自殺。偵探 Lucy、Gingold、Monaghan 丟掉網格，把流體切成會呼吸的粒子。

案件檔案	內容
發生時間	1977 年
發生地點	英國與澳洲雙線作戰
報案人	天體物理、潰壩海浪工程
偵探	Leon Lucy、Gingold、Monaghan
兇器	光滑核 \(W\) 、粒子求和、無網格
結論	SPH 誕生，天體與自由表面可算

案發現場#

1977 年的流體模擬是一間被網格勒死的密室。

有限差分靠網格吃飯，一遇星體撕裂、波浪破碎就扭曲纏結，重劃網格貴到破產。天體案子最慘：兩顆原恆星相撞，沒有容器沒有邊界，網格不知鋪在哪裡。

另有兩派僵持。Lagrangian 主張跟著流體走， Eulerian 主張原地看流量。有限元素在小變形大勝，卻在自由表面前卻步。

Lucy 想算星體分裂，Gingold 與 Monaghan 想算非球對稱星，不約而同想到：流體即一堆帶質量粒子，各揹密度壓力，用光滑核感應，何須網格。

偵查過程#

偵查從核近似開始，任意場 \(f\) 的光滑估計為：

\[ \langle f\rangle = \int fWdx \]

其中 \(W\) 為光滑核函數， \(dx\) 為體積元， \(\langle f\rangle\) 為核近似場。核取高斯或三次樣條，具歸一對稱緊支撐，半徑記為 \(h\) 。

連續積分化為粒子求和，密度估計最先落網：

\[ \rho_i = \sum_j m_j W_{ij} \]

其中 \(\rho_i\) 為粒子密度， \(m_j\) 為鄰粒子質量， \(W_{ij}\) 為兩粒子間核值。

動量方程化為壓力推擠：

\[ \frac{dv_i}{dt} = -\sum_j m_j \left(\frac{p_i}{\rho_i^2} + \frac{p_j}{\rho_j^2}\right)\nabla_i W_{ij} \]

其中 \(v_i\) 為粒子速度， \(p_i\) 為粒子壓力， \(\rho_i\) 為粒子密度。

核函數	支撐域	優點	缺點
高斯核	無限速衰	光滑無比	需截斷
三次樣條	2 倍 \(h\)	緊支撐標準	導數較抖
Quintic 核	3 倍 \(h\)	高精度低噪	鄰居暴增

關鍵有三：無網格即無扭曲，自由表面自動現形；對稱壓力項保證動量守恆；代價在鄰居搜索，需鏈表或樹加速。人工黏性用於捕捉激波，太黏抹平真相，太稀振盪殺人。

結案報告#

1977 年三篇論文一出，天體模擬率先淪陷，星系碰撞與恆星形成全用 SPH 重演，後來潰壩海嘯血液流動紛紛歸案。

遺產有三：無網格正名，粒子也能微分；多物理耦合，自重輻射磁流體皆可掛粒子身上；工程落地，自由表面本是網格墳場，如今是動畫海浪常客。兇手是對拓撲的執念。

證據與工具#

關鍵證物一：密度求和虛擬碼。

# x 為位置，m 為質量，h 為光滑長度
for i in particles:
    rho[i] = 0.0
    for j in neighbors(i, h):
        rho[i] += m[j] * kernel(norm(x[i]-x[j]), h)

關鍵證物二：參數選擇。

參數	符號	典型值	影響
光滑長度	\(h\)	間距 1.2 倍	精度鄰居數
鄰居數	\(N\)	30 至 80 個	噪聲代價
人工黏性	\(alpha\)	約為 1	激波穩定

延伸閱讀伏筆：SPH 粒子哲學將呼應分子動力學。

補充：程式實作#

網格自縊身亡，偵探把流體切成 64 粒會呼吸的證人，用光滑核逐一盤問密度。

對應程式：1977-sph_density.py

本程式重演本文核心密度估計 \(rho_i = sum_j m_j W_{ij}\) ，核取一維三次樣條，支撐半徑為 \(2h\) 。 64 粒子按正弦真值 \(rho(x) = 1 + A sin(2pi x)\) 以分位數反演佈置，參數取 \(A = 0.2\) 與 \(h = 0.03\) 。每粒質量 \(m = 1/64\) ，週期域上用最小映像算距離，緊支撐外一律視為陌生人、不予理會。驗收標準是 SPH 估計對真值的 L2 誤差小於 0.05，通過才算無網格微分出師。光滑長度 \(h\) 正是本案的鬆緊箍：太大抹平正弦波，太小只剩粒子噪聲。

執行指令：

python3 _code/1977-sph_density.py

實測口供（種子固定為 0，分層分位數佈置）：

N=64 A=0.2 h=0.03 m=0.015625
rho_true mean=1.0200 est mean=1.0206
max|err|=0.0040
VERIFY L2=0.00177 rel=0.00173 (<0.05? True)
PASS

解讀要點：

L2 誤差僅 0.00177，相對誤差 0.00173，遠低於 0.05 的結案門檻。
估計均值 1.0206 對上真值均值 1.0200，最大單點誤差僅 0.0040，波形幾乎重合。
鄰居數由 \(h = 0.03\) 自動決定，64 粒下每粒約握有數個鄰居，正是精度與代價的甜蜜點。
全向量化成對距離一次算完，免去鏈表搜尋，適合當庭快速示範。

讀者可改光滑長度 \(h\) （如改為 0.06 或 0.015）或振幅 \(A\) ，重跑看 L2 誤差與最大誤差如何一升一降。

完整程式如下：

# -*- coding: utf-8 -*-
"""1977 SPH 密度估計：一維三次樣條核
對應 wiki：Lucy (1977) / Gingold-Monaghan SPH；一維三次樣條核做密度估計。
64 粒子按正弦密度 rho(x)=1+A*sin(2πx), A=0.2 分層取樣（分位數），週期域 [0,1)。
m=1/64，h=0.03，週期最小映像。驗證 SPH vs 真值 L2 誤差 < 0.05。只用 numpy，固定種子。
"""
import numpy as np

np.random.seed(0)

N = 64
A = 0.2
h = 0.03
m = 1.0 / N


def rho_true(x):
    return 1.0 + A * np.sin(2.0 * np.pi * x)


def cdf(x):
    # C(x)=x + A*(1-cos(2πx))/(2π)
    return x + A * (1.0 - np.cos(2.0 * np.pi * x)) / (2.0 * np.pi)


# 分層分位數反演 CDF（單調）求粒子位置
targets = (np.arange(N) + 0.5) / N
pos = np.empty(N)
for i, t in enumerate(targets):
    lo, hi = 0.0, 1.0
    for _ in range(60):
        mid = 0.5 * (lo + hi)
        if cdf(mid) < t:
            lo = mid
        else:
            hi = mid
    pos[i] = 0.5 * (lo + hi)
pos = np.sort(pos)


def W1d(r, h):
    # 一維三次樣條核（含 1/h 歸一）
    q = np.abs(r) / h
    w = np.zeros_like(q)
    m1 = q < 1.0
    m2 = (~m1) & (q < 2.0)
    w[m1] = (2.0 / 3.0 - q[m1] ** 2 + 0.5 * q[m1] ** 3) / h
    w[m2] = ((2.0 - q[m2]) ** 3 / 6.0) / h
    return w


# 週期最小映像成對距離 (N,N)，向量化
dx = pos[:, None] - pos[None, :]
dx = dx - np.round(dx)  # 映到 [-0.5, 0.5)
K = W1d(dx, h)
rho_est = m * np.sum(K, axis=1)
rt = rho_true(pos)
err_l2 = float(np.sqrt(np.mean((rho_est - rt) ** 2)))
err_rel = float(err_l2 / np.mean(rt))
print(f"N={N} A={A} h={h} m={m:.6f}")
print(f"rho_true mean={np.mean(rt):.4f} est mean={np.mean(rho_est):.4f}")
print(f"max|err|={np.max(np.abs(rho_est - rt)):.4f}")
print(f"VERIFY L2={err_l2:.5f} rel={err_rel:.5f} (<0.05? {err_l2 < 0.05})")
assert err_l2 < 0.05, f"SPH L2 error too large: {err_l2}"
print("PASS")