中心差分解兩點(diǎn)邊值問(wèn)題

1、一題目用中心差分格式計(jì)算如下兩點(diǎn)邊值問(wèn)題已知其精確解為二理論作為模型，考慮兩點(diǎn)邊值問(wèn)題：（1.1）（1.2）假定是給定的常數(shù)。1. 建立差分格式（1）.區(qū)域網(wǎng)格剖分首先取個(gè)節(jié)點(diǎn)：將區(qū)間分成個(gè)小區(qū)間：于是得到區(qū)間的一個(gè)網(wǎng)格剖分。記，稱為網(wǎng)格最大步長(zhǎng)。用表示網(wǎng)格內(nèi)點(diǎn)，的集合，表示內(nèi)點(diǎn)和界點(diǎn)的集合。取相鄰節(jié)點(diǎn)的中點(diǎn)，稱為半整數(shù)點(diǎn)。則由節(jié)點(diǎn)又構(gòu)成的一個(gè)網(wǎng)格剖分，稱為對(duì)偶剖分。（2）.微分方程的離散，建立相應(yīng)差分格式用差商代替微商，將方程（1.1）在內(nèi)點(diǎn)離散化.注意對(duì)充分光滑的，由Taylor展式有（1.3） （1.5）由（1.5）減（1.4），并除以，得（1.6）令則由（1.3）（1.6）知，邊值問(wèn)題

2、的解滿足方程： （1.7） 其中 （1.8）為差分算子的截?cái)嗾`差，舍去，便得逼近邊值問(wèn)題（1.1）（1.2）的差分方程： （1.9）i=1,2,N-1，由方程（1.7）（1.9），截?cái)嗾`差可表示為 （1.10）當(dāng)網(wǎng)格均勻，即時(shí)差分方程（1.9）簡(jiǎn)化為 （1.11）這相當(dāng)于用一階中心差商，二階中心差商依次代替（1.1）的一階微商和二階微商的結(jié)果。這個(gè)方程就是中心差分格式。截?cái)嗾`差為： （1.12）所以截?cái)嗾`差按或的階為。在本題中， ，因?yàn)閞=0方程（1.11）的系數(shù)對(duì)角矩陣是三對(duì)角矩陣。我們可以用消元法或迭代法求解方程組（1.1）（1.2）式（1.11）用方程組展開(kāi)：寫(xiě)成矩陣形式為：2.收斂性分

3、析根據(jù)（1.10）我們引進(jìn)誤差則誤差函數(shù)滿足下列差分方程：于是收斂性及收斂速度的估計(jì)問(wèn)題，就歸結(jié)到通過(guò)右端（截?cái)嗾`差）估計(jì)誤差函數(shù)的問(wèn)題。由（1.12）我們知，有從而差分方程滿足相容條件。若引進(jìn)記號(hào)，設(shè)則可將（1.9）改寫(xiě)為將差分解表成 （2.1）其中滿足 （2.2）而滿足 （2.3）先估計(jì)，由 （2.4）據(jù)差分格林公式再利用柯西不等式，有常數(shù)使 （2.5）將不等式（2.6）用于（2.5）右端，則 （2.6）解差分方程（2.2，易得）從而這樣， （2.7）利用范數(shù)，從（2.7）推出 （2.8）因?yàn)橐虼?（2.9）聯(lián)結(jié)（2.1）（2.7）及（2.9）即得差分解的先驗(yàn)估計(jì)： （2.10）其中不等式

4、（2.10）說(shuō)明差分解連續(xù)依賴于右端和邊值，因此差分格式（1.11）關(guān)于右端及邊值穩(wěn)定.根據(jù)定理1.1 ： 若邊值問(wèn)題的解u充分光滑，差分方程按滿足相容條件且關(guān)于右端穩(wěn)定，則差分解按收斂到邊值問(wèn)題的解，且有和相同的收斂階。所以差分方程的解的收斂速度為。三程序代碼：clcclfclfsyms x;a=1; %區(qū)間界點(diǎn)b=2; %區(qū)間界點(diǎn)p=exp(x); %這是p函數(shù)q=sin(x)+1+x; %這是q函數(shù)f=-exp(x)*(2*x+1)+(sin(x)+1+x)*x*(x-1);%這是f函數(shù)r=0; %這是r函數(shù).N=10; %將區(qū)間劃分的等分,這里控制！h=(b-a)/N; %這里確定步長(zhǎng)

5、value_of_f=zeros(N-1,1);%這是fdiag_0=zeros(N-1,1);%確定A的對(duì)角元diag_1=zeros(N-2,1);%確定A的偏離對(duì)角的上對(duì)角元diag_2=zeros(N-2,1);%確定A的偏離對(duì)角的下對(duì)角元X=a:h:b;u_a=0; %邊界條件u_b=2; %邊界條件for j=2:N diag_0(j-1)=(subs(p,x,(X(j+1)+X(j)/2)+(subs(p,x,(X(j-1)+X(j)/2)/(h2)+(subs(q,x,X(j);end %獲取對(duì)角元素for j=3:N diag_2(j-2)=-(subs(p,x,(X(j-1

6、)+X(j)/2)/(h2)-subs(r,x,X(j)/(2*h);end %獲取A的第三條對(duì)角for j=2:N-1 diag_1(j-1)=-(subs(p,x,(X(j+1)+X(j)/2)/(h2)+subs(r,x,X(j)/(2*h);end %獲取A的第二條對(duì)角for j=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(2)+X(1)/2)/(h2);value_of_f(N-1)=value_of_f(N-1)+u_b*(subs(p,x,(X(N

7、)+X(N+1)/2)/(h2);A=diag(diag_0)+diag(diag_1,1)+diag(diag_2,-1);%組裝系數(shù)矩陣format longU=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-precise_value' ; %誤差deta_max=max(abs(deta);%最大誤差fprintf(

8、'最大的誤差是%fn',deta_max)plot(X(2:N),U,'b*',X(2:N),precise_value,'r-') %差分解與精確解對(duì)比表figure();plot(X(2:N),deta) %誤差圖結(jié)果： X的值步長(zhǎng)h2.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.440001.710000以下僅給出步長(zhǎng)為N=20,h=0.05的精確值和差分值圖與誤差圖，其