1、題目用中心差分格式計算如下兩點邊值問題dxduxesinx1xuxe2x1sinx1xxx1,1x2dxdxu(1)0,u(2)2已知其精確解為u(x)x(x1)理論作為模型，考慮兩點邊值問題：,ddu、duLu(p(x)-)r丁q(x)uf(x),axb(1.1)dxdxdxu(a),u(b)(1.2)假定pC1a,b,p(x)Pmin0,r,q,fCa,b,是給定的常數(shù)。建立差分格式(1).區(qū)域網(wǎng)格剖分首先取N1個節(jié)點：ax0為LxiLxnb,將區(qū)間Ia,b分成N個小區(qū)間：Ii:為1xxi,i1,2,LN于是得到區(qū)間I的一個網(wǎng)格剖分。記hxixi1,稱hmaxh為網(wǎng)格最大步長。用Ih表i示

2、網(wǎng)格內(nèi)點x1,x2,L,xN1的集合，Ih表小內(nèi)點和界點冷a,xNb的集合。1取相鄰下點xi1,xi的中點x1-(xi1xi)(i1,2,L,N),稱為半整數(shù)點。則由節(jié)點i22a冷xx3Lx1Lx1-i-N-2222xnb又構(gòu)成a,b的一個網(wǎng)格剖分，稱為對偶剖分。(2).微分方程的離散，建立相應(yīng)差分格式用差商代替微商，將方程(1.1)在內(nèi)點xi離散化.注意對充分光滑的u，由Taylor展式有(1.3)u(xi1)u(xii)duhihi1頃1hi1hid2u2、丁"O(h)p(x1)i2rduP.dxiu(xi)u(xi1)du_P.1dxi芝hi24iO(h3)P(x1)i2u(x

3、i1)u(xi)dupdxihi1由(1.5)減(1.4),并除以2=P(xi1)hihi1i22dun-(phihi1dxid,duhi(P)i-dxdxhihi1得一2一，苛u(x1)u(x)hi1rdu_Pdxihi.d223hirduP乩1O(h)24dx|芝(1.4)2n1dup3iO(h)24dx(1.5)p(x)件*i2hiP四iO(h2)12dx3ihid3u2、-p3iO(h)dx1)2duh1Ht(P)iFdxdx12(1.6)令P.1i-2滿足方程:p(xii),rir(xi),qi2q(x),fif(x),則由(1.3)(1.6)知，邊值問題的解u(x)LhU(Xi)2

4、rp1hih1i2u(x1)u(xi)hi1u(x)12u(xi1)hi(1.7)其中ri-一-一u(xi1)hihi1U(Xi1)qiU(Xi)fiR(u)Ri(u)(hi1d2dui(P)i1hi)(一4dxdx1d3uP3i12dx1d2u2】")O(h)(1.8)為差分算子Lh的截斷誤差，舍去R(u),便得逼近邊值問題(1.1)(1.2)的差分方程:、u(xi)土p1心1)心)hihi1i2hi1Piu(xi)u(xi)xhi(1.9)ru(xi1)u(xi1)qiu(xi)hih1fii=1,2，N-1,u°,un由方程(1.7)(1.9)，截斷誤差R(u)可表示

5、為(1.10)R(u)LhU(x)LhUiLh(u(Xi)u)當網(wǎng)格均勻，即hh(i1,2,L,N)時差分方程(1.9)簡化為LhUi1r，、-rp1Ui1(P1P1)Uiiiih222P.1山1i2Ui1Ui1r-qufi2h(1.11)1.1)的一階微商和二階微商的結(jié)果。這相當于用一階中心差商，二階中心差商依次代替(這個方程就是中心差分格式。截斷誤差為：R(u)1d2du1d3uh(一2(p)ip3i4dx2dx12dx3<罰i)O(h2)2dx(1.12)xe2x1sinx1xxx1所以截斷誤差按|Rh(u)|。或|Rh(u)|C的階為O(h2)。在本題中，p(x)ex,r(x)0

6、,q(x)sinx1x,f(x)0,2因為r=0方程(1.11)的系數(shù)對角矩陣是三對角矩陣。我們可以用消元法或迭代法求解方程組(1.1)(1.2)式(1.11)用方程組展開：1,、1rp3U22(P3p1)q1U12P1U0f1h2h22h2LL1hpk1Uk1pk)2qkUph2N1UN22(p1hN1p3)qN1UN1N一2pi)qmf1h2pi12ph2N1,2(p1p3)qN1UN1NN-寫成矩陣形式為:Pi)qh221h2P21-1祁p312(p5p3)。h2201Rp；LLLL0000uU2LLLM0LL01h2pN；+(pN3PN5)qN2221h2PN|UN2UN10LL001

7、1,C"ch2(PN1Pn3)*1h2Nf12Plh2MfN2fN12p1hN12收斂性分析根據(jù)(1.10)我們引進誤差eu(x)u則誤差函數(shù)eh(Xi)e滿足下列差分方程:啟R(u)i1,2,L,N1Ri(u)(截斷誤差)估計誤差函數(shù)eeoeN0的問題。由(1.12)我們知，有晚IR(u)|0從而差分方程滿足相容條件。若引進記號/、ViVi1(Vi)V,(Vi)xVi,xxi,xhVi1Viu，(V)Vhi1xi,xVi1Vhihi2(hihi1),h0,hN于是收斂性及收斂速度的估計問題，就歸結(jié)到通過右端設(shè)PminCo則可將(1.9)改寫為LhU(pi(Ui)qufiIxx2將差

8、分解Ui表成(2.1)UiUiUi,iIh(Ui)xx0,iIh,U0Uo,UNUn(2.2)而Ui滿足LhUifiLhUi,iIh,U。UNo(2.3)先估計Uh,由Co|(Uh)|x0(fhLhUh,Uh)ih(fh,Uh)Ih(LhUhUh).其中Ui滿足據(jù)差分格林公式2.4)(LhUh,Uh)ih(Ph(Uh),(Uh)ih(qhUh,Uh)ihxx再利用柯西不等式，有常數(shù)&使|(LhUh,Uh)|h|Ci|(Uh)|0C2|Uh|0|(Uh)|(2.5)xx將不等式(2.6)用于(2.5)右端，則1,GC2|(Uh)|o|fh|i-|(Uh)|o-|Uh|o(2.6)xx解差

9、分方程（2.2,易得）從而UnUibUo/(為aa)Uo1|(Uh)|o.|UnUo|x.(ba)|Uh|°.(ba)max|Ui|,(ba)(|U°|un|)iIh這樣，|u5川廠bapduol|uNI)(2.7)利用范數(shù)|uhIk,從(2.7)推出ll(Uh)Ho-IIMhJ|Uh|k(2.8)xCoCo因為IIUh|oll(Uh)lb2x因此l|Uh卅|Uh|0|(Uh)|oXba(1)ll(Uh)|o(2.9)2xba1CiCo(1)-|fh|i-lluhlh2CoCo聯(lián)結(jié)(2.1)(2.7)及(2.9)即得差分解的先驗估計：II"|k|Uh|k|UhM(

10、20)MJIMHM2(|u0|unI)其中1ba、M1(1c°2M2J(ba)(ba)'1(12c0不等式(2.10)說明差分解連續(xù)依賴于右端和邊值，因此差分格式(1.11)關(guān)于右端及邊值穩(wěn)定.根據(jù)定理1.1:若邊值問題的解u充分光滑，差分方程按舊"滿足相容條件且關(guān)于右端穩(wěn)定，則差分解山按|劇收斂到邊值問題的解，且有和|Rh(u)|R相同的收斂階。所以差分方程的解的收斂速度為0(7)。三.程序代碼：clcclfclfsymsx;a=1;%區(qū)間界點b=2;%區(qū)間界點p=exp(x);%這是p函數(shù)q=sin(x)+1+x;%這是q函數(shù)f=-exp(x)*(2*x+1)+

11、(sin(x)+1+x)*x*(x-1);%這是f函數(shù)r=0;%這是r函數(shù).N=10;%將區(qū)間劃分的等分，這里控制！h=(b-a)/N;%這里確定步長value_of_f=zeros(N-1,1);%這是fdiag_0=zeros(N-1,1);%確定A的對角元diag_1=zeros(N-2,1);%確定A的偏離對角的上對角元diag_2=zeros(N-2,1);%確定A的偏離對角的下對角元X=a:h:b;u_a=0;%邊界條件u_b=2;%邊界條件forj=2:Ndiag_0(j-1)=(subs(p,x,(X(j+1)+X(j)/2)+(subs(p,x,(X(j-1)+X(j)/2)

12、/(hA2)+(subs(q,(x，(X(j);end%獲取對角元素forj=3:Ndiag_2(j-2)=-(subs(p,x,(X(j-1)+X(j)/2)/(hA2)-subs(r,x,X(j)/(2*h);end%獲取A的第三條對角forj=2:N-1diag_1(j-1)=-(subs(p,x,(X(j+1)+X(j)/2)/(hA2)+subs(r,x,X(j)/(2*h);end%獲取A的第二條對角forj=2:N;value_of_f(j-1)=subs(f,x,X(j);end%獲取F值value_of_f(1)=value_of_f(1)+u_a*(subs(p,x,(X(

13、2)+X(1)/2)/(hA2);value_of_f(N-1)=value_of_f(N-1)+u_b*(subs(p,x,(X(N)+X(N+1)/2)/(hA2);A=diag(diag_0)+diag(diag_1,1)+diag(diag_2,-1);%組裝系數(shù)矩陣formatlongU=inv(A)*value_of_f%差分解%fprintf('%11.5f,U)fprintf('n');dx=X(2:N);precise_value=dx.*(dx-1)%精確解%fprintf('%11.5f,precise_value)deta=U-preci

14、se_value'%誤差deta_max=max(abs(deta);%最大誤差fprintf('最大的誤差是%fn',deta_max)plot(X(2:N),U,'b*',X(2:N),precise_value,'r-')%差分解與精確解對比表figure();plot(X(2:N),deta)%誤差圖結(jié)果：、'<X.的值步長卜、2.12.22.32.42.52.62.72.82.9最大誤差0.10.110150.240260.390330.560370.750380.960381.190311.440221.710120.0003800.050.110030.240060.390080.560090.750090.960081.190071.440051.710030.0000950.0250.110000.240010.390020.560020.750020.960021.190021.440011.710010.0000240.01250.110000.240000.390000.560010.750010.960001.190001.440001.710000.000006精確解0.110000.240000.390000.560000.750000.960001.190001