一、前言

在工程上, 許多考量的物理現象均以微分方程式描述之, 故求解微分方程式在工學院必修課程之工程數學中, 佔有很大比重。 實務上, 由於許多條件的限制, 要解決的問題通常很難用解析的方式求解, 而必須仰賴數值方法。 傳統的數值方法如有限差分法(Finite Difference Method)、有限元素法(Finite Element Method)等, 均必須要有網格才能計算, 而網格的品質大大影響到計算結果, 故人員必須長時間接受產生網格的訓練, 才有辦法在使用傳統數值方法時能得心應手。 在一些控制方程式具有基本解 (Fundamental solution) 的問題中, 可利用邊界元素法 (Boundary Element Method) 降一個維度。 三維問題, 只需在邊界曲面上佈設面元素, 而二維問題, 則只需在邊界曲線上佈設線元素, 可大幅減少計算量。 有關邊界元素法之介紹, 可參考陳與洪 (1992) 。 不過, 要在邊界上佈設元素還是需要一些經驗和熟練度。

相較於需要網格的傳統數值方法, 無網格法為近二、 三十年來新興的一種數值方法。 無網格法又有強式與弱式之分。強式法之數學表示式較為直接, 但易有數值不穩定的問題。 而弱式法為了避免數值不穩定, 必須大量使用數值積分, 故雖為目前無網格法之主流, 入門門檻反而比傳統需要網格的數值方法還高 ,。 各無網格法之詳細發展歷程及特點, 可參閱相關文獻 ,,。

強式無網格法可分為以徑向基底函數(Radial Basis Function, RBF)為基底的和以局部多項式 (Local Polynomial) 為基底的。 以徑向基底函數為基底之強式無網格法主要為全域型(Global type), 也有一些為局部型 (Local type)。 所謂全域或局部, 係指進行配點(Collocation)時, 取計算域裡的所有結點來做, 或是只取配點中心之鄰近結點來做。 而以局部多項式為基底之強式無網格法, 則必為局部型, 且須搭配加權最小二乘法 (Weighted-Least-Squares approach) 或移動最小二乘法 (Moving-Least-Squares approach)。 這類強式無網格法, 較具代表性的有Oñate et al. (1996) 之有限配點法 (Finite Point Method), 以及廣義有限差分法 (General Finite Difference Method)。

強式無網格法在做配點時, 就是一個結點只拿一條方程式來用, 故在邊界點上, 理應控制方程式和邊界條件都該被滿足, 但也只能讓邊界條件被滿足。 通常這樣得到的數值計算結果, 在函數值的部分會準確, 但導數值會不準。 這也是強式法較易有數值不穩定問題的主因。 Jin et al. (2005) 和 Wu and Chang (2011) 在他們提出的方法中, 讓控制方程式和邊界條件在做配點時都有被滿足, 結果避免了強式法之數值不穩定, 並大幅改善數值計算結果之準確性。

Wu and Tsay (2013) 受到 Jin et al. (2005) 和 Wu and Chang (2011) 的啟發, 把 Oñate et al. (1996) 之有限配點法拿來做改良, 而提出一個新的強式無網格法。 他們係利用懲罰權重來將控制方程式及邊界條件嵌入局部近似中, 故該強式法具有強健性。 他們提出這個方法時, 還沒用到「嵌入」這個詞, 而是稱它為「強健的(Robust)局部多項式配點法」, 也稱之為「改良型有限配點法」。 由於 有限配點法還有其他改良版, 而且不只這個方法具有強健性, 故Chiou and Wu (2019) 在他們的論文中, 為了做區隔, 並突顯其特徵, 就參考Jin et al. (2005) 的做法, 用「嵌入」這個詞來稱呼 Wu and Tsay (2013) 之數值方法。值得一提的是, 用懲罰權重來將控制方程式及邊界條件嵌入局部近似中, 為 Wu and Tsay (2013) 首創。

雖然 Wu and Tsay (2013) 之方法的數學表示式已較其他無網格法簡潔許多, 且已被應用於求解若干二維及三維的問題 $-$, 但至今有在使用該方法的人仍佔極少數。 為了將它推廣, 我們特別選一條常微分方程式, 依計算步驟詳細列算式, 希望讓本文的讀者能輕易理解, 並且後續能做進一步的應用。

二、加權最小二乘法局部多項式近似

此數值方法應用於二維及三維邊界值問題之數學式, 可參考相關文獻 $-$。 我們在本文以一維邊界值問題來介紹之。 在計算範圍 $x_L\le x\le x_R$ 中, 函數 $y(x)$ 滿足控制方程式和邊界條件 \begin{equation} \begin{split} &y''+p(x)y'+q(x)y=r(x),\qquad x_L\le x\le x_R,\\ &s(x)y'+t(x)y=u(x),\ \hbox{at}\ x=x_L\ \hbox{and}\ x=x_R. \end{split}\label{1} \end{equation}

我們將計算範圍離散成 $N$ 個結點, 然後目標就是求出 $y(x_j)$ $(j=1,2,\ldots,N)$ 之近似值, 亦即數值解。

在第 $j$ 個結點附近之函數 $y$ 可用局部多項式來近似 \begin{equation} y\Big|_{x\approx x_j}\approx \alpha_{j1}+\alpha_{j2}X+\alpha_{j3}X^2/2+\alpha_{j4}X^3/6 ,\label{2} \end{equation} 其中 $X=x-x_j$, 乃以 $x_j$ 為中心之局部坐標, 而 $\alpha_{j1}$、 $\alpha_{j2}$、 $\alpha_{j3}$、 $\alpha_{j4}$ 這些為待決定的係數, 亦代表函數 $y$ 在 $x=x_j$ 那個點之函數近似值, 以及其一階、二階、三階等導數之近似值。 式 \eqref{2} 可用較為精簡且通用的寫法表示 \begin{equation} y\Big|_{x\approx x_j}\approx \sum_{i=1}^m \alpha_{ji}f_i(X),\label{3} \end{equation} 其中, 就一維問題而言, 以三次多項來做局部近似, 則 $m=4$。 我們需要超過 $m$ 個結點來建立聯立方程組, 並利用加權最小二乘法, 形成之聯立方程組以矩陣式表示如下 \begin{equation} {A\alpha}_j={\beta}\label{4} \end{equation} 其中, \begin{equation} {A}=\left[\begin{array}{ccccc} a_{{1}1} & a_{{1}2} & \cdots & \cdots & a_{{1}m} \\ a_{{2}1}&a_{{2}2}& & &a_{{2}m}\\ \vdots&&a_{{k}i}&&\vdots\\ \vdots&&&\ddots&\vdots\\ a_{{n_{loc}}1}&a_{{n_{loc}}2}&\cdots&\cdots&a_{{n_{loc}}m}\\ \end{array}\right], \label{5} \end{equation} \begin{equation} a_{{k}i}=w_{{k}}f_i(x_{\underline{k}}-x_j),\label{6} \end{equation} \begin{equation} w_{{k}}=\sqrt{W_{j\underline{k}}},\label{7} \end{equation} \begin{equation} {\alpha}_j=\Big[\alpha_{j1} \cdots \alpha_{ji} \cdots \alpha_{jm}\Big]^T,\label{8} \end{equation} \begin{equation} {\beta}=\Big[w_1y_{\underline{1}} \cdots w_ky_{\underline{k}} \cdots w_{{n_{loc}}} y_{\underline{n_{loc}}}\Big]^T,\label{9} \end{equation} 其中, $y_{\underline{k}}=y(x_{\underline{k}})$； $k=1,\ldots,n_{loc}$, 畫底線乃為了強調它是局部編號, 每個結點 $x_l$ ($l=1,\ldots,N$), 在以 $x_j$ ($j=1,\ldots,N$) 為中心之局部範圍, 都有各自的局部編號； $n_{loc}$ 為局部範圍內結點之數量, 必須大於 $m$； $W_{j\underline{k}}$ 為權重因子, 其值最小為 0 最大為 1, 依據 $x_{\underline{k}}$ 到 $x_j$ 之距離 $r_{j\underline{k}}=|x_{\underline{k}}-x_j|$ 而定。 計算權重因子的公式有許多可用, 本文採用正規化高斯函數 \begin{equation} W_{j\underline{k}}=\left\{\begin{array}{lcl} \dfrac{\exp(-\varepsilon(r_{j\underline{k}}/\rho_j)^2)-\exp(-\varepsilon)}{1-\exp(-\varepsilon)},&\quad &r_{j\underline{k}}\lt \rho_j,\\[5pt] 0,&&r_{j\underline{k}}\ge\rho_j, \end{array}\right.\label{10} \end{equation} 其中 $\varepsilon$ 為形狀參數, 而 $\rho_j$ 表示局部範圍之大小。 在此強調, 每個局部範圍內, 各有其局部中心 $x_j$, 各有其局部範圍大小 $\rho_j$, 會包含到 $n_{loc}$ 個結點, 而 $n_{loc}$ 之值可設為定值, 也可隨局部中心 $x_j$ 之位置而異。 若要分析的問題是二維或三維的, 則局部範圍包含的到結點數通常不是定值, 就真的會隨局部中心之位置而異。

若式 \eqref{2} 之 $\alpha_{j1}$ 為函數 $y$ 之真值 $y_j$ 而非其近似值, 則除去式 \eqref{8} 之 $\alpha_{j1}$, 刪去式 \eqref{5} 中之第一行 (column) 及 $l=j$ 所對應的那一列 (row), 並除去式 \eqref{9} 中 $l=j$ 所對應的那一項 (注意 : 第 $l$ 號結點在第 $j$ 個局部範圍之局部編號為 $\underline{k}$) 且將式 \eqref{9} 中的 $y_{\underline{k}}$ 全部代換成 $y_{\underline{k}}-y_j$, 即為廣義有限差分法。

式 \eqref{4} 為沒有正解之線性系統, 其最小二乘解為 \begin{equation} {\alpha}_j=({A}^T{A})^{-1}{A}^T{\beta}.\label{11} \end{equation} 一樣再強調一次, 每個 $x_j$ 附近之局部近似, 也都各自有其對應的 ${A}$。

三、以範例說明傳統強式法

接下來我們選擇以下邊界值問題來當範例。 \begin{equation} y''-10y'=0,\qquad y(0)=1,\qquad y(1)=2.\label{12} \end{equation} 把 $0\le x\le 1$ 的範圍切成 5 等分, 加上頭尾, 總共有 6 個結點, 分別為 $x_1=0$、 $x_2=0.2$、 $x_3=0.4$、 $x_4=0.6$、 $x_5=0.8$、 $x_6=1$。 我們設定局部範圍為能夠涵蓋到至少 5 個結點之最小距離再乘以 1.01, 因此, 在這個範例中, $\rho_1=\rho_6=0.808$, $\rho_2=\rho_5=0.606$, $\rho_3=\rho_4=0.404$。 我們取 $\varepsilon=12$。 為了頁面整潔, 所有算式一律只顯示到小數點以下 6 位並四捨五入。 而讀者需注意, 在我們的運算中, 所有的數值運算均保有 14 位有效位數。

對於第一個結點, 即 $x_1$, 由式 \eqref{4} 可得 \begin{eqnarray} &&\hskip -20pt \left[\begin{array}{cccc} 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 \end{array}\right] \left[\begin{array}{c} \alpha_{11}\\ \alpha_{12}\\ \alpha_{13}\\ \alpha_{14} \end{array}\right]\nonumber\\ &=& \left[\begin{array}{ccccc} 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 \end{array}\right] \left[\begin{array}{c} y_{1}\\ y_{2}\\ y_{3}\\ y_{4}\\ y_{5} \end{array}\right],\label{13} \end{eqnarray} 然後其最小二乘解為 \begin{equation} \left[\begin{array}{c} \alpha_{11}\\ \alpha_{12}\\ \alpha_{13}\\ \alpha_{14} \end{array}\right]\!\!=\!\! \left[\begin{array}{rrrrr} 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 \end{array}\right]\! \left[\begin{array}{c} y_{1}\\ y_{2}\\ y_{3}\\ y_{4}\\ y_{5} \end{array}\right].\label{14} \end{equation} 同理, 在第二個結點, 即 $x_2$, 由式 \eqref{4} 可得 \begin{eqnarray} &&\hskip -20pt \left[\begin{array}{cccc} 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 \end{array}\right] \left[\begin{array}{c} \alpha_{21}\\ \alpha_{22}\\ \alpha_{23}\\ \alpha_{24} \end{array}\right]\nonumber\\ &=& \left[\begin{array}{ccccc} 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 \end{array}\right] \left[\begin{array}{c} y_{1}\\ y_{2}\\ y_{3}\\ y_{4}\\ y_{5} \end{array}\right],\label{15a} \end{eqnarray} 然後其最小二乘解為 \begin{equation} \left[\!\begin{array}{c} \alpha_{21}\\ \alpha_{22}\\ \alpha_{23}\\ \alpha_{24} \end{array}\!\right]\!\!\!=\!\!\! \left[\!\begin{array}{rrrrr} 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 \end{array}\!\right]\!\! \left[\!\begin{array}{c} y_{1}\\[-2pt] y_{2}\\[-2pt] y_{3}\\[-2pt] y_{4}\\[-2pt] y_{5} \end{array}\!\right]\!.\label{15} \end{equation} 依照相同的步驟, 我們可列出第三、 第四、 第五、 第六個結點之局部近似多項式各係數與鄰近點函數值之關聯。 為節省頁面空間, 不再詳列, 讀者可自行嘗試。 在此需注意, 以 $x_4$、 $x_5$ 或 $x_6$ 為中心時, 局部範圍內之結點均為 $x_2,\ldots,x_6$, 不包含 $x_1$。 讀者可依以上所述自行研判各結點在各局部範圍內之編號。 比如, 全域第 2 號結點在第 4 個局部範圍之局部編號為 1, 全域第 6 號結點在第 6 個局部範圍之局部編號為 5。

接下來, 我們將 $j=2,\ldots,5$ 之導數近似值當成真值, 代入控制方程式 \begin{equation} \left\{\begin{array}{l} y'=\alpha_{j2}\\ y''=\alpha_{j3} \end{array}\right.\ \Rightarrow y''-10y'=0\Rightarrow \ \alpha_{j3}-10\alpha_{j2}=0,\label{16} \end{equation} 再利用邊界條件 $y_1=1$ 和 $y_6=2$, 即可組成全域矩陣式如下： \begin{equation} \left[\begin{array}{rrrrrr} 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 \end{array}\right]\!\! \left[\begin{array}{c} y_{1}\\ y_{2}\\ y_{3}\\ y_{4}\\ y_{5}\\ y_{6} \end{array}\right]\!\!=\!\! \left[\begin{array}{c} 1\\ 0\\ 0\\ 0\\ 0\\ 2 \end{array}\right]\!. \label{17} \end{equation} 式 \eqref{16} 中, 各個點的 $\alpha_{j2}$ 和 $\alpha_{j3}$ 都是如同式 \eqref{14}、 式 \eqref{15} 之局部近似關係式那樣, 為 $y_1,\ldots,y_6$ 之線性組合。 式 \eqref{17} 之第一列和最後一列分別為邊界條件, 而第二至五列則為將各內部結點之局部近似關係式代入式 \eqref{16} 而得。 最後求解式 \eqref{17} 可得到 \begin{equation} \left[\begin{array}{c} y_{1}\\[-2pt] y_{2}\\[-2pt] y_{3}\\[-2pt] y_{4}\\[-2pt] y_{5}\\[-2pt] y_{6} \end{array}\right]\!\!=\!\! \left[\begin{array}{c} 1.000000\\[-2pt] 1.009076\\[-2pt] 1.013945\\[-2pt] 1.069128\\[-2pt] 1.178573\\[-2pt] 2.000000 \end{array}\right]. \label{18} \end{equation} 以上這個就是傳統的強式法。 在傳統的強式法中, 並沒有用到讓邊界點滿足控制方程式的步驟。

四、將控制方程式和邊界條件嵌入局部近似

基本上, 控制方程式在任一結點都要被滿足。 依據局部近似之基本定義, 式 \eqref{1} 之控制方程式可表示為 \begin{equation} {\alpha}_{j3}+p_j{\alpha}_{j2}+q_j{\alpha}_{j1}=r_j,\label{19} \end{equation} 其中, $p_j$、 $q_j$、 $r_j$ 表示函數 $p(x)$、 $q(x)$、 $r(x)$ 在 $x=x_j$ 處之值。 若是在邊界點, 則除了控制方程式, 還有邊界條件也該被滿足。 用局部近似表示為 \begin{equation}\begin{split} &s(x)y'(x)+t(x)y(x)=u(x)\quad \hbox{at}\ x=x_j\\ \Rightarrow\ &s_j\alpha_{j2}+t_j\alpha_{j1}=u_j, \label{20} \end{split}\end{equation} 其中, $s_j$、 $t_j$、 $u_j$ 表示函數 $s(x)$、 $t(x)$、 $u(x)$ 在 $x=x_j$ 處之值。

把邊界點將控制方程式與邊界條件都嵌入局部近似中, 則式 \eqref{4} 可俢改為 \begin{equation} \left[\begin{array}{c} {A}\\ {A}' \end{array}\right]{\alpha}_j= \left[\begin{array}{c} {\beta}\\ {\beta}' \end{array}\right], \label{21} \end{equation} 其中, 若 $x=x_j$ 為內部點, 則 \begin{eqnarray} {A}'&=&[ w'q_j\quad w'p_j\quad w'\quad 0 ],\label{22}\\ {\beta}'&=&[ w'r_j ].\label{23} \end{eqnarray} 若 $x=x_j$ 為邊界點, 則 \begin{eqnarray} {A}'&=&\left[\begin{array}{cccc} w'q_j&w'p_j& w'& 0 \\ w't_j&w's_j&0& 0 \end{array}\right],\label{24}\\ {\beta}'&=&\left[\begin{array}{c} w'r_j\\ w'u_j \end{array}\right],\label{25} \end{eqnarray} 其中 $w'$ 稱為懲罰權重, 其值要比 1 大 2 至 3 個階次(order)。 其效用就是要儘可能減少控制方程式及邊界條件之誤差, 而把誤差分配給式 \eqref{21} 聯立方程組中較遠離局部中心之結點所對應的方程式。 式 \eqref{21} 一樣是一個沒有正解的線性系統, 其最小二乘解為 \begin{equation} {\alpha}_j=\left(\left[\begin{array}{c} {A}\\ {A}' \end{array}\right]^T \left[\begin{array}{c} {A}\\ {A}' \end{array}\right] \right)^{-1} \left[\begin{array}{c} {A}\\ {A}' \end{array}\right]^T \left[\begin{array}{c} {\beta}\\ {\beta}' \end{array}\right]. \label{26} \end{equation}

五、以範例說明嵌入式強式法

一樣用式 \eqref{12} 之一維邊界值問題來做範例。 取 $w'=100$, 則由式 \eqref{21} , 我們在第一個結點可得到 \begin{eqnarray} &&\hskip -20pt \left[\begin{array}{rrrr} 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 \end{array}\right] \left[\begin{array}{c} \alpha_{11}\\ \alpha_{12}\\ \alpha_{13}\\ \alpha_{14} \end{array}\right]\label{27}\\ &=& \left[\begin{array}{rrrrrrr} \,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 \end{array}\right] \left[\begin{array}{c} y_{1}\\ y_{2}\\ y_{3}\\ y_{4}\\ y_{5}\\ 0\\ 1 \end{array}\right].\nonumber \end{eqnarray} 讀者應可發現, 式 \eqref{27} 和式 \eqref{13} 只差在式 \eqref{27} 的聯立方程組比式 \eqref{13} 的聯立方程組多了兩條方程式。 這兩條就是我們利用懲罰權重強迫讓控制方程式和邊界條件在這個結點被滿足而加入的。 然後由最小二乘法, 我們可得到其局部近似多項式各係數與鄰近點函數值之間關聯如下： \begin{eqnarray} \left[\begin{array}{c} \alpha_{11}\\ \alpha_{12}\\ \alpha_{13}\\ \alpha_{14} \end{array}\right] &=& \left[\begin{array}{rrrrr} 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 \end{array}\right] \left[\begin{array}{c} y_{1}\\ y_{2}\\ y_{3}\\ y_{4}\\ y_{5} \end{array}\right]\nonumber\\ &&+\left[\begin{array}{r} 0.999899\\ -3.222565\\ -32.225650\\ 247.067993 \end{array}\right].\label{28} \end{eqnarray} 在第二結點, 我們可由式 \eqref{21} 得到 \begin{eqnarray} &&\hskip -20pt \left[\begin{array}{rrrr} 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 \end{array}\right] \left[\begin{array}{c} \alpha_{21}\\ \alpha_{22}\\ \alpha_{23}\\ \alpha_{24} \end{array}\right]\label{29}\\ &=& \left[\begin{array}{rrrrrr} 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 \end{array}\right] \left[\begin{array}{c} y_{1}\\ y_{2}\\ y_{3}\\ y_{4}\\ y_{5}\\ 0 \end{array}\right].\nonumber \end{eqnarray} 然後由最小二乘法, 我們可得到其局部近似多項式各係數與鄰近點函數值之間關聯如下： \begin{equation} \left[\begin{array}{c} \alpha_{21}\\ \alpha_{22}\\ \alpha_{23}\\ \alpha_{24} \end{array}\right] \!\!=\!\! \left[\begin{array}{rrrrr} 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 \end{array}\right] \left[\begin{array}{c} y_{1}\\[-3pt] y_{2}\\[-3pt] y_{3}\\[-3pt] y_{4}\\[-3pt] y_{5} \end{array}\right]\!\!+\!\!\left[\begin{array}{c} 0\\ 0\\ 0\\ 0 \end{array}\right].\label{30} \end{equation} 依照相同的步驟, 我們可列出第三、第四、第五、第六個結點之局部近似多項式各係數與鄰近點函數值之關聯。 為節省空間, 不再詳列, 讀者可自行嘗試。 為了組成全域矩陣式而求得$y_1$到$y_6$等所有結點的函數值, 我們取式 \eqref{28} 的第一列, 將近似值當成真值, 而列出全域矩陣式的第一列 \begin{equation} \begin{split} &y_1\!=\![0.000100\ 0.000002\ 0.000001\ 0.000000\ 0.000000]\left[\begin{array}{c} y_{1}\\[-3pt] y_{2}\\[-3pt] y_{3}\\[-3pt] y_{4}\\[-3pt] y_{5} \end{array}\right]\!+\!0.999899\\ \Rightarrow&[-0.999900\ 0.000002\ -0.000001\ 0.000000\ 0.000000] \left[\begin{array}{c} y_{1}\\[-3pt] y_{2}\\[-3pt] y_{3}\\[-3pt] y_{4}\\[-3pt] y_{5} \end{array}\right]\!=\!-0.999899. \end{split}\label{31} \end{equation} 同理, 我們可弄出全域矩陣式的第二列至第六列, 最後組成全域矩陣式如下： \begin{eqnarray} &&\hskip -20pt \left[\begin{array}{rrrrrr} -0.999900&0.000002&-0.000001&0.000000&0.000000&0.000000\\[-2pt] 0.046926&-0.028219&-0.028034&0.009308&0.000018&0.000000\\[-2pt] -0.000015&0.000157&-0.000147&0.000000&0.000005&0.000000\\[-2pt] 0.000000&-0.000015&0.000157&-0.000147&0.000000&0.000005\\[-2pt] 0.000000&-0.000019&-0.015727&0.142112&-0.142150&0.015784\\[-2pt] 0.000000&0.000000&-0.000001&0.000002&0.000002&-0.999900 \end{array}\right] \left[\begin{array}{c} y_{1}\\[-2pt] y_{2}\\[-2pt] y_{3}\\[-2pt] y_{4}\\[-2pt] y_{5}\\[-2pt] y_{6} \end{array}\right]\nonumber\\ &=&\left[\begin{array}{r} -0.999899\\[-2pt] 0.000000\\[-2pt] 0.000000\\[-2pt] 0.000000\\[-2pt] 0.000000\\[-2pt] -1.999795 \end{array}\right],\label{32} \end{eqnarray} 並求解得到 \begin{equation} \left[\begin{array}{c} y_{1}\\[-2pt] y_{2}\\[-2pt] y_{3}\\[-2pt] y_{4}\\[-2pt] y_{5}\\[-2pt] y_{6} \end{array}\right]=\left[\begin{array}{r} 1.000000\\[-2pt] 1.004545\\[-2pt] 1.009960\\[-2pt] 1.043472\\[-2pt] 1.153395\\[-2pt] 1.999998 \end{array}\right].\label{33} \end{equation}

六、比較與討論

本文所選之邊界值問題, 其正解為 $2-(1-e^{10(x+1)})/(1-e^{-10})$。 若讀者熟悉有限差分法, 亦可利用二階中央差分列出算式並求解得到 $y_1=y_2=y_3=y_4=y_5=1$, 而 $y_6=2$。 相較之下, 不論是傳統的強式無網格法, 或是 Wu and Tsay (2013) 所提出將控制方程式及邊界條件嵌入局部近似多項式中的強式無網格法, 都比有限差分法具有優勢。 另計算傳統強式法與 Wu and Tsay (2013) 方法之誤差方均根, 分別得到 0.027896 和 0.013149, 顯示 Wu and Tsay (2013) 方法又較為準確。

讀者了解本方法之計算步驟之後, 可自行設定式 \eqref{10} 裡的 $\varepsilon$ 為不同值, 比如設為 10 或設為 15, 或者採用更大或較小的懲罰權重, 比如用 200 或 50, 看看結果如何。 也可用更多的結點將計算域離散, 或是採用不等距的佈點方式, 看看成效如何。 也可自行選擇另一個邊界值問題, 嘗試用本文所述之方法求近數值解看看。

參考文獻

---本文通訊作者吳南靖任教國立嘉義大學土木與水資源工程學系, 第一作者林彥廷投稿時為該學系大學生---