可聽出鼓的形狀嗎 ?
一一 對偶邊界元素法分析
陳正宗 · 郭柏伸 · 高怡絹
∗摘要: 「我們可以聽出鼓的形狀嗎?」, 這個問題一直是數學家、 物理學家與工程師皆 感興趣的議題。 在 1992年數學家 Gordon 等人從數學角度去驗證兩個不同的幾何形 狀, 卻有著共同的特徵頻率分佈, 來否定聽出鼓的形狀可能性。 本文主要利用對偶邊 界元素法來分析兩個不同幾何形狀的鼓皮, 藉由直接搜尋法來找出其特徵頻率; 並且 由數值實驗結果發現奇異式比超奇異式在多邊形領域更適合求解特徵值問題。 利用 對偶邊界元素法數值分析驗證兩個不同幾何形狀的鼓皮在固定端的邊界條件下, 有 著共同特徵頻率分佈, 並繪出其所對應的特徵模態圖與文獻中實驗數據或有限元素 法分析的數值結果來做比對, 結果比較均相當吻合。 因此, 由此數值結果再次否定了
「我們可以聽出鼓的形狀」 的猜想。
關鍵詞: 對偶邊界元素、 特徵頻率、 奇異式、 超奇異式。
1. 前言
從數學傳播期刊中, 我們看到一個有趣的議題 「我們可以聽出鼓的形狀嗎?」[1], 此議題也 在日常生活中常會遇到, 例如: 在炎熱的夏天, 買西瓜時會把耳朵放在西瓜旁, 並用手敲擊西瓜, 聽出西瓜的實虛(密度), 來判斷其甜度。 又從一維的觀點來看, 我們可以藉由振動頻率去判斷弦 的形狀, 也就是弦長。 進而推展到二維, 就可以去判斷鼓的形狀 [2, 3, 4]。 而本文中所探討的就 是二維的問題。 這個有趣的議題, 最早從 1910 年荷蘭物理學家 Lorentz 在 Gottingen 大學演 講時所提出來的, 其背景理論來自輻射理論。 在未被證實的報告中提到 Hilbert 認為 Lorentz 所提出的猜想在他有生之年不會被證明出來, 但過了一年後卻被自己的學生 Weyl 用 Hilbert 積分方程理論證明出我們可以聽出鼓的面積。 而在1954 年瑞典數學家 Pleijel 證明出可以聽出
∗感謝教育部與國科會大專生專題補助。
70
鼓的周長, 並且可以透過特徵頻率來決定鼓是否為圓形。 然而 Kac [5] 想嘗試證明 Pleijel 的 結果, 卻沒有成功。 此時 Kac 把鼓的外形假設為多邊形, 所證明出來的結果與 Pleijel 的驗證 結果相近。 而在 Kac 數學公式證明裡暗示這鼓有 N 個洞。 之後 1967 年 McKean 與 Singer 從 Kac 數學公式證明出我們可以聽出鼓是否有洞。 直到了 1992 年數學家 Gordon [6, 7] 等 人找到反例, 發現兩個不同的幾何形狀, 卻有著共同的特徵頻率分佈, 來否定聽出鼓的形狀可能 性。 這個有趣的議題, 讓我們體會到物理與數學內涵是如此的美妙。 本文我們將利用對偶邊界元 素法的奇異式與超奇異式再次論證兩種不同的幾何外形的鼓皮, 是否有著共同的特徵頻率分佈。
什麼是邊界元素法, 或許大家比較熟悉的是有限元素法, 因為邊界元素法發展的時間較有 限元素法晚。 這兩者主要不同乃是邊界元素法只需對問題的邊界做離散, 而不像有限元素須對 問題的領域作離散。 邊界元素法除了可以應用在聲學振動問題 [8]-[17], 也可以應用在流體力 學、 彈性波、 電磁波及土壤互制結構 · · · 等 [18]-[23], 範圍相當廣泛, 此方法比有限元素法 多出許多優勢, 讓現今許多工程師或是學者在解決工程上問題時, 提供另一個有力的模擬分析 工具。 而邊界元素法最基本的數理觀念就是工數中向量微積分與勢能理論, 透過輔助系統的基 本解和欲解系統可做兩個系統的功能互換。 以數學而言, 即所謂Green 第三恆等式, 進而建立 邊界積分方程式 [24]。 當處理一實際工程問題時, 使用邊界元素法時, 我們會先將邊界積分方程 中的邊界做切割, 切成 n 個元素, 而不像有限元素法需對領域做切割, 參見圖1 與圖 2。 將邊界 積分方程中的邊界物理量以插值近似, 再對邊界作積分, 以達到求解問題的目的。
圖 1. 有限元素法網格 — 對領域做切割 圖 2. 邊界元素法網格 — 僅對邊界做切割
2. 問題描述與對偶邊界元素法
2.1.
問題描述
鼓皮自由振動問題在頻率域裡控制方程為二維 Helmholtz 方程式, 可表示如下:
(∇2+ k2)u(x) = 0, x∈ D, (1) 其中 ∇2 為 Laplace 運算子, k 為波數, k = ω
c, 其中, ω 為角頻率, c 為波速, u(x) 為 鼓皮在 x 處的振動位移, x 為場點 (觀察點), D 為鼓皮的領域。 而問題的邊界條件為固定端
(Dirichlet 型), 表示如下:
u(x) = 0, x∈ B, (2) 其中 B 為鼓皮的邊界, 如圖 3 所示:
圖 3. 多邊形鼓皮領域
2.2.
對偶邊界元素法
根據洪宏基與陳正宗教授所提對偶邊界元素法 [25], 可以得到兩條域內點之對偶邊界積分 方程式, 表示如下:
2πu(x) = Z
B
T(s, x)u(s)dB(s) − Z
B
U(s, x)t(s)dB(s), x∈ D, (3) 2πt(x) =
Z
B
M(s, x)u(s)dB(s) − Z
B
L(s, x)t(s)dB(s), x∈ D, (4) 其中 (3) 及 (4) 式分別為奇異積分方程及超奇異積分方程, x 為場點(觀察點), s 為源點位置 向量, t(x) 為 u(x) 的法向導微量
t(x) = ∂u(x)
∂nx
, u(s) 與 t(s) 為邊界物理量, U、T 、L 與 M 為核函數, 其所造成的勢位分別對應勢論中的單層勢能、 雙層勢能及法向導微, 而在 (3) 式中的核函數 U (s, x) 為基本解 (自由空間格林函數) 並滿足以下的方程式:
(∇2+ k2)U(s, x) = 2πδ(x − s), (5) 其中 δ 為 Dirac-delta 函數, 而 U(s, x)、 T (s, x)、 L(s, x) 與 M(s, x) 之核函數, 其定義 分別如下所示:
U(s, x) = −iπH0(1)(kr)
2 , (6)
T(s, x) = ∂U(s, x)
∂ns
, (7)
L(s, x) = ∂U(s, x)
∂nx
, (8)
M(s, x) = ∂U2(s, x)
∂nx∂ns
, (9)
H0(1)(kr) 為第一類第零階 Hankel 函數, r 為場點與源點之間的距離, ns 為源點 (s) 的單位 法向量, nx 為場點 (x) 的單位法向量。 若將式 (3) 與 (4) 中 x 場點 (觀察點) 推到邊界 B 上, 可得到邊界點的對偶邊界積分方程式, 表示如下:
πu(x) = C.P.V.
Z
B
T(s, x)u(s)dB(s) − Z
B
U(s, x)t(s)dB(s), x∈ B, (10) πt(x) = H.P.V.
Z
B
M(s, x)u(s)dB(s) − C.P.V.
Z
B
L(s, x)t(s)dB(s), x ∈ B, (11) 其中 C.P.V. 為柯西主值與 H.P.V. 則為阿達馬主值 (亦可稱為阿達馬有限部分)。
基於邊界元素法的邊界離散觀念, 我們可將問題的邊界分割成 n 個常數元素後, 方程式 (10) 與 (11) 式可離散化後得到如下線性代數方程:
[U] {t} = [T]{u}, (12) [L] {t} = [M]{u}, (13) 其中, [U]、[T]、[L]、[M]為對應 U、T、L、M 四個核函數影響係數矩陣, 將 Dirichlet 型邊界條 件 (2) 式代入 (12) 與 (13) 式可得:
[U] {t} = {0}, (14)
[L] {t} = {0}, (15) 將影響係數矩陣 [U] 與 [L] 作奇異值分解, 可表示如下:
[U] =ΦU ΣU ΨUH
, (16)
[L] =ΦL ΣL ΨLH
. (17)
其中 [Φ] 為左酉矩陣, [Σ] =
σ1
σ2 σ3
. ..
σn
為奇異值對角矩陣 (σ1 < σ2 <
σ3 <· · · < σn), [Ψ] 為右酉矩陣, H 為共軛轉置 (Hermitian) 運算子。 最後, 將影響係數矩 陣 [U] 與 [L] 做奇異值分解後所得的最小奇異值 (σ1) 對波數 (k) 繪圖, 由圖上尋找掉下來 的位置就是特徵值的發生處, 此法稱為直接搜尋法。 得到的特徵值後, 在奇異值很小時 (理論為 零) 所對應的右酉向量, 即為邊界特徵向量。
3. 數值結果
為了否定可聽出鼓的形狀的猜想, 以下舉兩例說明之。
3.1.
算例一
本算例中兩種不同幾何外形皆由 7 個等腰直角三角形所組成 (類似七巧板), 其中三角板 三個角的角度分別為 45◦、 45◦ 與 90◦, 邊界條件為固定 (Dirichlet 型) 邊界條件 (u(x) = 0), 如圖 4 所示。 接著我們以對偶邊界元素法來求得兩種不同幾何外型的特徵值, 所以我們將形狀 1 與形狀 2 各別為使用奇異式與超奇異式所得之最小奇異值對波數作圖, 此時圖 5 與圖 6 所掉下 來的位置即為特徵值發生位置。
表 1. 兩種不同幾何外型的前15個特徵值
編號 有限元素法 Driscoll [26] 對偶邊界元素法(180 個元素, 本文) 相對
(966 個元素) 奇異式 超奇異式 誤差值(%)
λ kD =√
λ 形狀 1 (kB) 形狀 2 形狀 1 (kC) 形狀 2 奇異式 超奇異式 1 2.53794399980 1.59309 1.59 1.59 1.59 1.59 0.19 0.19 2 3.65550971352 1.91194 1.91 1.91 1.91 1.91 0.10 0.10 3 5.17555935622 2.27499 2.27 2.27 2.27 2.27 0.22 0.22 4 6.53755744376 2.55686 2.56 2.56 2.55 2.56 0.12 0.27 5 7.24807786256 2.69223 2.69 2.69 2.69 2.69 0.08 0.08 6 9.20929499840 3.03468 3.03 3.03 3.03 3.03 0.15 0.15 7 10.59698569130 3.2553 3.25 3.25 3.25 3.25 0.16 0.16 8 11.54139539560 3.39726 3.4 3.4 3.4 3.4 0.08 0.08 9 12.33700550140 3.51241 3.51 3.51 3.51 3.51 0.07 0.07 10 13.05365405570 3.61298 3.61 3.61 3.61 3.61 0.08 0.08 11 14.31386246430 3.78337 3.78 3.78 3.78 3.78 0.09 0.09 12 15.87130262000 3.98388 3.98 3.98 3.98 3.98 0.10 0.10 13 16.94175168800 4.11604 4.12 4.12 4.11 4.11 0.10 0.15 14 17.66511843680 4.20299 4.2 4.2 4.2 4.2 0.07 0.07 15 18.98106738770 4.35673 4.35 4.35 4.35 4.35 0.15 0.15
表 1 為前 15 個特徵值所發生的位置, 包括有限元素法與邊界元素法的結果。 誤差公式定義
如下:
εr = |kD− kB|
|kD| × 100% (奇異式), εr = |kD − kC|
|kD| × 100% (超奇異式), 其中, εr 為相對誤差值, kD 為 Driscoll [26] 以有限元素法求得的數值解, 且因形狀1 與形狀 2 的有限元素法之解已收斂到相同的值 (966 個元素), 故僅以一欄表示, kB 為對偶邊界元素法對 形狀 1 的奇異式數值解, kC 為對偶邊界元素法對形狀 1 的超奇異式數值解。
另外, 本文亦發現當使用對偶邊界元素法分析此兩種不同幾何形狀的特徵值時, 發現當每 邊切割 10個常數元素的收斂效果不是很好。 若每邊切割20個常數元素, 則可達到不錯的收斂效 果, 因此在此算例中我們將採取每邊切割 20 個常數元素。 最後, 此算例的結果中可看出對偶邊 界元素法的奇異式比超奇異式收斂效果佳。 為了比較我們所得的結果是否準確, 我們將形狀1與 形狀 2的前八個特徵模態圖與文獻中 [26] 找到的結果做比較, 並整理於表2, 且由模態圖可看出 我們所做出的結果與文獻的結果一致吻合。
圖 4. 兩種不同幾何外型的鼓皮 (左為形狀1, 右為形狀2)
圖 5. 以奇異式分析兩種不同幾何外型 (圖4所示) 鼓皮的特徵頻率 – 最小奇異值對波數作圖
表 2. 兩種不同幾何外形前8個特徵模態
PP PP
PP PP
PP P 特徵模態
方法 有限元素法Driscoll [26] 對偶邊界元素法(本文)
形狀 1 形狀 2 形狀 1 形狀 2
k= 1.59309 k= 1.59
k= 1.91194 k= 1.91
k= 2.27499 k= 2.27
k= 2.55686 k= 2.56
k= 2.69223 k= 2.69
k= 3.03468 k= 3.03
k= 3.25530 k= 3.25
k= 3.39726 k= 3.40
圖 6. 以超奇異式分析兩種不同幾何外型 (圖 4 所示) 鼓皮的特徵頻率 – 最小奇異值對波數作圖
3.2.
算例二
本算例兩種不同幾何外形則由 7 個不同於算例一的三角形所組成的 (類似如七巧板), 其三 角板中的三個角的角度分別為 30◦、 70◦ 與 80◦, 邊界條件一樣為固定 (Dirichlet 型) 邊界條 件 (u(x) = 0), 如圖 7 所示。 接著我們以對偶邊界元素法來求得兩種不同幾何外型的特徵值, 所以我們將形狀 3 與形狀 4, 各別為使用奇異式與超奇異式所得之最小奇異值對波數作圖, 如前 所述圖 8 與圖 9 中所掉下來的位置即為特徵值發生位置。 接著, 我們也像算例一一樣, 將此算例 中的形狀 3 的前 10 個奇異式的特徵值 kB 與超奇異式的特徵值跟有限元素法 kD 的前 10 個特 徵值找出 (因形狀 3 與形狀 4 特徵值在 2672 個元素時已收斂到同一值, 故僅以一欄表示), 並算 出相對誤差且以此作比較, 如下表 3。 而此算例幾形狀較算例一複雜, 我們也發現此算例兩個不 同幾何外形所對應的特徵頻率所發生的位置均相同。 同樣地, 在此算例結果一樣可以觀察到對 偶邊界元素法的奇異式比超奇異式收斂效果佳。 另外, 在文獻 [26] 中 (表2 中的有限元素法模 態圖), 像是模態2、 模態5及模態7並未提供模態圖, 所以我們用 MATLAB 的 ToolBox 來求 解並補足當初 Driscoll 沒提供的模態圖, 我們並將其列於表 4 作一整體比較。
圖 7. 兩種不同幾何外型的鼓皮 (左為形狀3, 右為形狀4)
表 3. 兩種不同幾何外型的前10個特徵值
編號 有限元素法 Driscoll [26] 對偶邊界元素法(180 個元素, 本文) 相對
(2672 個元素) 奇異式 超奇異式 誤差值(%)
λ kD =√
λ 形狀 3 (kB) 形狀 4 形狀 3 (kC) 形狀 4 奇異式 超奇異式 1 5.63126379 2.37303 2.37 2.37 2.35 2.35 0.13 0.97 2 7.18148848 2.67983 2.68 2.68 2.68 2.68 0.01 0.01 3 12.7905748 3.57639 3.57 3.57 3.57 3.56 0.18 0.18 4 13.0935554 3.6185 3.61 3.61 3.59 3.6 0.23 0.79 5 17.0680091 4.13134 4.13 4.13 4.13 4.13 0.03 0.03 6 18.8537757 4.34209 4.34 4.34 4.34 4.34 0.05 0.05 7 19.8509471 4.45544 4.45 4.46 4.45 4.45 0.12 0.12 8 24.1803291 4.91735 4.92 4.92 4.92 4.92 0.05 0.05 9 27.5379471 5.24766 5.25 5.25 5.25 5.25 0.04 0.04 10 30.0098327 5.47812 5.47 5.47 5.47 5.46 0.15 0.15
圖 8. 以奇異式分析兩種不同幾何外型圖 7 所示鼓皮的特徵頻率 – 最小奇異值對波數作圖
圖 9. 以超奇異式分析兩種不同幾何外型圖 7 所示鼓皮的特徵頻率 – 最小奇異值對波數作圖
表 4. 兩種不同幾何外形前 7 個特徵模態
特徵模態 有限元素法 Driscoll [26] 對偶邊界元素法(本文) 有限元素法(MATLAB)
形狀 3 形狀 4 形狀 3 形狀 4 形狀 3 形狀 4
k= 2.37303 k= 2.37 k= 2.374
k= 2.67983 k= 2.68 k= 2.679
k= 3.57639 k= 3.57 k= 3.576
k= 3.61850 k= 3.61 k= 3.619
k= 4.1314 k= 4.13 k= 4.132
k= 4.34209 k= 4.34 k= 4.343
k= 4.45544 k= 4.45 k= 4.455
4. 結論
本文主要利用對偶邊界元素法來分析兩個不同幾何形狀鼓皮的特徵值, 進而來重新檢視可 聽出鼓的形狀的猜想。 透過直接搜尋法檢視最小奇異值對波數作圖來找出特徵值。 由圖中可知,
最小奇異值掉落時所對應的 k 值即為特徵值。 並且在數值實驗中發現奇異式比超奇異式在多邊 形領域更適合求解特徵值問題。 利用對偶邊界元素法驗證兩個不同幾何形狀在固定端位移為零 的邊界條件下, 其特徵譜竟是一樣的, 並與 Driscoll 的數值結果來做比對, 此比對結果的相對 誤差值均在 0.23% 以下 (奇異式的結果)。 另外特徵模態和文獻 [26] 的結果比較也均相當吻 合。 總結可知, 可用邊界元素法來否定可聽出鼓的形狀的猜想。 源於不同外形的鼓皮, 可用邊界 元素法分析出一模一樣的特徵譜。
本研究主要鎖定二維單連通多邊領域的 Dirichlet 問題進行探討。 當此成果推導到二維含 孔洞或三維問題時, 亦或 Neumann 條件結果會如何, 實為一開放問題可以繼續探討。
參考資料
1. 蔡亞倫, 我們可以聽出鼓的形狀嗎? 數學傳播, 30 卷4 期, pp. 34-47, 2005。
2. B. Cipra, You can’t hear the shape of a drum, Science, Vol. 255, No. 5052, pp. 1642- 1643, 1992.
3. S. Sridhar and A. Kudrolli, Experiments on not “hearing the shape” of drums, Physical Review Letters, Vol. 72, No. 14, pp. 2175-2178, 1994.
4. S. J. Chapman, Drums that sound the same, The American Mathematical Monthly, Vol. 102, No. 2, pp. 124-138, 1995.
5. M. Kac, Can one hear the shape of a drum? The American Mathematical Monthly, Vol.
73, No. 4, pp. 1-23, 1966.
6. C. Gordon, D. Webb, and S. Wolpert, One cannot hear the shape of a drum, Bulletin of the American Mathematical Society, Vol. 27, No. 1, pp. 134-138, 1992.
7. C. Gordon, D. Webb, and S. Wolpert, Isospectral plane domains and surfaces via Rie- mannian orbifolds, Inventiones Mathematicae, Vol. 110, No. 1, pp.1-22, 1992.
8. P. N. Swarztrauber, On the numerical solution of the Dirichlet problem for a region of general shape, SIAM Journal on Numerical Analysis, Vol. 9, No. 2, pp. 300-306, 1972.
9. H. Urakawa, Bounded domains which are isospectral but not congruent, Annales Sci- entifiques de l’ ´Ecole Normale Superieure, Vol. 15, pp. 441-456, 1982.
10. S. R. Kuo, J. T. Chen and C. X. Huang, Analytical study and numerical experiments for true and spurious eigensolutions of a circular cavity using the real-part dual BEM, International Journal for Numerical Methods in Engineering, Vol. 48, pp. 1401-1422, 2000.
11. J. T. Chen, J. H. Lin, S. R. Kuo and S. W. Chyuan, Boundary element analysis for the Helmholtz eigenvalue problems with a multiply connected domain, Proceedings of The Royal Society Series A, Vol. 457, pp. 2521-2546, 2001.
12. J. T. Chen, J. W. Lee and Y. C. Cheng, On the spurious eigensolutions for the real-part boundary element method, Engineering Analysis with Boundary Elements, Vol. 33, pp.
342-355, 2009.
13. J. T. Chen, J. W. Lee and S. Y. Leu, Analytical investigation for true and spurious eigensolutions of multiply-connected membranes containing elliptical boundaries using
the dual BIEM, International Journal of Solids and Structures, Vol. 48, pp. 729-744, 2011.
14. I. L. Chen, Y. T. Lee, P. S. Kuo and J. T. Chen, On the true and spurious eigenvalues by using the real or imagainary part of the method of fundamental solution, International Journal Computational Methods, Vol. 10, No. 2, pp. 1-13, 2013.
15. 郭柏伸, 奇異值分解法架構在邊界積分方程與基本解法中的真假特徵值研究, 國立臺灣海洋大學碩 士機電與機械研究所碩士學位論文, 2012。
16. 郭柏伸、 陳義麟、 趙玉星、 陳正宗, 我們可以用對偶邊界元素法分析來聽出鼓的形狀嗎? 中華民國 力學學會第三十五屆全國力學會議, 2011。
17. J. T. Chen, J. W. Lee, I. L. Chen and P. S. Kuo, On the null and nonzero fields for true and spurious eigenvalues of annular and confocal elliptical membranes, Engineering Analysis with Boundary Elements, Vol. 37, No. 1, pp. 42-59, 2013.
18. A. Pleijel, A study of certain Green’s funcrions with applications in the theory of vibrating membranes, Arkiv f¨or Matematik, Vol. 2, No. 6, pp. 553-569, 1952.
19. S. Sridhar, Experimental observation of scarred eigenfunctions of chaotic microwave cavities, Physical Review Letters, Vol. 67, No. 7, pp. 785-788, 1991.
20. P. Buser, J. Conway, P. Doyle and K. D. Semmler, Some planar isospectral domains, Internat Math Res Notices, pp. 1-13, 1994.
21. H. Wu, D. W. L. Sprung, and J. Martorell, Numerical investigation of isospectral cavi- ties built from triangles, Physical Review E, Vol. 51, No. 1, pp. 703-708, 1995.
22. L. Ramdas Ram-Mohan, Finite Element and Boundary Element Applications in Quan- tum Mechanics, Oxford University Press, New York; 2002.
23. 李家瑋, 零場邊界積分方程法求解含圓與橢圓邊界之自由振動與水波 Helmholtz 問題, 國立台灣 海洋大學河海工程學系碩士學位論文, 2010。
24. J. R. Kuttler and V. G. Sigillito, Eigenvalues of the Laplacian in two dimensions, SIAM Review, Vol. 26, No. 2, pp. 163-193, 1984.
25. 陳正宗、 洪宏基, 邊界元素法, 第二版, 新世界出版社, 1992。
26. T. A. Driscoll, Eigenmodes of isospectral drums, SIAM Review, Vol. 39, No. 1, pp.
1-17, 1997.
—本文第一作者陳正宗任教國立台灣海洋大學河海工程學系與機械與機電工程學系, 第二作者 郭柏伸為國立台灣海洋大學機械與機電工程學系碩士研究生, 第三作者高怡絹為國立台灣海洋 大學河海工程學系大學部學生。—
