1、第8章 常微分方程邊值問題的數(shù)值解法8.1 引 言第7章介紹了求解常微分方程初值問題的常用的數(shù)值方法；本章將介紹常微分方程的邊值問題的數(shù)值方法。只含邊界條件(boundary-value condition)作為定解條件的常微分方程求解問題稱為常微分方程的邊值問題(boundary-value problem). 為簡明起見，我們以二階邊值問題為例介紹常用的數(shù)值方法。一般的二階常微分方程邊值問題(boundary-value problems for second-order ordinary differential equations)為, (8.1.1)其邊界條件為下列三種情況之一：(1

2、) 第一類邊界條件 (the first-type boundary conditions)：(2) 第二類邊界條件 (the second-type boundary conditions)：(3) 第三類邊界條件 (the third-type boundary conditions)： 定理8.1.1 設(shè)(8.1.1)中的函數(shù)及其偏導(dǎo)數(shù), 在上連續(xù). 若(1) 對所有，有；(2) 存在常數(shù)，對所有，有,則邊值問題(8.1.1)有唯一解。推論 若線性邊值問題 (8.1.2)滿足(1) 和在上連續(xù)；(2) 在上, ,則邊值問題(8.1.1)有唯一解。求邊值問題的近似解，有三類基本方法：(1)

3、 差分法(difference method)，也就是用差商代替微分方程及邊界條件中的導(dǎo)數(shù)，最終化為代數(shù)方程求解；(2) 有限元法(finite element method)；(3) 把邊值問題轉(zhuǎn)化為初值問題，然后用求初值問題的方法求解。8.2 差分法8.2.1 一類特殊類型二階線性常微分方程的邊值問題的差分法設(shè)二階線性常微分方程的邊值問題為其中在上連續(xù)，且.用差分法解微分方程邊值問題的過程是：(i) 把求解區(qū)間分成若干個等距或不等距的小區(qū)間，稱之為單元；(ii) 構(gòu)造逼近微分方程邊值問題的差分格式. 構(gòu)造差分格式的方法有差分法, 積分插值法及變分插值法；本節(jié)采用差分法構(gòu)造差分格式；(iii

4、) 討論差分解存在的唯一性、收斂性及穩(wěn)定性；最后求解差分方程.現(xiàn)在來建立相應(yīng)于二階線性常微分方程的邊值問題(8.2.1), (8.2.2)的差分方程.( i ) 把區(qū)間等分，即得到區(qū)間的一個網(wǎng)格剖分：,其中分點(diǎn)，并稱之為網(wǎng)格節(jié)點(diǎn)(grid nodes)；步長.( ii ) 將二階常微分方程(8.2.2)在節(jié)點(diǎn)處離散化：在內(nèi)部節(jié)點(diǎn)處用數(shù)值微分公式 (8.2.3)代替方程(8.2.2)中，得, (8.2.4)其中.當(dāng)充分小時，略去式(8.2.4)中的，便得到方程(8.2.1)的近似方程, (8.2.5)其中, 分別是的近似值, 稱式(8.2.5)為差分方程(difference equation)

5、，而稱為差分方程(8.2.5)逼近方程(8.2.2)的截?cái)嗾`差(truncation error). 邊界條件(8.7.2)寫成 (8.2.6)于是方程(8.2.5), (8.2.6)合在一起就是關(guān)于個未知量，以及個方程式的線性方程組： (8.2.7)這個方程組就稱為逼近邊值問題(8.2.1), (8.2.2)的差分方程組(system of difference equations)或差分格式(difference scheme)，寫成矩陣形式. (8.2.8)用第2章介紹的解三對角方程組的追趕法求解差分方程組(8.2.7)或(8.2.8), 其解稱為邊值問題(8.2.1), (8.2.2)

6、的差分解(difference solution). 由于(8.2.5)是用二階中心差商代替方程(8.2.1)中的二階微商得到的，所以也稱式(8.2.7)為中心差分格式(centered-difference scheme).( iii ) 討論差分方程組(8.2.7)或(8.2.8)的解是否收斂到邊值問題(8.2.1), (8.2.2)的解，估計(jì)誤差.對于差分方程組(8.2.7)，我們自然關(guān)心它是否有唯一解；此外，當(dāng)網(wǎng)格無限加密，或當(dāng)時，差分解是否收斂到微分方程的解. 為此介紹下列極值原理：定理8.2.1 (極值原理) 設(shè)是給定的一組不全相等的數(shù)，設(shè) . (8.2.9)(1) 若, 則中非負(fù)

7、的最大值只能是或；(2) 若, 則中非正的最小值只能是或.證 只證(1) 的情形，而(2) 的情形可類似證明. 用反證法. 記，假設(shè), 且在中達(dá)到. 因?yàn)椴蝗嗟龋钥偪梢哉业侥硞€，使，而和中至少有一個是小于的. 此時因?yàn)?，所? 這與假設(shè)矛盾，故只能是或. 證畢！推論 差分方程組(8.2.7)或(8.2.8)的解存在且唯一.證明 只要證明齊次方程組 (8.2.10)只有零解就可以了. 由定理8.7.1知，上述齊次方程組的解的非負(fù)的最大值和非正的最小值只能是或. 而，于是 證畢！利用定理8.2.1還可以證明差分解的收斂性及誤差估計(jì). 這里只給出結(jié)果：定理8.2.2 設(shè)是差分方程組(8.2.7

8、)的解，而是邊值問題(8.2.1), (8.2.2)的解在上的值，其中. 則有 (8.2.11)其中.顯然當(dāng)時，. 這表明當(dāng)時，差分方程組(8.2.7)或(8.2.8)的解收斂到原邊值問題(8.7.1), (8.7.2)的解.例8.2.1 取步長，用差分法解邊值問題 并將結(jié)果與精確解進(jìn)行比較.解 因?yàn)椋? 由式(8.2.7)得差分格式 , , 其結(jié)果列于表8.2.1.表8.2.1準(zhǔn)確值010010.10. 0.20.20. 0.30.30. 0.40.40. 0.50.50. 0.60.60. 0.70.70. 0.80.80. 0.90.90. 0.101.000從表8.2.1可以看出, 差

9、分方法的計(jì)算結(jié)果的精度還是比較高的. 若要得到更精確的數(shù)值解，可用縮小步長的方法來實(shí)現(xiàn).8.2.2 一般二階線性常微分方程邊值問題的差分法對一般的二階微分方程邊值問題 (8.2.12)假定其解存在唯一. 為求解的近似值，類似于前面的做法，( i ) 把區(qū)間等分，即得到區(qū)間的一個網(wǎng)格剖分：,其中分點(diǎn)，步長.( ii ) 對式(8.2.12)中的二階導(dǎo)數(shù)仍用數(shù)值微分公式代替，而對一階導(dǎo)數(shù)，為了保證略去的逼近誤差為，則用3點(diǎn)數(shù)值微分公式；另外為了保證內(nèi)插，在2個端點(diǎn)所用的3點(diǎn)數(shù)值微分公式與內(nèi)網(wǎng)格點(diǎn)所用的公式不同，即 (8.2.13)略去誤差，并用的近似值代替，便得到差分方程組 (8.2.14)其中,

10、 是的近似值. 整理得 (8.2.15)解差分方程組(8.2.15)，便得邊值問題(8.2.12)的差分解.特別地, 若，則式(8.2.12)中的邊界條件是第一類邊值條件：此時方程組(7.7.16)為 (8.2.16)方程組(8.2.16)是三對角方程組，用第2章介紹的解三對角方程組的追趕法求解差分方程組(8.2.16)，便得邊值問題(8.2.12)的差分解.( iii ) 討論差分方程組(8.2.16)的解是否收斂到微分方程的解，估計(jì)誤差. 這里就不再詳細(xì)介紹.例8.2.2 取步長，用差分法求下列邊值問題的近似解，并將結(jié)果與精確解進(jìn)行比較.精確解是.解 因?yàn)椋? 由式(8.2.17)得差分格

11、式 , , 其結(jié)果列于表8.2.2.表8.2.2準(zhǔn)確值00-0.3-0.31/16-0.-0.22/16-0.-0.33/16-0.-0.44/16-0.-0.55/16-0.-0.66/16-0.-0.77/16-0.-0.8/2-0.-0.8.3 有限元法有限元法(finite element method)是求解微分方程定解問題的有效方法之一，它特別適用在幾何、物理上比較復(fù)雜的問題. 有限元法首先成功地應(yīng)用于結(jié)構(gòu)力學(xué)和固體力學(xué)，以后又應(yīng)用于流體力學(xué)、物理學(xué)和其他工程科學(xué). 為簡明起見，本節(jié)以線性兩點(diǎn)邊值問題為例介紹有限元法.考慮線性兩點(diǎn)邊值問題其中, . 此微分方程描述了長度為的可變交叉

12、截面(表示為)的橫梁在應(yīng)力和下的偏差.8.3.1 等價性定理 記, 引進(jìn)積分. (8.3.3)任取，就有一個積分值與之對應(yīng)，因此是一個泛函(functional)，即函數(shù)的函數(shù). 因?yàn)檫@里是的二次函數(shù)，因此稱為二次泛函.對泛函(8.3.3)有如下變分問題(variation problem)：求函數(shù)，使得對任意, 均有, (8.3.4)即在處達(dá)到極小, 并稱為變分問題(8.3.4)的解.可以證明：定理8.3.1(等價性定理) 是邊值問題(8.3.1), (8.3.2)的解的充分必要條件是使泛函在上達(dá)到極小，即是變分問題(8.3.4)在上的解.證 (充分性) 設(shè)是變分問題的解；即使泛函在上達(dá)到極

13、小，證明必是邊值問題(8.3.1), (8.3.2)的解.設(shè)是任意一個滿足的函數(shù)，則函數(shù)，其中為參數(shù). 因?yàn)槭沟眠_(dá)到極小，所以，即積分作為的函數(shù)，在處取極小值，故. (8.3.5)計(jì)算上式，得利用分部積分法計(jì)算積分代入式(8.3.6)，得因?yàn)槭侨我夂瘮?shù)，所以必有. (8.3.8)否則，若在上某點(diǎn)處有，不妨設(shè)，則由函數(shù)的連續(xù)性知，在包含的某一區(qū)間上有.作顯然，且，但，這與式(8.3.7)矛盾. 于是式(8.3.8)成立，即變分問題(8.3.4)的解滿足微分方程(8.3.1), 且故它是邊值問題(8.3.1), (8.3.2)的解.(必要性) 設(shè)是邊值問題(8.3.1), (8.3.2)的解，證明

14、是變分問題(8.3.4)的解；即：使泛函在上達(dá)到極小.因?yàn)闈M足方程(8.3.1)，所以. 設(shè)任意，則函數(shù)滿足條件，且. 于是因?yàn)椋援?dāng)時，, 即.只有當(dāng)時，. 這說明使泛函在上達(dá)到極小. 證畢！定理8.3.2 邊值問題(8.3.1), (8.3.2)存在唯一解.證明 用反證法. 若都是邊值問題(8.3.1), (8.3.2)的解，且不相等，則由定理8.3.1知，它們都使泛函在上達(dá)到極小，因而 且 ，矛盾！因此邊值問題(8.3.1), (8.3.2)的解是唯一的. 由邊值問題解的唯一性，不難推出邊值問題(8.3.1), (8.3.2)解的存在性(留給讀者自行推導(dǎo)). 8.3.2 有限元法等價性

15、定理說明，邊值問題(8.3.1), (8.3.2)的解可化為變分問題(8.3.4)的求解問題. 有限元法就是求變分問題近似解的一種有效方法. 下面給出其解題過程：第1步 對求解區(qū)間進(jìn)行網(wǎng)格剖分區(qū)間稱為單元，長度稱為步長，. 若，則稱上述網(wǎng)格剖分為均勻剖分. 給定剖分后，泛函(8.3.3)可以寫成. (8.3.9)第2步 構(gòu)造試探函數(shù)空間。為了計(jì)算積分(8.3.9)，最簡單的近似方法是將分段線性函數(shù)的集合作為試探函數(shù)空間。設(shè)分別是邊值問題(8.3.1), (8.3.2)的解在節(jié)點(diǎn)處的值，用分段線性插值 (8.3.10)近似，并稱式(8.3.10)為單元形狀函數(shù)(element shape fun

16、ction).為了將線性插值(8.3.10)標(biāo)準(zhǔn)化，令，則將變到軸上的單元. 記，則式(8.3.10)可寫成 (8.3.11)第3步 建立有限元方程組. 將式(8.3.10)代入泛函(8.3.9)，有.由式(8.3.11)知 (8.3.12)式(8.3.12)右端第1個求和號內(nèi)的第項(xiàng)(對應(yīng)第個單元)是關(guān)于的二次型，可以寫成, (8.3.13)其中，, (8.3.14)稱為單元剛度矩陣(element stiffness matrix)， (8.3.15)而式(8.3.12)的第2個求和號內(nèi)的第項(xiàng)可以寫成 (8.3.16)其中，于是式(8.3.12)中求和號內(nèi)的項(xiàng)可以寫成. (8.3.17)再令

17、及矩陣則. 于是式(8.3.17)又可寫成 . (8.3.18)把式(8.3.18)代入式(8.3.12)右端求和號內(nèi)，得. (8.3.19)記，(8.3.20), (8.3.21)則式(8.3.19)化為 (8.3.22)并稱為總剛度矩陣(total stiffness matrix)，稱為右端向量. 因?yàn)槭鞘谷O小值的函數(shù)，所求的自然使式(8.3.22)的右邊取極小值，所以應(yīng)有. (8.3.23)記, 則式(8.8.23)為或 (8.3.24)得方程組. (8.3.25)因?yàn)槭且阎?，不能任意選取，所以不能要求式(8.3.23)對也成立. 因此方程組(8.3.24)或(8.3.25)中應(yīng)當(dāng)

18、去掉首末2個方程，并把其他方程中含有的改為已知量，所得方程組就是未知量滿足的代數(shù)方程組. (8.3.26)方程組(8.3.25)或(8.3.26)稱為有限元方程組(finite element system of equations). 用第2章介紹的解三對角方程組的追趕法求解有限元方程組(8.3.26)，可解出，即得變分問題(8.3.4)的近似解，也就是邊值問題(8.3.1), (8.3.2)解的近似值. 這種方法稱為有限單元法(finite element method)或有限元法.例8.3.1 用有限元法求下列邊值問題的近似解，并將結(jié)果與精確解進(jìn)行比較.取，精確解是.解 因?yàn)椋? 由式(

19、8.3.26)得有限元方程組 其結(jié)果列于表8.3.1.表8.3.1準(zhǔn)確值000010.2-0.1644-0.1620.4-0.2443-0.2430.6-0.2434-0.2440.8-0.1620-0.165100上面雖然是對邊值問題(8.3.1), (8.3.2)介紹的有限元解法，但其解題步驟對一般的微分方程定解問題也是適用的. 對所給微分方程定解問題，首先找出相應(yīng)微分方程的變分形式，然后進(jìn)行下列步驟：第1步 對定義區(qū)域(或定義區(qū)間)進(jìn)行網(wǎng)格剖分, 將定義區(qū)域(或定義區(qū)間)剖分為若干個小單元(一維是區(qū)間；多維是區(qū)域，如矩形、三角形等)；第2步 構(gòu)造試探空間; 即選擇適當(dāng)?shù)牟逯岛瘮?shù)類.第3步

20、 建立有限元方程組. 計(jì)算單元剛度矩陣及右端向量，再形成總剛度矩陣及總右端項(xiàng)，寫出有限元方程組；結(jié)合定解條件，求解有限元方程組.注 從形式上看，用有限元法解微分方程定解問題很繁，但是有限元法的求解步驟規(guī)范，便于在計(jì)算機(jī)上實(shí)現(xiàn)，并且總剛度矩陣是三對角對稱正定矩陣，因此有限元方程組可用第2章介紹的解三對角方程組的追趕法求解. 有限元法最主要的優(yōu)點(diǎn)是可以處理相當(dāng)復(fù)雜的區(qū)域上的邊值問題。8.4 打靶法解常微分方程邊值問題的打靶法(shooting method)，也稱為嘗試法，其基本思想是把邊值問題轉(zhuǎn)化為初值問題來解：首先作出一些只滿足一端邊值條件的解，然后再從這些解中找出適合另一端邊值條件的解. 下

21、面以二階常微分方程帶第一類邊界條件的邊值問題為例介紹常微分方程邊值問題的打靶法.7.6節(jié)曾介紹過二階常微分方程初值問題(7.6.10) 的求解方法. 將上式中的改為，改為，得 (8.4.3)設(shè)初值問題(8.4.3)的解為，顯然依賴于，即. 為了求解邊值問題(8.4.1), (8.4.2)，必須求，使之滿足. 下面介紹2種方法來求.方法1 根據(jù)實(shí)際問題情況選一個，解初值問題 (8.4.4)得到的解仍記為. 若或(為事先給定的精度)，則就是邊值問題(8.4.1), (8.4.2)的解. 否則，根據(jù)與的誤差，將修改為；例如取, 再解初值問題 (8.4.5)得到解. 若滿足或，則就是邊值問題(8.4.

22、1), (8.4.2)的解. 否則，再將適當(dāng)修改為; 例如可用作線性插值求：,然后解初值問題 (8.4.6)的解. 如此繼續(xù)下去，直到達(dá)到精度要求為止.方法2 求另一種方法是求函數(shù)的一個零點(diǎn). 對于每一個自變量，通過解初值問題(8.4.4)，可解出. 計(jì)算, (8.4.7)然后采用第3章介紹的求方程根的方法求的零點(diǎn). 首先尋找，使，則在區(qū)間或內(nèi)至少有的一個零點(diǎn). 然后可用二分法求. 也可用Newton迭代公式求的近似值.幾何解釋：微分方程邊值問題(8.4.1), (8.4.2)的解是一條通過兩點(diǎn)的曲線(如圖8.4.1). 初值問題(8.4.4)的解是一條通過點(diǎn)、斜率為的曲線(如圖8.4.1).

23、 初值問題(8.4.5)的解是一條通過點(diǎn)、斜率為的曲線(如圖8.4.1). 這有點(diǎn)象射擊者從定點(diǎn)向目標(biāo)射擊一樣. 根據(jù)經(jīng)驗(yàn)以某一角度(斜率為)試射一次. 如果與目標(biāo)相差太大，調(diào)整射擊角度(斜率為)，再進(jìn)行射擊. 如此繼續(xù)進(jìn)行下去，直到擊中目標(biāo)，或擊中的點(diǎn)與的誤差在允許的范圍之內(nèi)。圖8.4.1參考文獻(xiàn)1提供的資料還討論了選擇初始值的重要性，并介紹了多重打靶法，這里就不作詳細(xì)介紹，有興趣的讀者可參看相關(guān)資料.8.5 算法程序8.5.1 用差分法解二階線性常微分方程的邊值問題%用差分法解二階線性常微分方程的邊值問題%a,b是區(qū)間, Step是步長, Alpha, Beta是初值% f, q分別是式(

24、8.2.1)中的f (x), q (x)function Y = DiffMethod(f, q, a, b, Step, Alpha, Beta)N = (b-a) / Step;X = a : Step : b;A = zeros( N-1 );for i = 1 : N-1 A(i,i) = -1 * ( 2+feval(q, X(i+1) * Step2 ); if i = N-1 A(i,i+1) = 1; A(i+1,i) = 1; endendB(1) = Step2 * feval(f, X(2) - Alpha;B(2:N-2) = Step2 * feval(f, X(3:

25、N-1);B(N-1) = Step2 * feval(f, X(N) - Beta;B = B;L,U = lu( A );Y = U ( L B ) ;for i = 1 : length(Y) sprintf(%10.7f,Y(i)endend例8.5.1 取，用差分法求下列邊值問題的近似解解 在MATLAB命令窗口輸入：f = inline(-4.*x); q= inline(4); DiffMethod(f, q, 0, 1, 0.1, 0, 2)回車，可的結(jié)果。8.5.2 用差分法解一般二階常微分方程的邊值問題%用差分法解一般二階常微分方程的邊值問題% a,b是區(qū)間, Step是步

26、長, Alpha, Beta是初值% f, p, q 分別是式(8.2.12)中的f (x), p (x), q (x) function Y = DiffMethod_2(f, p, q, a, b, Step, Alpha, Beta)N = (b-a) / Step;X = a : Step : b;A = zeros( N-1 );A(1, 1) = -2 * ( 2 - Step2 * feval(p, X(2) );A(1,2) = 2 + Step * feval(f, X(3);for i = 2 : N-2 A(i,i) = -2 * ( 2 - feval(p, X(i+1

27、) * Step2 ); A(i,i-1)= 2 - Step * feval(f, X(i+1); A(i-1,i) = 2 + Step * feval(f, X(i+1);endA(N-1,N-2) = 2 - Step * feval(f, X(N);A(N-1,N-1) = -2 * ( 2 - feval(p, X(N) * Step2 );B(1) = 2*Step2 * feval(q, X(2) - ( 2 - Step * feval(f, X(2) ) * Alpha;B(2:N-2) = 2 * Step2 * feval(q, X(3:N-1);B(N-1) = 2*

28、Step2 * feval(q, X(N) - ( 2 + Step * feval(f, X(N) ) * Beta;B = B;L,U = lu( A );Y = U ( L B ) ;for i = 1 : length(Y) sprintf(%10.7f,Y(i)endend例8.5.2 取，用差分法求下列邊值問題的近似解解 在MATLAB命令窗口輸入：f = inline(cos(x); p = inline(-1); q = inline(-2);DiffMethod_2(f, p, q, 0, pi/2, pi/16, -0.3, -0.1)回車，可得結(jié)果。8.5.3 用有限元法

29、解二階常微分方程的邊值問題%用有限元法解二階常微分方程的邊值問題%a,b是區(qū)間, h是步長, Alpha, Beta是初值% f, p, q 分別是式(8.3.1)中的f (x), p (x), q (x) function Y = FiniElem(f, p, q, a , b, Step, Alpha, Beta)N = length(Step) - 1;X(1) = a;for i = 1 : N+1 X(i+1) = X(i) + Step(i);endsyms t;ff = -1/Step(N+1) * feval(f, X(N)+t*Step(N+1) + . Step(N+1)*

30、 feval(p, X(N)+t*Step(N+1) * (1-t)*t;Knnn = int(ff, t, 0, 1);syms t;ff = -1/Step(2) * feval(f, X(1)+t*Step(2) + . Step(2)*feval(p, X(1)+t*Step(2)*(1-t)*t;K101 = int(ff, t, 0, 1);for i = 2 : N syms t; ff = Step(i+1) * feval(q, X(i)+t*Step(i+1) * (1-t); B1(i) = int(ff, t, 0, 1);endfor i = 1 : N-1 syms t; ff = Step(i+1) * feval(q, X(i)+t*Step(i+1) * t; B2(i) = int(ff, t, 0, 1);endB(1) = B2(1) + B1(2) - Alpha*K101;B(N-1) = B2(N-1) + B1(N) - Beta*Knnn;for i = 2: N-2 B(i) = B2(i) + B1(i+1);endB=B;A=zeros(N-1);for i