中心差分解兩點邊值問題

-.z一．題目用中心差分格式計算如下兩點邊值問題其準確解為二．理論作為模型，考慮兩點邊值問題：…………〔1.1〕…………〔1.2〕假定是給定的常數。建立差分格式〔1〕.區(qū)域網格剖分首先取個節(jié)點：將區(qū)間分成個小區(qū)間：于是得到區(qū)間的一個網格剖分。記，稱為網格最大步長。用表示網格點，，，的集合，表示點和界點的集合。取相鄰節(jié)點的中點，稱為半整數點。則由節(jié)點又構成的一個網格剖分，稱為對偶剖分。〔2〕.微分方程的離散，建立相應差分格式用差商代替微商，將方程〔1.1〕在點離散化.注意對充分光滑的，由Taylor展式有………〔1.3〕………〔1.5〕由〔1.5〕減〔1.4〕，并除以，得…………〔1.6〕令則由〔1.3〕〔1.6〕知，邊值問題的解滿足方程：…………〔1.7〕其中…………〔1.8〕為差分算子的截斷誤差，舍去，便得逼近邊值問題〔1.1〕〔1.2〕的差分方程：…………〔1.9〕i=1,2,…,N-1，由方程〔1.7〕〔1.9〕，截斷誤差可表示為…………〔1.10〕當網格均勻，即時差分方程〔1.9〕簡化為…………〔1.11〕這相當于用一階中心差商，二階中心差商依次代替〔1.1〕的一階微商和二階微商的結果。這個方程就是中心差分格式。截斷誤差為：…………〔1.12〕所以截斷誤差按或的階為。在此題中，，，，，因為r=0方程〔1.11〕的系數對角矩陣是三對角矩陣。我們可以用消元法或迭代法求解方程組〔1.1〕〔1.2〕式〔1.11〕用方程組展開：寫成矩陣形式為：2.收斂性分析根據〔1.10〕我們引進誤差則誤差函數滿足以下差分方程：于是收斂性及收斂速度的估計問題，就歸結到通過右端〔截斷誤差〕估計誤差函數的問題。由〔1.12〕我們知，有從而差分方程滿足相容條件。假設引進記號，，，，，設則可將〔1.9〕改寫為將差分解表成…………〔2.1〕其中滿足…………〔2.2〕而滿足…………〔2.3〕先估計，由…………〔2.4〕據差分格林公式再利用柯西不等式，有常數使…………〔2.5〕將不等式〔2.6〕用于〔2.5〕右端，則…………〔2.6〕解差分方程〔2.2，易得〕從而這樣，…………〔2.7〕利用數，從〔2.7〕推出…………〔2.8〕因為因此…………〔2.9〕聯結〔2.1〕〔2.7〕及〔2.9〕即得差分解的先驗估計：…………〔2.10〕其中不等式〔2.10〕說明差分解連續(xù)依賴于右端和邊值，因此差分格式〔1.11〕關于右端及邊值穩(wěn)定.根據定理1.1：假設邊值問題的解u充分光滑，差分方程按滿足相容條件且關于右端穩(wěn)定，則差分解按收斂到邊值問題的解，且有和一樣的收斂階。所以差分方程的解的收斂速度為。三．程序代碼：clcclfclfsyms*;a=1;%區(qū)間界點b=2;%區(qū)間界點p=e*p(*);%這是p函數q=sin(*)+1+*;%這是q函數f=-e*p(*)*(2**+1)+(sin(*)+1+*)***(*-1);%這是f函數r=0;%這是r函數.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的偏離對角的下對角元*=a:h:b;u_a=0;%邊界條件u_b=2;%邊界條件forj=2:Ndiag_0(j-1)=((subs(p,{*},{(*(j+1)+*(j))/2}))+(subs(p,{*},{(*(j-1)+*(j))/2})))/(h^2)+(subs(q,{*},{*(j)}));end%獲取對角元素forj=3:Ndiag_2(j-2)=-((subs(p,{*},{(*(j-1)+*(j))/2})))/(h^2)-subs(r,{*},{*(j)})/(2*h);end%獲取A的第三條對角forj=2:N-1diag_1(j-1)=-((subs(p,{*},{(*(j+1)+*(j))/2})))/(h^2)+subs(r,{*},{*(j)})/(2*h);end%獲取A的第二條對角forj=2:N;value_of_f(j-1)=subs(f,{*},{*(j)});end%獲取F值value_of_f(1)=value_of_f(1)+u_a*(subs(p,{*},{(*(2)+*(1))/2}))/(h^2);value_of_f(N-1)=value_of_f(N-1)+u_b*(subs(p,{*},{(*(N)+*(N+1))/2}))/(h^2);A=diag(diag_0)+diag(diag_1,1)+diag(diag_2,-1);%組裝系數矩陣formatlongU=inv(A)*value_of_f%差分解%fprintf('%11.5f',U)fprintf('\n');d*=*(2:N);precise_value=d*.*(d*-1)%準確解%fprintf('%11.5f',precise_value)deta=U-precise_value';%誤差deta_ma*=ma*(abs(deta));%最大誤差fprintf('最大的誤差是%f\n',deta_ma*)plot(*(2:N),U,'b*',*(2:N),precise_value,'r--')%差分解與準確解比照表figure();plot(*(2:N),deta)%誤差圖結果：*的值步長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以下僅給出步長為N=20,h=0.