1、偏微分方程數值解法課 程 設 計題 目： 六點對稱差分格式解熱傳導方程的初邊值問題姓 名： 王曉霜 學 院： 理學院 專 業： 信息與計算科學 班 級： 0911012 學 號： 091101218 指導老師：翟方曼2012年12月14日一、題目用六點對稱差分格式計算如下熱傳導方程的初邊值問題已知其精確解為二、理論1考慮的問題考慮一維模型熱傳導方程(1.1) ，其中為常數。是給定的連續函數。（1.1）的定解問題分兩類： 第一，初值問題（Cauchy 問題）：求足夠光滑的函數，滿足方程（1.1）和初始條件：（1.2） , 第二，初邊值問題（也稱混合問題）：求足夠光滑的函數，滿足方程（1.1）和初

2、始條件： , 及邊值條件 ， 假定和在相應的區域光滑，并且于，兩點滿足相容條件,則上述問題有唯一的充分光滑的解。現在考慮邊值問題（1.1），（1.3）的差分逼近取 為空間步長，為時間步長，其中，是自然數，, ； , 將矩形域分割成矩形網格。其中 表示網格節點；表示網格內點（位于開矩形中的網格節點）的集合；表示位于閉矩形中的網格節點的集合；表示-網格邊界點的集合。表示定義在網點處的待求近似解，。注意到在節點處的微商和差商之間的下列關系（）：2區域網格剖分取空間步長和時間步長，其中都是正整數。用兩族平行直線和將矩形域分割成矩形網格，網格節點為。以表示網格內點集合，即位于開矩形的網點集合；表示所有位

3、于閉矩形的網點集合；是網格界點集合。其次，用表示定義在網點的函數，3建立相應差分格式數值分析中，Crank-Nicolson方法是有限差分方法中的一種，用于數值求解熱方程以及形式類似的偏微分方程。它在時間方向上是隱式的二階方法，數值穩定。該方法誕生于20世紀，由John Crank與Phyllis Nicolson發展。向前差分格式 ， =0 向后差分格式 ， =0將向前差分格式和向后差分格式做算術平均，得到的差分格式稱之為六點對稱格式，也稱為Grank-Nicholson格式： ， =0進一步， +=+按層計算：首先，取，則利用初值和邊值=0，來確定出第一層的，即求解方程組：+=+，=0。求

4、出，在由,取，可利用，解出，。如此下去，即可逐層算出所有，。若記三、截斷誤差=。注意：=又兩式相加而+故有。四、穩定性與收斂性 拋物方程的兩層差分格式可以統一寫成向量形式:(1.7) 其中，和是階矩陣。我們假定可逆，即（1.7）是唯一可解的。對于顯格式，等于單位矩陣。三層格式可以通過引入新變量化成兩層格式。假設差分解的初始值（其實可以是任一層的值）有誤差，以后各層計算沒有誤差，讓我們來考察初始誤差對以后各層的影響。令和分別是以和為初始值由差分格式（1.7）得到的兩組差分解，則滿足（1.8） 因此，按初值穩定應該意味著。這就導致如下定義： 假設，我們稱差分格式（2.1）按初值穩定，如果存在正常數

5、和，使得以下不等式成立：（1.8） ， 這里是上的某一個范數，例如 類似地，假設，我們稱差分格式（2.1）按右端穩定，如果存在正常數和，使得以下不等式成立：（1.8） ， 可以證明，差分格式若按初值穩定，則一定按右端穩定。因此，這時我們簡單地稱差分格式穩定。前面討論的向前差分格式當網比時穩定，當時不穩定。這就意味著給定空間步長以后，時間步長必須足夠小，才能保證穩定。而向后差分格式和Grank-Nicholson格式（1.6）則對任何網比都是穩定的，時間步長可以取得大一些，從而提高運算效率。如果某個差分格式的截斷誤差當和趨于0時隨之趨于0，則稱這個差分格式是相容的。可以證明：若差分格式是相容的和

6、穩定的，則它是收斂的，并且差分解與微分解之間誤差的階等于截斷誤差的階。因此，對任何網比，向后差分格式（1.6）有收斂階。五、結論對于擴散方程（包括許多其他方程），可以證明Crank-Nicolson方法無條件穩定。但是，如果時間步長與空間步長平方的比值過大（一般地，大于1/2），近似解中將存在虛假的振蕩或衰減。基于這個原因，當要求大時間步或高空間分辨率的時候，往往會采用數值精確較差的后向歐拉方法進行計算，這樣即可以保證穩定，又避免了解的偽振蕩。由本題可以總結出，拋物型方程的六點對稱差分法所得的數值解能夠較好地逼近方程的精確 解，且區域剖分得越細，即步長越小，數值解與精確解的誤差就越小，數值解越

7、逼近精確解。六、附錄取，則，滿足穩定性條件取，則，滿足穩定性條件另取，則，亦滿足穩定性條件程序代碼format longa=2;l=1;T=1;N=5;M=100;h=l/N;to=T/M;r=(a*to)/h2;for j=1:N+1 x(j)=(j-1)*h; for k=1:M+1 t(k)=(k-1)*to; u(j,k)=exp(x(j)+2*t(k); endendu %求解精確解for j=1:N+1 x(j)=(j-1)*h; us(j,1)=exp(x(j);endfor k=1:M+1 t(k)=(k-1)*to; us(1,k)=exp(2*t(k); us(N+1,k)

8、=exp(1+2*t(k);endfor k=2:M+1 for j=2:N us(j,k)=r*us(j-1,k-1)+(1-2*r)*u(j,k-1)+r*us(j+1,k-1); endendus %求解數值解for k=1:M+1 for j=1:N+1 R(j,k)=abs(u(j,k)-us(j,k); endendR %計算誤差Rmax=max(max(R) %求誤差的最大值精確解與數值解的比較：x=0:0.1:1;hold onplot(x,u(:,M+1),'b');plot(x,us(:,M+1),'y');title('t=1，h=1/10，=1/400時精確解和數值解的比較')text(0.05,21,'藍：精確解');text(0.05,20,'黃：數值解');hold off取不同步長時的誤差比較：x=0:1/5:1;y=0:1/10:1;z=0:1/20:1;hold onplot(x,R(:,M+1),'b');hold o