第十章 配重殘值法
配重殘值法(Method of Weighted Residuals, MWR)的基本觀念，是利用一組多項式表示近似於問題的正確解，將正確解與近似解的誤差殘值，乘以適當的配重函數後，使其總和或積分值最小化，以達到最接近正確解的目標。工程界採用MWR配重殘值法解決工程分析問題，約自1960年代開始，目前已有許多專書 [ 1，2，3 ] 詳細討論這種方法，其應用廣見於流體力學、熱傳遞、質量傳遞、反應工程等各領域，是一種相當有效率的數學工具。配重殘值法中，又以數學配置法最易於使用，近年來已成為解微分方程式的數值方法中最廣泛被使用的方法之一。
第一節 配重殘值法基本原理
為了說明配重殘值法的基本想法與處理方法，我們首先以一個最簡單的邊界值常微分方程式為例加以說明。
(10-1.1)
假設
及
在0＜x＜1的區間內都是連續函數，則以上方程式的解也必然是一個連續函數。根據本書第三章的說明，任何連續函數都可以利用一個多項式來代表﹔因此，方程式(10-1.1)的解，可以利用一無限級數的和表示之。
其中ui(x)稱為基礎函數(Basis Function)。由於x的範圍為0＜x＜1，當i值增大時，xi即快速變小﹔因此，為了簡化起見，考慮當i > N後xi即可忽略， 
可以假設為利用一個(N＋1)次多項式近似之。
(10-1.2)
由於
為方程式(10-1.1)的近似解，因此，
也要滿足微分方程式的兩個邊界條件，BC1及BC2；將
代入邊界條件BC1中，得到：
BC1 
(10-1.3)
將
代入邊界條件BC2中，得到：
BC2 
或將此方程式整理後得到
(10-1.4)
將方程式(10-1.3)及(10-1.4)代回方程式(10-1.2)，則θ的近似解表示式可以改寫成為
(10-1.5)
此方程式中含有C2，C3，……，CN＋1總共N個未定係數，因此，需要再利用N個滿足微分方程式(10-1.1)的條件來決定這些係數。
為了簡化說明起見，首先考慮N＝1的情況，亦即
，則方程式(10-1.5)可簡化成
(10-1.6)
以上步驟為建立微分方程式近似解表示式的方法。其次要考慮如何決定未定係數C2，C3，……，CN＋1。為了簡化說明起見，將微分方程式寫成L{x，y}＝0的形式﹔若近似多項式(10-1.2)完全滿足原微分方程式，則代入微分方程式將得到L{x，
}＝0。但由於如前面所說明
為近似解，將
代入微分方程式L{x，y} 中，結果不一定能得到零。因此，定義
代入微分方程式所得到的殘值為R＝L{x，
}。
將方程式(10-1.6)代入原方程式，得到殘值為
(10-1.7)
若
為正確解，則殘值R＝0。但由方程式(10-1.7)可以知道，由於x位於0與1間(0＜x＜1)，因此，我們若選定某一特定C2值，必能使殘值R在0＜x＜1間的某一點時為零，而在其他x位置的殘值R則可能大於或小於零。可是需注意由於
只是一個近似解，如果殘值R在每一x位置都等於零，當然是最好的解答，否則也可以退而求其次，要求殘值R在所考慮區間0＜x＜1內的某一種平均值為零，亦即，讓殘值R的配重積分為零，亦不失為一良好的近似解。這種觀念以數學式表示，即為
(10-1.8)
亦即，只要給定配重函數
，利用方程式(10-1.8)執行積分後，即可求解未定係數Ci。然後，代回
的表示式(10-1.2)，即可求得微分方程式的近似解
。以上所說明的就是配重殘值法的基本想法與步驟，加以整理如下：
(1) 利用一個多項式
作為微分方程式L{x，y}＝0的近似解；
(2) 使多項式滿足微分方程式的邊界條件，並簡化
表示式﹔
(3) 代入微分方程式L{x，y}＝0，求出殘值的表示式R＝L{x，
}﹔
(4) 定義配重函數
；
(5) 使配重殘值的積分式為零，即
﹔可以建立N個未定係數Ci的聯立方程式；
(6) 解聯立方程式，找出近似多項式的未定係數Ci﹔
(7) 代回原近似多項式
，建立微分方程式的近似解﹔
(8) 利用適當內插法，計算微分方程式在特定位置的解。
配重殘值法(Methods of Weighted Residuals)的基本想法如前面的介紹，但配重函數
有多種不同的選擇。依所選擇配重函數
的不同，對應的配重殘值法即有不同的名稱及不同的特性，以下分別簡要介紹之。
[A] 配置法(Collocation Method)
配置法的配重函數
定義為
(10-1.9)
其中
稱為配置點(Collocation points)。
稱為Dirac Delta函數，假設殘值R＝L{x，
}函數在所考慮的控制體積V之外不存在，則配重殘值積分式可以寫成
這種方法表示在配置點位置上的殘值一定為零。當N值增大時，殘值R(x)會在越多的點上為零，使得近似解趨近於實際解。
蘭若斯(Lanczos, C. 1938)選擇柴比雪夫多項式(Chebyshev Polynomial)作為基礎函數ui(x)，柴比雪夫多項式是一種正交多項式。蘭若斯並利用柴比西夫多項式的根作為配置點，這種方法即稱為正交配置法(Orthogonal Collocation)。
[B] 巴諾夫—葛勒金法(Bubnov-Galerkin Method)
在配重殘值法中最有名的方法，可能就是巴諾夫--葛勒金法。在這種方法中，定義配重函數
為基礎函數uj(x)對未定係數Cj的導函數，以數學式表示為
(10-1.10)
在這種方法中，uj(x)為一組基礎函數之一，使得在所考慮控制體積內，任何函數都可以利用這組基礎函數表示之，
。因此，當N趨近於無窮大時，近似解在所考慮的空間內即能代表正確解。
[C] 高斯—李根德最小平方法(Gauss-Legendre Least Square Method)
高斯—李根德最小平方法的基本觀念，是利用殘值的平方積分式為零，找出讓殘值最小化的未定係數。
。亦即，
由以上關係式，可以建立N組代數方程式。這種方法極具有數學意義，理論上是強迫讓誤差值得平方和達最小化。由上一方程式，比較配重函數之定義，可知這種方法的配重函數
為
(10-1.11)
最小平方法在許多工程應用上常被使用，但由於使用殘值的平方積分式，也很容易產生很困擾的代數方程式。
[D] 慣量法(Method of Meoments)
慣量法最先被應用在非線性滲透(Non-linear Diffusion)及流體力學中的層流邊界層(Laminar Boundary Layer)問題的解析上，基本上要求殘值的連續階次慣量為零。亦即配重函數為
(10-1.12)
第一階近似慣量法的結果，與考慮整個區間視為一個子區間的區間法相同，這種近似法在層流邊界層問題的解析上亦稱為馮卡門及薄豪森近似法(von
Karman & Pohlhausen Similarity)。
[E] 副區間法(Subdomain Method)
副區間法是將原來的區間(0，1)分割成N個副區間，並定義配重函數為
(10-1.13)
例10-1
配重殘值法
試利用各種不同的配重殘值法解邊界值常微分方程式(BVP-ODE)
註：正確解為
[解]
[A] N＝1；
；殘值表示式為
根據MWR的基本原理：
，以下以各種不同配重函數，建立未定係數C2的值。
(i) 配置法：配重函數為W＝δ(x－x1)
令配置點x1＝0.5﹔得R＝2 C2－0.5＝0，即C2＝0.25，近似函數為
若令x1＝0.6﹔則得C2＝0.3。近似函數為
注意：配置法所得結果會因配置點位置的選擇而異。
(ii) 巴諾夫－葛勒金法：配重函數為
故
； 
。
(iii) 最小平方法：配重函數為
故
； 
。
(iv) 慣量法：配重函數為
得到
； 
。
(v) 副區間法：配重函數為W(x)＝1；0＜x＜1
得到
； 
。
[B] N＝2；考慮二階近似函數，即取N=2﹔由方程式(10-1.5)得二階近似解為
代回原微分方程式，得到殘值為
(i) 配置法：配重函數為
；j＝1，2
令配置點為
及
﹔則得配置點殘值為
解之，得
，且
，代回原近似函數方程式，得到微分方程式的二階近似解為
與正確解相同。
(ii) 巴諾夫－葛勒金法：配重函數為
配重積分式可以分成兩組
解以上二聯立方程式，得到
，
。亦即近似解為
亦與正確解完全相同。
表10-1.1 配置法與葛勒金法準確度之比較
|
x |
正確解 |
N=1 |
N=2 |
||
|
配置法* |
葛勒金法 |
配置法 |
葛勒金法 |
||
|
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 |
0.9168 0.8347 0.7545 0.6773 0.6042 0.5360 0.4738 0.4187 0.3715 |
0.908 0.822 0.742 0.668 0.600 0.538 0.482 0.432 0.388 |
0.908 0.822 0.742 0.668 0.600 0.538 0.482 0.432 0.388 |
0.9168 0.8347 0.7545 0.6773 0.6042 0.5360 0.4738 0.4187 0.3715 |
0.9168 0.8347 0.7545 0.6773 0.6042 0.5360 0.4738 0.4187 0.3715 |
＊ 配置法N＝1時，配置點選用x＝0.6﹔若配置點選用x＝0.5，則結果略差。
由本例可發現各種配重殘值法均能很快地獲得此問題的近似解。但一般而言較低階的近似解，以葛勒金法能獲得較佳的結果。高階近似解則各方法所得結果都相當準確，但由於配置法只需解一組代數聯立方程式，即可求得係數
，而其他方法則均須先作積分處理，因此，N值較大時，配置法最具使用潛力，而且程式設計也最容易。因此，本章其餘各節將著重於各種數學配置法的介紹。
第二節 配重殘值法的應用
為了說明配重殘值法的基本使用方法，其次再以圓柱狀觸媒粒子的徑向質傳滲透及等溫n次非可逆反應模式為例，說明配重殘值法的應用方法。
(10-2.1)
BC1 x＝0時， 
＝0
BC2 x＝1時， y＝1
其中y＝C／Cs為無因次濃度，x＝r／R為無因次半徑，
稱為希笠數(Thiele modulus)。觸媒粒子的有效度係數為
(10-2.2)
[A] 最低階配重殘值法y1(x)
當n=1時，考慮滿足兩個邊界條件即方程式(10-2.1)的最低階近似函數為y1(x)=1+a1(1-x2)，代入方程式(10-2.1)得到殘值表示式為
(10-2.3)
根據配重殘值法的定義
或 
(10-2.4)
配置法：配重函數為
；令x1=1/2，代入方程式(10-2.4)得到
，故得到未定係數a1值及一階近似函數y1(x)為
(10-2.5)
觸媒粒子的有效度係數為
。
葛勒金法：配重函數為
，或
；代入方程式(10-2.4)得到
，積分後得到未定係數a1值及一階近似函數y1(x)為
(10-2.6)
觸媒粒子的有效度係數為
。
副區間法：配重函數為
；
。代入方程式(10-2.4)得到
，積分後得到未定係數a1值及一階近似函數y1(x)為
(10-2.7)
觸媒粒子的有效度係數為
。
慣量法：配重函數為
；代入方程式(10-2.4)得到
，積分後得到未定係數a1值及一階近似函數y1(x)為
(10-2.8)
觸媒粒子的有效度係數為
。
最小平方法：配重函數為殘值對未定係數之微分，
，或
；將配重函數代入方程式(10-2.4)得到
，積分後得到未定係數a1值及一階近似函數y1(x)為
(10-2.9)
觸媒粒子的有效度係數為
。
方法比較: 當
時，比較各種方法所得結果，如表10-2.1所示。可以發現葛勒金法結果最佳。
表10-2.1 一階配重殘值法所得結果之比較
[B] 高階配重殘值法yN(x)
當n=1時，考慮滿足兩個邊界條件及方程式(10-2.1)的N階近似函數為
，其中
為滿足微分方程式邊界條件的函數，
為滿足微分方程式均勻邊界條件(homogeneous boundary conditions)的函數。就方程式(10-2.1)而言，
(10-2.10)
或將N階近似函數寫成下式，其中令
(10-2.11)
方程式(10-2.1)可以改寫成
(10-2.12)
將方程式(10-2.11)代入方程式(10-2.12)，經整理後，得到
(10-2.13)
根據配重殘值法的定義
即 
或再寫成u的函數為
，j=1,2,…,N (10-2.14)
將方程式(10-2.13)代入方程式(10-2.14)，積分結果可以寫成以下的矩陣形式，
(10-2.15)
其中矩陣A，B，及C分別因配重函數Wj(u)的定義而會有改變。解方程式(10-2.15)得到未定係數陣列為
(10-2.16)
將此結果代入方程式(10-2.11)，即可得到微分方程式的近似解。同時，也可以求得觸媒的有效度係數為
(10-2.17)
配置法：配重函數為
；代入方程式(10-2.13)及方程式(10-2.14)，得到
或
與方程式(10-2.15)比較，得到矩陣A，B，及C分別為
(10-2.18)
若考慮使用等間距配置點，由於u0=0，uN+1=1，則
，代回方程式(10-2.17)，得到
(10-2.19)
葛勒金法：配重函數為
；代入方程式(10-2.13)及方程式(10-2.14)，得到
，積分展開後，得到矩陣A，B，及C分別為
(10-2.20)
副區間法：配重函數為
；其他區間
。代入方程式(10-2.14)，得到
。將方程式(10-2.13)代入，積分後得到得到矩陣A，B，及C分別為
(10-2.21)
若考慮使用等間距區間，由於u0=0，uN+1=1，則
。
慣量法：配重函數為
；代入方程式(10-2.13)及方程式(10-2.14)，得到
積分後得到矩陣A，B，及C分別為
(10-2.22)
最小平方法：配重函數為殘值對未定係數之微分，
，或
；將此配重函數代入方程式(10-2.14)，並將方程式(10-2.13)代入，積分後可以得到矩陣A，B，及C。由於積分冗長，不予列出。
計算方法與使用策略
配重殘值法的計算方法，可以歸納如下:
(1) 利用一個多項式
作為微分方程式L{y,u}＝0的近似解；
(2) 代入微分方程式L{y,u}＝0，求出殘值的表示式R＝L{yN,u
}-L{y,u}﹔
(3) 選擇使用的配重殘值法，定義配重函數
；
(4) 使殘值的配重積分值為零
(5) 建立N個未定係數陣列a的聯立方程式；
(6) 利用數值方法解聯立方程式，找出近似多項式的未定係數陣列a﹔
(7) 代回原近似多項式
，建立微分方程式的近似解﹔
(8) 利用適當內插法，計算微分方程式在特定位置的解。
以上所述之計算邏輯如圖10-2.12所示；利用Excel或高斯消去法，都可以很輕易地建立所需的程式。利用Excel所建立程式如圖10-2.2所示。程式編寫細節可參考本書所附光碟。
圖10-2.1 配重殘值法的邏輯圖
圖10-2.2 (a) 利用Excel編寫配重殘值法的範例，配置法。將微分方程式的求解，變成簡單的矩陣運算。利用Excel即可輕易地完成任務，N=3所得結果與正確解相當一致。
圖10-2.2 (b)
利用Excel編寫配重殘值法的範例，配置法。
N=5所得結果。
第三節 數學配置法
上一節中，我們利用簡單的線性問題為例，說明了各種MWR配重殘值法及其特性，本節中則將針對配置法作較詳盡的討論。使用配置法的時候，近似函數及配置點的選擇方法大致可分成三類：
(i) 內部配置法(Interior Collocation)
所選用的近似函數滿足邊界條件，配置點選擇位於微分方程式所考慮的區間內，以決定未定係數。
(ii) 邊界配置法(Boundary Collocation)
所選用的近似函數滿足原微分方程式，而配置點則在邊界上選取，以決定未定係數。
(iii) 混和配置法(Mixed Collocation)
所選用的近似函數為任意函數，其未定係數由所考慮區間內及邊界上選定的配置點來決定。
配置法常見於邊界值常微分方程式、特徵值問題(Eigenvalue Problem)及偏微分的求解。以下再利用幾個例題加以說明。
例10-2
一次元熱傳導(K≠常數)
一平板厚度為L，兩側溫度分別為T0及T1。其能量平衡方程式為
﹔0＜x＜L (10-3.1)
BC1. T(0)＝T0
BC2. T(L)＝T1
此平板的熱傳導率與溫度關係為k(T)＝k0＋α(T－T0)，試求其溫度分布。
[解] 令無因次溫度及無因次長度分別為
則原微分方程式可以改寫成
﹔0＜ξ＜1 (10-3.2)
θ(0)＝1
θ(1)＝1
其中a＝α(T－T0)／k0﹔為了簡化說明起見，令a＝1，則滿足邊界條件的基礎函數，可以寫成
(10-3.3)
若取N＝2，其殘值為
(10-3.4)
配置點取為
及
，代入殘值表示式，並讓殘值在配置點上為0，則得
C1＝－0.5992
C2＝0.1916
故得θ的二階近似解為 
或溫度分布近似解為
例10-3
非線性微分方程式
試利用數學配置法解非線性微分方程式。
﹔0＜θ＜1 (10-3.5)
BC1 θ(0)＝0
BC2 θ(1)＝0
[解] 滿足邊界條件的基礎函數可為
(10-3.6)
考慮最簡單的情況，N＝1，即
。則殘值為
(10-3.7)
取配置點為x＝0.5﹔令R在配置點上為0，得到
利用牛頓法求得上式的解為C1＝0.44712。故θ1為
與正確解之比較如下表，結果顯示用最簡單的近似解，效果還不錯呢！
|
x |
0.1 |
0.2 |
0.3 |
0.4 |
0.5 |
|
正確解 |
-0.0414 |
-0.0733 |
-0.0958 |
-0.1092 |
-0.1137 |
|
θ1 |
-0.0402 |
-0.0715 |
-0.0939 |
-0.1073 |
-0.1178 |
|
誤差% |
2.8 |
2.4 |
1.99 |
1.73 |
1.69 |
第四節 正交配置法
配置法使用方便，但其準確度卻決定於基礎函數的選擇、N值的大小及配置點的選擇。1938年蘭若斯(Lanczos)提議利用正交級數(Orthogonal Polynomials)作為基礎函數，並以正交級數的根作為配置點﹔後來再經過魏德森及史都華(Villadsen and Stewart)的改進，目前已成為被使用得最廣泛的配置法，通稱為「正交配置法」(Orthogonal Collocation Method)。
正交配置法基本上是以正交級數作為基礎函數﹔如果基礎函數的最高次數為N，則利用N＋1次的正交級數的根作為配置點。由於正交級數具有許多易於處理的數學特性，因此，正交配置法基本原理雖與前述配置法完全一致，但使用上卻可善用正交級數的特性，使處理方法更為簡化，更適於利用計算機來求解。
假設一微分方程式的近似解，可以利用一組正交級數P(x)的線性和來表示，
(10-4.1)
在以上的方程式中總共有N＋2個未定係數
，但近似解yN的下註標只寫成N，是由於邊界值常微分方程式(BVP ODE)的兩個邊界條件已經可以決定兩個
值，只需要再利用N個內部配置點即可以決定所有
值，故寫成yN；其中，N代表所需的內部配置點數。
正交級數P(x)可以寫成級數的一般表示式
(10-4.2)
其正交性質為
﹔k＝0，1，2，……，m－1 (10-4.3)
其中(a，b)為邊界值常微分方程式的定義範圍，通常都先利用座標變換轉移成(0，1)，以便處理。
將方程式(10-4.2)代入方程式(10-4.1)中，可將基礎函數yN改寫成以下的方式，以利計算處理：
(10-4.4)
假設可以利用適當方法得到配置點，則在配置點xj位置上的函數值yN為
﹔j＝1，2，……，N＋2 (10-4.5)
x1＝0， xN＋2＝1
而在配置點xj位置上，yN的一次及二次導函數，亦可分別由方程式(10-4.4)對x作微分，得到
(10-4.6)
(10-4.7)
令矩陣
，
及
分別為
(10-4.8)
令陣列
，
，則方程式(10-4.5)，(10-4.6)及(10-4.7)分別改寫成矩陣符號後，得到
(10-4.9)
(10-4.10)
(10-4.11)
根據方程式(10-4.8)可知，若已知配置點位置，則矩陣
，
及
都可以很快的計算出來。而由方程式(10-4.9)，可以將未定義的係數陣列
表示成配置點上的函數值
的關係
(10-4.12)
將上式代回方程式(10-4.10)及方程式(10-4.11)中，可以分別得到
(10-4.13)
(10-4.14)
根據方程式(10-4.8)可知，若選定一組配置點
，則
、
及
均可直接求得，而上二方程式中
、
也都可以利用矩陣操作很快的求得。因此，可將微分方程式中的一次及二次導函數項，都表示成配置點位置的函數值
的線性函數，而將原微分方程式轉化成一組聯立代數方程式。矩陣A及B即稱為微分操作矩陣。利用微分操作矩陣，可以將微分方程式改寫成函數值
的線性函數，亦即變成一組線性聯立方程式。利用高斯消去法或適當數值方法，即可求得在配置點處的函數值
。
以上所說明的正交配置法基本邏輯，可以歸納如下：
(1)
利用正交多項式的解，得到一組配置點
；
(2)
利用方程式(10-4.8)及配置點
，建立矩陣
，
及
；
(3)
建立一次微分操作矩陣
，其中
；
(4)
建立二次微分操作矩陣
，其中
；
(5)
利用邊界條件建立y1及yN+2的表示式(在本章第五節中詳細說明)；
(6)
將微分操作矩陣代入原微分方程式，建立函數值
的聯立方程式；
(7)
利用高斯消去法或其他適當數值方法，求解函數值
的聯立方程式；
(8)
利用適當內插法，建立特定點之函數值y。
例10-4
利用配置法解邊界值常微分方程式
BC1 y(0)＝0
BC2 y(1)＝1
[解] 為了說明起見，我們只考慮最簡單的情況N＝1，則N＋2個配置點分別為
，
及
。根據方程式(10-4.8)得
，
及
分別為
矩陣
的反矩陣
可利用高斯約旦法求得為
因此，微分操作矩陣可分別根據其定義求得。
代回原微分方程式得到
﹔j＝1，2，……，N＋2
由於y1＝0，yN＋2＝y3＝1，只剩下y2未知，因此，只考慮N＝1且j＝2的情況，由上式得
將y1及y3代入，得
，由與y2＞0，因此，可以得到y2＝0.57916。
利用這種處理方式，我們可獲得各配置點位置的y值。若需其他位置的y值，則可利用本書第三章所介紹的插值法求之。需要求y的積分值，則可利用本書第六章所介紹的高斯積分法(Gauss Quadrature)求之。
第五節 正交配置法邊界條件的處理
微分方程式的邊界條件通常可分為三大類，本節將以實例來說明使用正交配置法時，邊界條件的處理方法。
考慮一常微分方程式為例
(10-5.1)
其正確解為
(10-5.2)
利用上一節所介紹的處理方法，方程式(10-5.1)可以被改寫成
﹔j＝1，2，……，N+1 (10-5.3)
或將y1及yN+1移到方程式右側，改寫成
(10-5.4)
其中δji稱為Kronecker Delta，定義為
(10-5.5)
方程式(10-5.4)中，y1及yN＋2通常利用邊界條件決定之，以下分別利用各種不同的邊界條件加以說明。
[A] 邊界條件(I ─ I)
兩側邊界條件均為第一類(Dirichlet)邊界條件：
BC1 y(0)＝0﹔ x＝0
BC2 y(1)＝1﹔ x＝1 (10-5.6)
上式可改寫成y1＝0，yN＋2＝1，代入方程式(10-5.4)得到
或利用矩陣符號簡寫成
或進一步簡化，並寫成以下的矩陣方程式
(10-5.7)
其中
方程式(10-5.7)代表一組線性聯立方程式，可利用高斯消去法解之。
[B] 邊界條件(I ─ II)
原點側邊界條件為第一類邊界條件，在x=1位置之邊界條件為第二類邊界條件：
BC1 y(0)＝0﹔ x＝0
BC2 
﹔ x＝1 (10-5.8)
若以yi表示，則得到
由上式求得yN＋2為
代回方程式(10-5.4)，得到
j＝2，3，……，N＋1
或寫成
(10-5.9)
其中
[C] 邊界條件(I ─ III)
原點側邊界條件為第一類邊界條件，在x=1位置之邊界條件為第三類邊界條件：
BC1 y(0)＝0﹔ x＝0
BC2 
﹔ x＝1 (10-5.10)
BC2可利用一次導函數操作矩陣表示為
整理之，並利用BC1(y1＝0)。可以得到
代入方程式(10-5.4)得到
或以矩陣符號表示為
(10-5.11)
其中
[D] 邊界條件(II ─ I)
原點側邊界條件為第二類邊界條件，在x=1位置之邊界條件為第一類邊界條件：
BC1 
﹔ x＝0
BC2 y＝1﹔ x＝1 (10-5.12)
仿以上處理，得到
或
(10-5.13)
[E] 邊界條件(II ─ II)
兩側邊界條件均為第二類邊界條件：
BC1 
﹔ x＝0
BC2 
﹔ x＝1 (10-5.14)
仿以上處理，得y1及yN＋2分別為
或將上式中括號內的表示式以
表示，並寫成
代入方程式(10-5.4)得到
(10-5.15)
[F] 邊界條件(II ─ III)
原點側邊界條件為第二類邊界條件，在x=1位置之邊界條件為第三類邊界條件：
BC1 
﹔ x＝0
BC2 
﹔ x＝1 (10-5.16)
BC1及BC2分別可用微分操作矩陣寫成：
或整理成
解以上聯立方程式，得y1及yN＋2分別為
將上式代入方程式(10-5.4)得到
(10-5.17)
其中
[G] 邊界條件(III ─ III)
兩側邊界條件均為第三類邊界條件：
BC1 
﹔ x＝0
BC2 
﹔ x＝1 (10-5.18)
利用微分操作矩陣
，即
，可將邊界條件改寫成
整理成y1及yN＋2的聯立線性方程式為
由以上二方程式，可聯立解得y1及yN＋2分別為
或以符號代表，簡寫成
將上式代入方程式(10-5.4)，可以得到
(10-5.19)
其中
以上係根據常微分方程式
為例，所做的演導與說明，讀者針對特定微分方程式求解時，可根據相同作法作適當處理。
第六節 雅可必多項式與操作矩陣
正交配置法係採用正交多項式作為基礎函數，利用正交多項式的特性，簡化數學配置法的處理程序。常用的正交多項式包括李根德多項式(Legendre Polynomial)及本節所介紹的雅可必多項式(Jacobi Polynomial)等。以下針對雅可必多項式的特性及如何利用其特性撰寫程式，作簡要說明。
[Ａ] 雅可必多項式的基本定義
雅可必多項式是一種典型的正交多項式，其基本定義為
﹔j＝0，1，2，……，N－1 (10-6.1)
或表示成
(10-6.2)
有關雅可必多項式的詳細說明，請參考工程數學書籍。
[Ｂ] 雅可必多項式的羅德芮格表示法
以羅德芮格方程式(Rodrigues’ Formula)表示，雅可必多項式可以寫成
(10-6.3)
[C] 雅可必多項式的循序計算公式表示法
雅可必多項式的循序計算公式(Recurrence formula)為
pN = [x - gN(N,a,b)] pN-1
– hN(N,a,b) pN-2 ;
p0 = 0 (10-6.4)
; 
h1 = 0 ; 
[D] 常見的雅可必多項式
常見的雅可必多項式如下表所示：
|
a |
b |
N=0 |
N=1 |
N=2 |
N=3 |
|
0 |
0 |
1 |
2x-1 |
6x2-6x+1 |
20x3-30x2+12x-1 |
|
1 |
0 |
1 |
3x-1 |
10x2-8x+1 |
35x3-45x2+15x-1 |
|
2 |
0 |
1 |
4x-1 |
15x2-10x+1 |
56x3-63x2+18x-1 |
|
0 |
1 |
1 |
|
|
|
|
1 |
1 |
1 |
2x-1 |
6x2-6x+1 |
14x3-21x2+9x-1 |
|
2 |
1 |
1 |
|
|
|
例10-5
雅可必多項式之計算
試求
的值。
[解]
gN=1/2
[E] 雅可必多項式與配置點
由雅可必多項式的循序計算公式(Recurrence formula)pN = [x - gN(N,a,b)] pN-1
– hN(N,a,b) pN-2 ，
p0 = 0﹔利用牛頓法求雅可必多項式的根作為配置點xi時，假設已經知道x1、x2、…、xk共k個解，要求第k+1個解時，可以利用輔助函數
，將已知的解x1、x2、…、xk由多項式中先行去除掉。再利用牛頓法解輔助函數GN-K，得到迭代式為
xk+1,i= xk+1,i-1 –
d(x) (10-6.5)
其中
執行時，令xk+1,0 為略大於 xk 的數字，即xk+1,0 = xk + e﹔其中e為一個很小的值，例如10-4。雅可必多項式利用本節 [C] 的方程式計算，於計算求得xk+1後，讓k增加1，再重複以上步驟。利用這種方法所建立的副程式如下所示。
|
'
========================== '
ROOTS OF JACOBI POLYNOMIAL '
========================== ‘ Diff1(I) = G(I)
‘ Diff2(I) = H(I) ‘ AlphaPlusBeta
= Alpha + Beta AlphaMinusBeta
= Beta - Alpha AlphaBeta
= Beta * Alpha Diff1(1)
= (AlphaMinusBeta / (AlphaPlusBeta + 2) + 1) / 2 Diff2(1)
= 0 ‘
Calculate G(I) & H(I) If
N >= 2 Then For I = 2 To N
IM1 = I - 1 Z
= AlphaPlusBeta + 2 * IM1 Diff1(I)
= (AlphaPlusBeta * AlphaMinusBeta / Z / (Z + 2) + 1) / 2
If I = 2 Then
Diff2(I) = (AlphaPlusBeta + AlphaBeta + IM1) / Z / Z / (Z + 1)
Else
Z = Z * Z
Y = IM1 * (AlphaPlusBeta + IM1)
Y = Y * (AlphaBeta + Y)
Diff2(I) = Y / Z / (Z - 1)
End If Next I End
If '
Newton’s Method to find ROOT X
= 0 For
I = 1 To N Do
XD = 0
XN = 1
XE = 0
XM = 0
For J = 1 To N
XP = (Diff1(J) - X) * XN - Diff2(J) * XD
XQ = (Diff1(J) - X) * XM - Diff2(J) * XE - XN
XD = XN
XE = XM
XN = XP
XM = XQ
Next J
ZC = 1 Z
= XN / XM
If I > 1 Then
For J = 2 To I
ZC = ZC - Z / (X - Root(J - 1))
Next J
End If Z
= Z / ZC X
= X - Z Loop While (Abs(Z) >
0.000000001) Root(I) = X X = X + 0.0001 Next
I ‘
To include ROOT at x= 0 and x=1 If
(N1 = 1) Then For I = 1 To N J
= N + 1 - I
Root(J + 1) = Root(J) Next I Root(1) = 0 End
If If
(N2 = 1) Then Root(ND) = 1 ‘
Print ROOTs for the Polynomial Print Print
"** COLLOCATION POINTS:" Print For
I = 1 To ND Print Format(Root(I),
" 0.000000E+00 ") Next
I |
[F] 雅可必多項式的導函數
假設利用 [E] 中所介紹的牛頓法，可以得到雅可必多項式的N+1個根，xj，j=1,2,….,N+1，則可將雅可必多項式寫成
(10-6.6)
或改寫成循序計算式
p0(x)
= 1
pj(x)
= (x-xj) pj-1(x) ﹔ j = 1, 2, ……, N+1 (10-6.7)
將上式微分，可以得到
(10-6.8)
(10-6.9)
(10-6.10)
由於
，因此，利用以上方程式，可以很快的求得在各配置點位置上的導函數值。寫成副程式範例如下。
|
‘
Prepare Derivatives of the Polynomial For
I = 1 To ND X = Root(I) Diff1(I) = 1 Diff2(I) = 0 Diff3(I) = 0 For J = 1 To ND
If J <> I Then
Y = X - Root(J)
Diff3(I) = Y * Diff3(I) + 3 * Diff2(I)
Diff2(I) = Y * Diff2(I) + 2 * Diff1(I) Diff1(I) = Y *
Diff1(I)
End If Next J Next
I |
[G] 雅可必多項式與拉格蘭奇內插法
參考本書第三章數值內插法的說明，若已知(xi，yi)，i=1,2,…..,N+1，則yN(x)可以利用拉格蘭奇內插法(Lagrangian Interpolation)寫成
(10-6.11)
利用這種方法所建立的內插法副程式如下所示。
|
'
======================== '
LAGRANGIAN INTERPOLATION '
======================== POL
= 1 For
I = 1 To ND YVA = X1 - Root(I) XP(I) = 0 If YVA = 0 Then XP(I) = 1 POL = POL * YVA Next
I If
POL <> 0 Then For I = 1 To ND XP(I)
= POL / Diff1(I) / (X1 - Root(I)) Next I End
If |
[H] 雅可必多項式與微分操作矩陣
根據拉格蘭奇內插方程式(10-6.11)，取微分，可以得到y(x)的一次及二次導函數分別為
(10-6.12)
(10-6.13)
其中
其中雅可必多項式的導函數
可以利用方程式(10-6.8)、(10-6.9)及(10-6.10)求之。利用這種方法所建立的副程式如下所示。
|
Sub
DefMatrix(ID, ND, N1, N2, Diff1, Diff2, Diff3, Root, Vect, A, B) '
============================= '
DERIVATIVE OPERATION MATRICES '
============================= For
I = 1 To ND ID = 1 Call OpMatrix(I, ID, ND,
N1, N2, Diff1, Diff2, Diff3, Root, Vect) For J = 1 To ND
A(I, J) = Vect(J) Next J ID = 2 Call OpMatrix(I, ID, ND, N1,
N2, Diff1, Diff2, Diff3, Root, Vect) For J = 1 To ND
B(I, J) = Vect(J) Next J Next
I End
Sub Sub
OpMatrix(I, ID, ND, N1, N2, Diff1, Diff2, Diff3, Root, Vect) '
================ '
OPERATION MATRIX '
================ '
--ENTRY POINT For J = 1 To ND
If J = Index Then
If ID = 1 Then
Vect(J) = Diff2(Index) / Diff1(Index) / 2
Else
Vect(J) = Diff3(Index) / Diff1(Index) / 3
End If
Else
Y = Root(Index) - Root(J)
Vect(J) = Diff1(Index) / Diff1(J) / Y
If ID = 2 Then Vect(J) = Vect(J) * (Diff2(Index) / Diff1(Index) - 2 /
Y)
End If Next J End Sub |
[I] 雅可必多項式與高斯數值積分
積分操作可以採用本書第六章所介紹之高斯積分法。
(10-6.14)
其中，求出積分配重函數wi，即可計算積分結果。高斯雅可必配重函數為
W(x)
= xb(1-x)a
(10-6.15)
(10-6.16)
若x=０為積分配置點，則x1 = 1，否則x1 = 0。
若x=１為積分配置點，則x2 = 1，否則x2 = 0。
[J] 雷道與羅伯特數值積分
除了高斯積分法以外，雷道與羅伯特(Radau and Lobatto)積分法也是較常被使用的數值積分方法。其配重函數為
(10-6.17)
其中第一及第二式稱為雷道數值積分配重，第三式稱為羅伯特數值積分配重。
第七節 線性問題的程式規劃
根據以上諸節的介紹，我們可擬定出利用正交配置法解微分方程式的幾項步驟：
（i）
紙上作業：
l
由邊界條件縯導出y1及yN＋2的表示式(如本章第五節的說明)。
l
利用微分操作矩陣
及
，將微分方程式轉換成
並寫下
及
的表示方法。
（ii）
程式規劃：
l
決定使用的N值(通常選用4∼10間即可得到相當正確的解)。
l
決定使用的正交級數種類。
l
設計程式找出正交多項式PN(x)的根xj，j＝1，2，……，N，作為配置點。
l
設計程式求出微分操作矩陣
及
。
l
建立
及
，並將微分方程式改寫成聯立方程式。
l
利用高斯消去法解線性聯立方程式
，求得配置點上的函數值
。
l
利用內插法以適當x間距列印y值。
第八節 非線性常微分方程式
非線性常微分方程式基本上分成三大類：(1)線性常微分方程式，邊界條件非線性﹔(2)非線性常微分方程式，邊界條件為線性﹔(3)非線性常微分方程式，非線性邊界條件。但仿以上各節處理後，這三類非線性常微分方程式均可變成一組聯立的非線性代數方程式，利用本書所介紹的牛頓拉福森法或割線法即可仿以上各節求得
。以下利用實例說明之。
參考文獻
1.
Finlayson, B. A., “The Method of
Weighted Residuals and Variational Principles” Academic Press, (1972).
2.
Villadsen, J., and M. L. Michelsen;
“Solution of Differential Equation Models by Polynomial Models”, Prentice-Hall,
(1978).
3.
Finalyson, B. A., “Nonlinear Analysis
in Chemical Engineering”. (1980).
4.
Davis, M. E., “Numerical Methods and
Modeling for Chemical Engineers”, John Wiley, New York, (1984).
5.
Aris, R. “The Mathematical Theory of
Diffusion and Reaction in Permeable Catalysts”. Oxford, Clarendon Press (1975).
6.
Carberry, J. J., “Chemical and
Catalytic Reaction Engineering” McGraw-Hill, New York, (1976).
7.
Weisz, P. B., and J. S. Hicks, “The
Behavior of Porous Catalyst Particles in View of Internal Mass and Heat
Diffusion Effects”, Chem. Eng. Sci., 17, 265 (1962).
