以一個簡單的一維邊界值問題 介紹嵌入式無網格法
林彥廷 · 吳南靖
一、 前言
在工程上, 許多考量的物理現象均以微分方程式描述之, 故求解微分方程式在工學院必 修課程之工程數學中, 佔有很大比重。 實務上, 由於許多條件的限制, 要解決的問題通常很難 用解析的方式求解, 而必須仰賴數值方法。 傳統的數值方法如有限差分法 (Finite Difference Method)、 有限元素法 (Finite Element Method) 等, 均必須要有網格才能計算, 而網格的 品質大大影響到計算結果, 故人員必須長時間接受產生網格的訓練, 才有辦法在使用傳統數值 方法時能得心應手。 在一些控制方程式具有基本解 (Fundamental solution) 的問題中, 可利 用邊界元素法 (Boundary Element Method) 降一個維度。 三維問題, 只需在邊界曲面上佈 設面元素, 而二維問題, 則只需在邊界曲線上佈設線元素, 可大幅減少計算量。 有關邊界元素法 之介紹, 可參考陳與洪 (1992) [1]。 不過, 要在邊界上佈設元素還是需要一些經驗和熟練度。
相較於需要網格的傳統數值方法, 無網格法為近二、 三十年來新興的一種數值方法。 無網 格法又有強式與弱式之分。 強式法之數學表示式較為直接, 但易有數值不穩定的問題。 而弱式 法為了避免數值不穩定, 必須大量使用數值積分, 故雖為目前無網格法之主流, 入門門檻反而 比傳統需要網格的數值方法還高 [2, 3]。 各無網格法之詳細發展歷程及特點, 可參閱相關文獻 [2, 3, 4]。
強式無網格法可分為以徑向基底函數(Radial Basis Function, RBF) 為基底的和以局 部多項式 (Local Polynomial) 為基底的。 以徑向基底函數為基底之強式無網格法主要為全域 型 (Global type), 也有一些為局部型 (Local type)。 所謂全域或局部, 係指進行配點 (Collo- cation) 時, 取計算域裡的所有結點來做, 或是只取配點中心之鄰近結點來做。 而以局部多項式 為基底之強式無網格法, 則必為局部型, 且須搭配加權最小二乘法 (Weighted-Least-Squares approach) 或移動最小二乘法 (Moving-Least-Squares approach)。 這類強式無網格法, 較 具代表性的有 O˜nate et al. (1996) [5] 之有限配點法 (Finite Point Method), 以及廣義有 限差分法 (General Finite Difference Method)。
62
強式無網格法在做配點時, 就是一個結點只拿一條方程式來用, 故在邊界點上, 理應控制 方程式和邊界條件都該被滿足, 但也只能讓邊界條件被滿足。 通常這樣得到的數值計算結果, 在 函數值的部分會準確, 但導數值會不準。 這也是強式法較易有數值不穩定問題的主因。 Jin et al. (2005) [6] 和 Wu and Chang (2011) [7] 在他們提出的方法中, 讓控制方程式和邊界條件 在做配點時都有被滿足, 結果避免了強式法之數值不穩定, 並大幅改善數值計算結果之準確性。
Wu and Tsay (2013) [8] 受到 Jin et al. (2005) [6] 和 Wu and Chang (2011) [7]
的啟發, 把 O˜nate et al. (1996) [5] 之有限配點法拿來做改良, 而提出一個新的強式無網格 法。 他們係利用懲罰權重來將控制方程式及邊界條件嵌入局部近似中, 故該強式法具有強健性。
他們提出這個方法時, 還沒用到 「嵌入」 這個詞, 而是稱它為 「強健的(Robust) 局部多項式配 點法」, 也稱之為 「改良型有限配點法」。 由於有限配點法還有其他改良版, 而且不只這個方法具 有強健性, 故 Chiou and Wu (2019) [9] 在他們的論文中, 為了做區隔, 並突顯其特徵, 就參 考Jin et al. (2005) [6] 的做法, 用 「嵌入」 這個詞來稱呼 Wu and Tsay (2013) [8] 之數值方 法。 值得一提的是, 用懲罰權重來將控制方程式及邊界條件嵌入局部近似中, 為 Wu and Tsay (2013) [8]首創。
雖然 Wu and Tsay (2013) [8] 之方法的數學表示式已較其他無網格法簡潔許多, 且已 被應用於求解若干二維及三維的問題 [8]−[12], 但至今有在使用該方法的人仍佔極少數。 為了 將它推廣, 我們特別選一條常微分方程式, 依計算步驟詳細列算式, 希望讓本文的讀者能輕易理 解, 並且後續能做進一步的應用。
二、 加權最小二乘法局部多項式近似
此數值方法應用於二維及三維邊界值問題之數學式, 可參考相關文獻 [8]−[12]。 我們在本 文以一維邊界值問題來介紹之。 在計算範圍 xL ≤ x ≤ xR 中, 函數 y(x) 滿足控制方程式和 邊界條件
y′′+ p(x)y′+ q(x)y = r(x), xL≤ x ≤ xR,
s(x)y′+ t(x)y = u(x), at x = xL and x = xR. (1) 我們將計算範圍離散成 N 個結點, 然後目標就是求出 y(xj) (j = 1, 2, . . . , N) 之近似 值, 亦即數值解。
在第 j 個結點附近之函數 y 可用局部多項式來近似 y
x≈xj
≈ αj1+ αj2X + αj3X2/2 + αj4X3/6, (2) 其中 X = x − xj, 乃以 xj 為中心之局部坐標, 而 αj1、 αj2、 αj3、 αj4 這些為待決定的係數, 亦代表函數 y 在 x = xj 那個點之函數近似值, 以及其一階、 二階、 三階等導數之近似值。 式
(2) 可用較為精簡且通用的寫法表示 y
x≈xj
≈
m
X
i=1
αjifi(X), (3)
其中, 就一維問題而言, 以三次多項來做局部近似, 則 m = 4。 我們需要超過 m 個結點來建立 聯立方程組, 並利用加權最小二乘法, 形成之聯立方程組以矩陣式表示如下
Aαj = β (4)
其中,
A=
a11 a12 · · · a1m
a21 a22 a2m
... aki ... ... . .. ... anloc1anloc2 · · · anlocm
, (5)
aki= wkfi(xk− xj), (6)
wk =pWjk, (7)
αj =h
αj1 · · · αji · · · αjm
iT
, (8)
β=h
w1y1 · · · wkyk · · · wnlocynloc
iT
, (9)
其中, yk = y(xk); k = 1, . . . , nloc, 畫底線乃為了強調它是局部編號, 每個結點 xl (l = 1, . . . , N), 在以 xj (j = 1, . . . , N) 為中心之局部範圍, 都有各自的局部編號; nloc 為局部範 圍內結點之數量, 必須大於 m; Wjk 為權重因子, 其值最小為 0 最大為 1, 依據 xk 到 xj 之距 離 rjk= |xk− xj| 而定。 計算權重因子的公式有許多可用, 本文採用正規化高斯函數
Wjk =
exp(−ε(rjk/ρj)2) − exp(−ε)
1 − exp(−ε) , rjk< ρj,
0, rjk≥ ρj,
(10)
其中 ε 為形狀參數, 而 ρj 表示局部範圍之大小。 在此強調, 每個局部範圍內, 各有其局部中心 xj, 各有其局部範圍大小 ρj, 會包含到 nloc 個結點, 而 nloc 之值可設為定值, 也可隨局部中 心 xj 之位置而異。 若要分析的問題是二維或三維的, 則局部範圍包含的到結點數通常不是定 值, 就真的會隨局部中心之位置而異。
若式 (2) 之 αj1 為函數 y 之真值 yj 而非其近似值, 則除去式 (8) 之 αj1, 刪去式 (5) 中 之第一行 (column) 及 l = j 所對應的那一列 (row), 並除去式 (9) 中 l = j 所對應的那一 項 (注意 : 第 l 號結點在第 j 個局部範圍之局部編號為 k) 且將式 (9) 中的 yk 全部代換成 yk− yj, 即為廣義有限差分法。
式 (4) 為沒有正解之線性系統, 其最小二乘解為
αj = (ATA)−1ATβ. (11) 一樣再強調一次, 每個 xj 附近之局部近似, 也都各自有其對應的 A。
三、 以範例說明傳統強式法
接下來我們選擇以下邊界值問題來當範例。
y′′− 10y′ = 0, y(0) = 1, y(1) = 2. (12) 把 0 ≤ x ≤ 1 的範圍切成 5 等分, 加上頭尾, 總共有 6 個結點, 分別為 x1 = 0、 x2 = 0.2、
x3 = 0.4、 x4 = 0.6、 x5 = 0.8、 x6 = 1。 我們設定局部範圍為能夠涵蓋到至少 5 個結 點之最小距離再乘以 1.01, 因此, 在這個範例中, ρ1 = ρ6 = 0.808, ρ2 = ρ5 = 0.606, ρ3 = ρ4 = 0.404。 我們取 ε = 12。 為了頁面整潔, 所有算式一律只顯示到小數點以下 6 位並 四捨五入。 而讀者需注意, 在我們的運算中, 所有的數值運算均保有 14 位有效位數。
對於第一個結點, 即 x1, 由式 (4) 可得
1.000000 0.000000 0.000000 0.000000 0.692384 0.138477 0.013848 0.000923 0.229811 0.091924 0.018385 0.002451 0.036487 0.021892 0.006568 0.001314 0.001280 0.001024 0.000410 0.000109
α11
α12
α13
α14
=
1.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.692384 0.000000 0.000000 0.000000 0.000000 0.000000 0.229811 0.000000 0.000000 0.000000 0.000000 0.000000 0.036487 0.000000 0.000000 0.000000 0.000000 0.000000 0.001280
y1
y2
y3
y4
y5
, (13)
然後其最小二乘解為
α11
α12
α13
α14
=
0.999998 0.000006 −0.000010 0.000006 −0.000002
−9.157044 14.961509 −7.442263 1.628175 0.009623 49.859417 −124.437666 99.156499 −24.437666 −0.140584
−124.323431 372.293724 −370.940585 122.293724 0.676569
y1
y2
y3
y4
y5
.
(14)
同理, 在第二個結點, 即 x2, 由式 (4) 可得
0.520202 −0.104040 0.010404 −0.000694 1.000000 0.000000 0.000000 0.000000 0.520202 0.104040 0.010404 0.000694 0.073190 0.029276 0.005855 0.000781 0.001280 0.000768 0.000230 0.000046
α21
α22
α23
α24
=
0.520202 0.000000 0.000000 0.000000 0.000000 0.000000 1.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.520202 0.000000 0.000000 0.000000 0.000000 0.000000 0.073190 0.000000 0.000000 0.000000 0.000000 0.000000 0.001280
y1
y2
y3
y4
y5
, (15)
然後其最小二乘解為
α21
α22 α23
α24
=
0.000007 0.999974 0.000039 −0.000026 0.000007
−1.667868 −2.495193 4.992790 −0.828527 −0.001202 24.998619 −49.994478 24.991717 0.005522 −0.001381
−124.831051 374.324206 −373.986309 124.324206 0.168949
y1
y2
y3
y4
y5
. (16)
依照相同的步驟, 我們可列出第三、 第四、 第五、 第六個結點之局部近似多項式各係數與鄰近點 函數值之關聯。 為節省頁面空間, 不再詳列, 讀者可自行嘗試。 在此需注意, 以 x4、 x5 或 x6 為 中心時, 局部範圍內之結點均為 x2, . . . , x6, 不包含 x1。 讀者可依以上所述自行研判各結點在 各局部範圍內之編號。 比如, 全域第 2 號結點在第 4 個局部範圍之局部編號為 1, 全域第 6 號 結點在第 6 個局部範圍之局部編號為 5。
接下來, 我們將 j = 2, . . . , 5 之導數近似值當成真值, 代入控制方程式 (y′ = αj2
y′′ = αj3
⇒ y′′− 10y′= 0 ⇒ αj3− 10αj2= 0, (17) 再利用邊界條件 y1 = 1 和 y6 = 2, 即可組成全域矩陣式如下:
1.000000 0.000000 0.000000 0.000000 0.000000 0.000000 41.677303 −25.042543 −24.936185 8.290790 0.010636 0.000000
−4.163319 58.319944 −49.979916 −8.346723 4.170014 0.000000 0.000000 −4.163319 58.319944 −49.979916 −8.346723 4.170014 0.000000 −0.013397 −8.279746 74.919618 −74.946412 8.319936 0.000000 0.000000 0.000000 0.000000 0.000000 1.000000
y1 y2
y3
y4
y5
y6
=
1 0 0 0 0 2
.
(18)
式 (17) 中, 各個點的 αj2 和 αj3 都是如同式 (14)、 式 (16) 之局部近似關係式那樣, 為 y1, . . . , y6 之線性組合。 式 (18) 之第一列和最後一列分別為邊界條件, 而第二至五列則為將各 內部結點之局部近似關係式代入式 (17) 而得。 最後求解式 (18) 可得到
y1
y2
y3
y4 y5 y6
=
1.000000 1.009076 1.013945 1.069128 1.178573 2.000000
. (19)
以上這個就是傳統的強式法。 在傳統的強式法中, 並沒有用到讓邊界點滿足控制方程式的步驟。
四、 將控制方程式和邊界條件嵌入局部近似
基本上, 控制方程式在任一結點都要被滿足。 依據局部近似之基本定義, 式 (1) 之控制方 程式可表示為
αj3+ pjαj2+ qjαj1 = rj, (20) 其中, pj、 qj、 rj 表示函數 p(x)、 q(x)、 r(x) 在 x = xj 處之值。 若是在邊界點, 則除了控制方 程式, 還有邊界條件也該被滿足。 用局部近似表示為
s(x)y′(x) + t(x)y(x) = u(x) at x = xj
⇒ sjαj2+ tjαj1 = uj, (21) 其中, sj、 tj、 uj 表示函數 s(x)、 t(x)、 u(x) 在 x = xj 處之值。
把邊界點將控制方程式與邊界條件都嵌入局部近似中, 則式 (4) 可俢改為
"
A A′
# αj =
"
β β′
#
, (22)
其中, 若 x = xj 為內部點, 則
A′= [ w′qj w′pj w′ 0 ], (23)
β′= [ w′rj ]. (24)
若 x = xj 為邊界點, 則
A′=
"
w′qjw′pj w′ 0 w′tj w′sj 0 0
#
, (25)
β′=
"
w′rj
w′uj
#
, (26)
其中 w′ 稱為懲罰權重, 其值要比 1 大 2 至 3 個階次(order)。 其效用就是要儘可能減少控制 方程式及邊界條件之誤差, 而把誤差分配給式 (22) 聯立方程組中較遠離局部中心之結點所對 應的方程式。 式 (22) 一樣是一個沒有正解的線性系統, 其最小二乘解為
αj =
"
A A′
#T "
A A′
#
−1"
A A′
#T "
β β′
#
. (27)
五、 以範例說明嵌入式強式法
一樣用式 (12) 之一維邊界值問題來做範例。 取 w′ = 100, 則由式 (22) , 我們在第一個 結點可得到
1.000000 0.000000 0.000000 0.000000 0.692384 0.138477 0.013848 0.000923 0.229811 0.091924 0.018385 0.002451 0.036487 0.021892 0.006568 0.001314 0.001280 0.001024 0.000410 0.000109 0.000000 − 1000.000000 100.000000 0.000000 100.000000 0.000000 0.000000 0.000000
α11
α12
α13 α14
(28)
=
1.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.692384 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.229811 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.036487 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.001280 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 100.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 100.000000
y1 y2 y3
y4
y5
0 1
.
讀者應可發現, 式 (28) 和式 (13) 只差在式 (28) 的聯立方程組比式 (13) 的聯立方程組多了兩 條方程式。 這兩條就是我們利用懲罰權重強迫讓控制方程式和邊界條件在這個結點被滿足而加
入的。 然後由最小二乘法, 我們可得到其局部近似多項式各係數與鄰近點函數值之間關聯如下:
α11
α12
α13
α14
=
0.000100 0.000002 −0.000001 0.000000 0.000000
−0.000322 3.462731 −0.159202 −0.080323 −0.000319
−0.003223 34.627311 −1.592019 −0.803233 −0.003185 0.024707 −347.410951 84.906268 15.356922 0.055062
y1 y2 y3
y4
y5
+
0.999899
−3.222565
−32.225650 247.067993
. (29)
在第二結點, 我們可由式 (22) 得到
0.520202 −0.104040 0.010404 − 0.000694 1.000000 0.000000 0.000000 0.000000 0.520202 0.104040 0.010404 0.000694 0.073190 0.029276 0.005855 0.000781 0.001280 0.000768 0.000230 0.000046 0.000000 −1000.000000 100.000000 0.000000
α21
α22
α23 α24
(30)
=
0.520202 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.520202 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.073190 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.001280 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 100.000000
y1
y2
y3
y4
y5
0
.
然後由最小二乘法, 我們可得到其局部近似多項式各係數與鄰近點函數值之間關聯如下:
α21
α22
α23
α24
=
0.046926 0.971781 −0.028034 0.009308 0.000018 1.975491 −4.684370 2.812911 −0.103760 −0.000272 19.754909 −46.843700 28.129109 −1.037599 −0.002719
−498.382621 598.779261 −150.484541 50.014276 0.073620
y1
y2
y3 y4
y5
+
0 0 0 0
.
(31) 依照相同的步驟, 我們可列出第三、 第四、 第五、 第六個結點之局部近似多項式各係數與鄰近點 函數值之關聯。 為節省空間, 不再詳列, 讀者可自行嘗試。 為了組成全域矩陣式而求得y1到y6等
所有結點的函數值, 我們取式 (29) 的第一列, 將近似值當成真值, 而列出全域矩陣式的第一列
y1= [0.000100 0.000002 0.000001 0.000000 0.000000]
y1
y2
y3 y4
y5
+0.999899
⇒[−0.999900 0.000002 − 0.000001 0.000000 0.000000]
y1
y2
y3
y4
y5
= −0.999899.
(32)
同理, 我們可弄出全域矩陣式的第二列至第六列, 最後組成全域矩陣式如下:
−0.999900 0.000002 −0.000001 0.000000 0.000000 0.000000 0.046926 −0.028219 −0.028034 0.009308 0.000018 0.000000
−0.000015 0.000157 −0.000147 0.000000 0.000005 0.000000 0.000000 −0.000015 0.000157 −0.000147 0.000000 0.000005 0.000000 −0.000019 −0.015727 0.142112 −0.142150 0.015784 0.000000 0.000000 −0.000001 0.000002 0.000002 −0.999900
y1
y2
y3
y4
y5 y6
=
−0.999899 0.000000 0.000000 0.000000 0.000000
−1.999795
, (33)
並求解得到
y1
y2
y3
y4
y5
y6
=
1.000000 1.004545 1.009960 1.043472 1.153395 1.999998
. (34)
六、 比較與討論
本文所選之邊界值問題, 其正解為 2 − (1 − e10(x+1))/(1 − e−10)。 若讀者熟悉有限差分 法, 亦可利用二階中央差分列出算式並求解得到 y1 = y2 = y3 = y4 = y5 = 1, 而 y6 = 2。
相較之下, 不論是傳統的強式無網格法, 或是 Wu and Tsay (2013) [8] 所提出將控制方程式 及邊界條件嵌入局部近似多項式中的強式無網格法, 都比有限差分法具有優勢。 另計算傳統強 式法與 Wu and Tsay (2013) [8] 方法之誤差方均根, 分別得到 0.027896 和 0.013149, 顯示 Wu and Tsay (2013) [8] 方法又較為準確。
讀者了解本方法之計算步驟之後, 可自行設定式 (10) 裡的 ε 為不同值, 比如設為 10 或 設為 15, 或者採用更大或較小的懲罰權重, 比如用 200 或 50, 看看結果如何。 也可用更多的 結點將計算域離散, 或是採用不等距的佈點方式, 看看成效如何。 也可自行選擇另一個邊界值問 題, 嘗試用本文所述之方法求近數值解看看。
參考文獻
1. 陳正宗, 洪宏基。 邊界元素法。 臺北市 : 新世界出版社, 1992。
2. Liu, G. R. and Gu, Y. T., An Introduction to Meshfree Methods and Their Programming, Springer: The Netherlands, 2005.
3. Nguyen, V. P., Rabczuk, T., Bordas, S., and Duflot, M., Meshless methods: a review and computer implementation aspects, Mathematics and Computers in Simulation, 79, 763-813, 2008.
4. Chen, J.-S., Hillman, M. and Chi, S.-W., Meshfree Methods: Progress Made after 20 Years, Journal of Engineering Mechanics, 143, 04017001, 2017.
5. O˜nate, E., Idelsohn, S., Zienkiewicz, O. C., and Taylor, R. L. A finite point method in computational mechanics. Applications to convective transport and fluid flow, Inter- national Journal for Numerical Methods in Engineering, 39, 3839-3866, 1996.
6. Jin, X. Z., Li, G., and Aluru, N. R., New approximations and collocation schemes in the finite cloud method, Computers and Structures, 83, 1366-1385, 2005.
7. Wu, N.-J. and Chang, K.-A., Simulation of free-surface waves in liquid sloshing using a domain type-meshless method, International Journal for Numerical Methods in Fluids, 67, 269-288, 2011.
8. Wu, N.-J. and Tsay, T.-K., A robust local polynomial collocation method, International Journal for Numerical Methods in Engineering, 93, 355-375, 2013.
9. Chiou, Y.-C. and Wu, N.-J., A GE/BC imbedded local polynomial collocation method for two dimensional multivariable problems, Engineering Analysis with Boundary Ele- ments, 100, 185-194, 2019.
10. Wu, N.-J., Tsay, T.-K., and Chen, Y.-Y., Generation of stable solitary waves by a piston-type wave maker, Wave Motion, 51, 240-255, 2014.
11. Wu, N.-J., Hsiao, S.-C., and Wu H.-L., Mesh-free simulation of liquid sloshing subjected to harmonic excitations, Engineering Analysis with Boundary Elements, 64, 90-100, 2016.
12. Hsiao, S.-C., Shih, M.-Y. and Wu, N.-J., Simulation of Propagation and Run-Up of Three Dimensional Landslide-Induced Waves Using a Meshless Method, Water, 10(5), 552, 2018.
—本文通訊作者吳南靖任教國立嘉義大學土木與水資源工程學系, 第一作者林彥廷投稿時為該 學系大學生—
