第8章 常微分方程邊值問(wèn)題的數(shù)值解法

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

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

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

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

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

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

7、的最大值只能是或；(2) 若, 則中非正的最小值只能是或.證 只證(1) 的情形，而(2) 的情形可類(lèi)似證明. 用反證法. 記，假設(shè), 且在中達(dá)到. 因?yàn)椴蝗嗟?，所以總可以找到某個(gè)，使，而和中至少有一個(gè)是小于的. 此時(shí)因?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、)的解，而是邊值問(wèn)題(8.2.1), (8.2.2)的解在上的值，其中. 則有 (8.2.11)其中.顯然當(dāng)時(shí)，. 這表明當(dāng)時(shí)，差分方程組(8.2.7)或(8.2.8)的解收斂到原邊值問(wèn)題(8.7.1), (8.7.2)的解.例8.2.1 取步長(zhǎng)，用差分法解邊值問(wèn)題 并將結(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ù)值解，可用縮小步長(zhǎng)的方法來(lái)實(shí)現(xiàn).8.2.2 一般二階線(xiàn)性常微分方程邊值問(wèn)題的差分法對(duì)一般的二階微分方程邊值問(wèn)題 (8.2.12)假定其解存在唯一. 為求解的近似值，類(lèi)似于前面的做法，( i ) 把區(qū)間等分，即得到區(qū)間的一個(gè)網(wǎng)格剖分：,其中分點(diǎn)，步長(zhǎng).( ii ) 對(duì)式(8.2.12)中的二階導(dǎo)數(shù)仍用數(shù)值微分公式代替，而對(duì)一階導(dǎo)數(shù)，為了保證略去的逼近誤差為，則用3點(diǎn)數(shù)值微分公式；另外為了保證內(nèi)插，在2個(gè)端點(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)，便得邊值問(wèn)題(8.2.12)的差分解.特別地, 若，則式(8.2.12)中的邊界條件是第一類(lèi)邊值條件：此時(shí)方程組(7.7.16)為 (8.2.16)方程組(8.2.16)是三對(duì)角方程組，用第2章介紹的解三對(duì)角方程組的追趕法求解差分方程組(8.2.16)，便得邊值問(wèn)題(8.2.12)的差分解.( iii ) 討論差分方程組(8.2.16)的解是否收斂到微分方程的解，估計(jì)誤差. 這里就不再詳細(xì)介紹.例8.2.2 取步長(zhǎng)，用差分法求下列邊值問(wèn)題的近似解，并將結(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)是求解微分方程定解問(wèn)題的有效方法之一，它特別適用在幾何、物理上比較復(fù)雜的問(wèn)題. 有限元法首先成功地應(yīng)用于結(jié)構(gòu)力學(xué)和固體力學(xué)，以后又應(yīng)用于流體力學(xué)、物理學(xué)和其他工程科學(xué). 為簡(jiǎn)明起見(jiàn)，本節(jié)以線(xiàn)性?xún)牲c(diǎn)邊值問(wèn)題為例介紹有限元法.考慮線(xiàn)性?xún)牲c(diǎn)邊值問(wèn)題其中, . 此微分方程描述了長(zhǎng)度為的可變交叉

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

13、小，證明必是邊值問(wèn)題(8.3.1), (8.3.2)的解.設(shè)是任意一個(gè)滿(mǎn)足的函數(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)成立，即變分問(wèn)題(8.3.4)的解滿(mǎn)足微分方程(8.3.1), 且故它是邊值問(wèn)題(8.3.1), (8.3.2)的解.(必要性) 設(shè)是邊值問(wèn)題(8.3.1), (8.3.2)的解，證明

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

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

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

17、及矩陣則. 于是式(8.3.17)又可寫(xiě)成 . (8.3.18)把式(8.3.18)代入式(8.3.12)右端求和號(hào)內(nèi)，得. (8.3.19)記，(8.3.20), (8.3.21)則式(8.3.19)化為 (8.3.22)并稱(chēng)為總剛度矩陣(total stiffness matrix)，稱(chēng)為右端向量. 因?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)對(duì)也成立. 因此方程組(8.3.24)或(8.3.25)中應(yīng)當(dāng)

18、去掉首末2個(gè)方程，并把其他方程中含有的改為已知量，所得方程組就是未知量滿(mǎn)足的代數(shù)方程組. (8.3.26)方程組(8.3.25)或(8.3.26)稱(chēng)為有限元方程組(finite element system of equations). 用第2章介紹的解三對(duì)角方程組的追趕法求解有限元方程組(8.3.26)，可解出，即得變分問(wèn)題(8.3.4)的近似解，也就是邊值問(wèn)題(8.3.1), (8.3.2)解的近似值. 這種方法稱(chēng)為有限單元法(finite element method)或有限元法.例8.3.1 用有限元法求下列邊值問(wèn)題的近似解，并將結(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上面雖然是對(duì)邊值問(wèn)題(8.3.1), (8.3.2)介紹的有限元解法，但其解題步驟對(duì)一般的微分方程定解問(wèn)題也是適用的. 對(duì)所給微分方程定解問(wèn)題，首先找出相應(yīng)微分方程的變分形式，然后進(jìn)行下列步驟：第1步 對(duì)定義區(qū)域(或定義區(qū)間)進(jìn)行網(wǎng)格剖分, 將定義區(qū)域(或定義區(qū)間)剖分為若干個(gè)小單元(一維是區(qū)間；多維是區(qū)域，如矩形、三角形等)；第2步 構(gòu)造試探空間; 即選擇適當(dāng)?shù)牟逯岛瘮?shù)類(lèi).第3步

20、 建立有限元方程組. 計(jì)算單元?jiǎng)偠染仃嚰坝叶讼蛄浚傩纬煽倓偠染仃嚰翱傆叶隧?xiàng)，寫(xiě)出有限元方程組；結(jié)合定解條件，求解有限元方程組.注 從形式上看，用有限元法解微分方程定解問(wèn)題很繁，但是有限元法的求解步驟規(guī)范，便于在計(jì)算機(jī)上實(shí)現(xiàn)，并且總剛度矩陣是三對(duì)角對(duì)稱(chēng)正定矩陣，因此有限元方程組可用第2章介紹的解三對(duì)角方程組的追趕法求解. 有限元法最主要的優(yōu)點(diǎn)是可以處理相當(dāng)復(fù)雜的區(qū)域上的邊值問(wèn)題。8.4 打靶法解常微分方程邊值問(wèn)題的打靶法(shooting method)，也稱(chēng)為嘗試法，其基本思想是把邊值問(wèn)題轉(zhuǎn)化為初值問(wèn)題來(lái)解：首先作出一些只滿(mǎn)足一端邊值條件的解，然后再?gòu)倪@些解中找出適合另一端邊值條件的解. 下

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

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

23、 初值問(wèn)題(8.4.5)的解是一條通過(guò)點(diǎn)、斜率為的曲線(xià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 用差分法解二階線(xiàn)性常微分方程的邊值問(wèn)題%用差分法解二階線(xiàn)性常微分方程的邊值問(wèn)題%a,b是區(qū)間, Step是步長(zhǎng), 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 取，用差分法求下列邊值問(wèn)題的近似解解 在MATLAB命令窗口輸入：f = inline(-4.*x); q= inline(4); DiffMethod(f, q, 0, 1, 0.1, 0, 2)回車(chē)，可的結(jié)果。8.5.2 用差分法解一般二階常微分方程的邊值問(wèn)題%用差分法解一般二階常微分方程的邊值問(wèn)題% a,b是區(qū)間, Step是步

26、長(zhǎng), 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 取，用差分法求下列邊值問(wèn)題的近似解解 在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)回車(chē)，可得結(jié)果。8.5.3 用有限元法

29、解二階常微分方程的邊值問(wèn)題%用有限元法解二階常微分方程的邊值問(wèn)題%a,b是區(qū)間, h是步長(zhǎng), 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