偏微分方程數(shù)值解法

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

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

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

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

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

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

7、逼近精確解。六、附錄取，則，滿足穩(wěn)定性條件取，則，滿足穩(wěn)定性條件另取，則，亦滿足穩(wěn)定性條件程序代碼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 %求解數(shù)值解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) %求誤差的最大值精確解與數(shù)值解的比較：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時精確解和數(shù)值解的比較')text(0.05,21,'藍：精確解');text(0.05,20,'黃：數(shù)值解');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