常微分方程邊值問題的數(shù)值解法

1、常微分方程論文題目： 常微分方程邊值問題的數(shù)值解法 組 長(zhǎng)： 數(shù)學(xué)132文洲 組 員： 數(shù)學(xué)131王琦 數(shù)學(xué)132姚瑤 信息132郭斌 院 （系）： 理學(xué)院 指導(dǎo)教師： 岳宗敏 時(shí) 間： 2015年6月9日 常微分方程邊值問題的數(shù)值解法摘要：作為一類定解問題，補(bǔ)充條件由以自變量取某些值時(shí)，未知函數(shù)及其導(dǎo)數(shù)的值而定，稱其為邊值條件。許多物理和數(shù)學(xué)問題都?xì)w結(jié)為邊值問題。本文介紹邊值問題的待定常數(shù)法和格林函數(shù)。關(guān)鍵詞：邊值問題 待定常數(shù)法 格林函數(shù)Abstract: as a kind of definite solution problems, the supplementary conditio

2、ns by took the certain values in the independent variable, the value of the unknown function and its derivative, referred to as boundary value conditions. Many physical and mathematical problems boil down to boundary value problems. In this paper, the boundary value problem of the method of undeterm

3、ined constants and greens function.Keywords: boundary value problem Method of undetermined constants Greens function11.1 引言在很多實(shí)際問題中都會(huì)遇到求解常微分方程邊值問題. 考慮如下形式的二階常微分方程, , (11.1.1) 在如下三種邊界條件下的定解問題:第一種邊界條件: , (11.1.2) 第二種邊界條件: , (11.1.2)第三種邊界條件: , (11.1.13)其中. 常微分方程邊值問題有很多不同解法, 這里只介紹打靶方法和有限差分方法.11.2 打靶算法 將

4、邊值問題轉(zhuǎn)化成初值問題來求解，即根據(jù)邊界條件(11.1.2)，也就是說，反復(fù)是調(diào)整初始時(shí)刻的斜率值，使得初值問題的積分曲線能“命中”。例11.1 利用打靶法求解兩點(diǎn)邊值問題解：(1) 為應(yīng)用打靶方法，需要假定初值來先求解初值問題，取初始參數(shù) ， 解 令，則上述初值問題變?yōu)椋?由標(biāo)準(zhǔn)的R-K 方法，取，可得。（2）再令， 解， 取。仍由標(biāo)準(zhǔn)的R-K 方法，可得。（3）令再解， 取，得。（4）令，得（5）類似有，因?yàn)椋?可以認(rèn)為即為所求。11.3 有限差分方法原理：微分方程數(shù)值解的方法?；舅枷胧前堰B續(xù)的定解區(qū)域用有限個(gè)離散點(diǎn)構(gòu)成的網(wǎng)格來代替， 這些離散點(diǎn)稱作網(wǎng)格的節(jié)點(diǎn)；把連續(xù)定解區(qū)域上的連續(xù)變量

5、的函數(shù)用在網(wǎng)格上定義的離散變量函數(shù)來近似；把原方程和定解條件中的微商用差商來近似， 積分用積分和來近似，于是原微分方程和定解條件就近似地代之以代數(shù)方程組，即有限差分方程組 ， 解此方程組就可以得到原問題在離散點(diǎn)上的近似解。然后再利用插值方法便可以從離散解得到定解問題在整個(gè)區(qū)域上的近似解。有限差分方法常用來求解常微分方程邊值問題, 下面就線性微分方程為例進(jìn)行討論. 設(shè)微分方程具有形式:, , (11.3.1) , . (11.3.2)不失一般性, 將求解區(qū)間等份, 分點(diǎn)為, , , 步長(zhǎng), 這里 稱為結(jié)點(diǎn).在內(nèi)部結(jié)點(diǎn)上, 記 , , , , 并分別用一階和二階中心差分來代替一階導(dǎo)數(shù)和二階導(dǎo)數(shù),

6、,.將上面各式代入方程(11.3.1), 略去項(xiàng), 得到關(guān)系式:,將上式改寫成, , (11.3.3)在上式中, , , (11.3.4) 利用邊界條件 , 可以得到如下形式的方程組: (11.3.5)其中 (11.3.6) (11.3.7) 系數(shù)矩陣(11.3.6)為n-1階三對(duì)角陣, 當(dāng)系數(shù)矩陣非奇異時(shí),方程組(11.3.5) 有唯一解,可以利用追趕法求解.對(duì)于第二種邊界條件: , (11.3.8)和第三種邊界條件: , (11.3.9)其中. 由于邊界條件中包含了導(dǎo)數(shù), 故邊界條件也必須用差商近似來表示. 因?yàn)椴荒芾脜^(qū)間以外的點(diǎn),所以在引進(jìn)(即 或 (即 ) 的近似式時(shí),就不能利用中心

7、差商. 如果只要求誤差是, 則可以用最簡(jiǎn)單的近似公式, (11.3.10)要使誤差達(dá)到, 可以利用Newton等距插值公式,得到如下近似公式:, (11.3.11), (11.3.12)將這些近似式代入邊值條件中, 再和方程(11.3.4)聯(lián)立, 就可以得到對(duì)應(yīng)的差分方程組. 例11.2 利用差分法求解兩點(diǎn)邊值問題解：將方程離散后寫成矩陣形式取， 可以得到三對(duì)角方程，下表11-1給出了計(jì)算結(jié)果。表 11-1取時(shí)計(jì)算結(jié)果，近似值近似值精確值0.11074610744107430.21170111696116950.31287712871128690.41429414286142840.51597

8、615968159650.61795817950179470.72028520276202740.82301323006230040.926219262152621411.4 應(yīng)用實(shí)例與Matlab11.4.1 MATLAB 關(guān)于常微分方程邊值問題數(shù)值解法的命令常微分方程初、邊值問題的符號(hào)解法 命令形式1：dsolve(，eqution，) 功能：求常微分方程eqution的解。 命令形式2：dsolve(，eqution，condl，cond2，var，) 功能：求常微分方程eqution的滿足初始條件的特解。其中eqution是求解的微分方程或微分方程組，condl，cond2是初始條件或

9、邊值，var是自變量.例11.3求解兩點(diǎn)邊值問題：.解 Matlab命令為 y=dsolve(xD2y-3Dy=x2，y(1)=0，y(5)=0,x)執(zhí)行結(jié)果： y -13x3+125468+31468x4 11.4.2 應(yīng)用實(shí)例當(dāng)流體以均勻的流速縱掠一平壁時(shí)，由于流體粘性的影響，在靠近壁面鄰近會(huì)形成很大的速度梯度，這就是速度邊界層用來描述速度邊界層的動(dòng)量方程經(jīng)過相似變換可以化為如下常微分方程兩點(diǎn)邊界值問題(著名的布拉修斯邊界層方程) (11.3.13)邊界條件為 (11.3.14)這里表示壁面噴注(或吸入). 根據(jù)給定的噴注速度值，(常數(shù))表示吸入， (常數(shù))表示噴注.為了利用打靶法求解該問

10、題的解，假設(shè)（這里是一個(gè)待定常數(shù)，需要在計(jì)算中估計(jì)). 表11-1 給出了數(shù)值計(jì)算方法。當(dāng)時(shí)取、和.這里是一個(gè)常數(shù)(估計(jì)值)，這樣就能根據(jù)式 (11.3.13) - (11.3.14), 求得.接著，令具有某一增量，分別計(jì)算、和，繼續(xù)進(jìn)行下去一直到邊界外邊緣(在平壁的情況下它相當(dāng)于). 此時(shí)，應(yīng)該得到. 如果這個(gè)要求不能滿足，就要修改原來的估計(jì)值，直至使這個(gè)邊界條件得到滿足為止.，表11-2 給出了數(shù)值計(jì)算方法。表11-2 布拉修斯方程式的數(shù)值計(jì)算方法0N0式(1.31)0.1N+0.1*00+0.1*+0.1*式(1.31)0.2N+0.1*0+0.1*(0+0.1*)0+0.1*+0.1*

11、+0.1*)+0.1*+0.1*式(1.31)50.991.0在程序中為了能夠方便的變換系數(shù)，我們把方程中系數(shù)設(shè)為epsilo，從而把方程改寫為： (11.3.13) function itBlasius(n,m,minbeta,maxbeta,dbeta,dita,maxita,tol)found=0;fprintf(n白拉修斯方程數(shù)值計(jì)算迭代程序nn正在計(jì)算.nn);for beta=minbeta:dbeta:maxbeta ita=Blasius(n,m,beta,dita,maxita,tol,1,1,1,1,0); if ita-1 fprintf(n滿足白拉修斯收斂條件的beta

12、已經(jīng)找到，beta=%fnn,beta); found=found+1; break endendif found=0 fprintf(n在給定的范圍 (%f.%f) 內(nèi)無法找到滿足白拉修斯收斂條件的beta.nn,minbeta,maxbeta);endMatlab 命令：itBlasius(0,0,0,0.5,0.00001,0.1,5,0.001)beta=0.得出每一步的迭代結(jié)果：Matlab 命令：blasius(0,0,0.,0.1,5,0.001,1,1,1,1,1)beta=0.表11-2給出了當(dāng)情況在區(qū)間0，4的部分計(jì)算結(jié)果.Matlab程序function itaconv=

13、Blasius(n,m,beta,dita,maxita,tol,options)format long;epsilo=0.5;if options(5)=1 fprintf(nbeta=%fnttitattfttfttfttfnn,beta);endfound=0;x0=n;y0=m;z0=beta;w0=-x0*z0*epsilo;for ita=0 : dita : maxita switch options(1) case 1 multiplier=0.1; otherwise fprintf(n出錯(cuò)了！nn); break end switch options(2) case 1 x

14、=x0+multiplier*y0; otherwise fprintf(n出錯(cuò)了！nn); break end switch options(3) case 1 y=y0+multiplier*z0;%推薦 otherwise fprintf(n出錯(cuò)了！nn); break end switch options(4) case 1 z=z0+multiplier*w0;%推薦 otherwise fprintf(n出錯(cuò)了！nn); break end w=-x*z*epsilo; if options(5)=1 fprintf(t%2.5ft%2.5ft%2.5ft%2.5ft%2.5fn,

15、ita,x0,y0,z0,w0); end if abs(1-y)tol found=1;itaconv=ita; break end x0=x;y0=y;z0=z;w0=w;endif found itaconv=-1;end下面表11-3 給出了當(dāng)情況在區(qū)間0，5的部分計(jì)算結(jié)果表11-3：當(dāng)部分計(jì)算結(jié)果. 0000.332060.20.006640.066410.331990.40.026560.132770.331470.60.059740.198940.330080.80.106110.264710.3273910.165570.329790.323011.20.237950.395780.316591.40.322980.456270.307871.60.420320.516760.296671.80.529520.574770.2829320.650030.629770.266752.20.78120.681320.248352.40.92230.728990.228092.61.072520.7724