拋物形擴散方程的有限差分法及數(shù)值實例

1、偏微分方程數(shù)值解 所在學(xué)院： 數(shù)學(xué)與統(tǒng)計學(xué)院 課題名稱：拋物形擴散方程的有限差分法及數(shù)值實例 學(xué)生姓名： 向聘 拋物形擴散方程的有限差分法及數(shù)值實例1.1拋物型擴散方程拋物型偏微分方程是一類重要的偏微分方程?？紤]一維熱傳導(dǎo)方程： (1.1.1)其中是常數(shù)，是給定的連續(xù)函數(shù)。按照初邊值條件的不同給法，可將(1.1.1)的定解分為兩類：第一，初值問題（Cauchy 問題）：求足夠光滑的函數(shù)，滿足方程(1.1.1)和初始條件： , (1.1.2)第二，初邊值問題（也稱混合問題）：求足夠光滑的函數(shù)，滿足方程(1.1.1)和初始條件： , (1.1.3)及邊值條件 ， (1.1.4)假定和在相應(yīng)的區(qū)域光

2、滑，并且于，兩點滿足相容條件,則上述問題有唯一的充分光滑的解。1.2拋物線擴散方程的求解下面考慮如下熱傳導(dǎo)方程 (1.2.1)其中，,(常數(shù))是擴散系數(shù)。取為空間步長，為時間步長，其中，是自然數(shù)，用兩族平行直線, 和, 將矩形域分割成矩形網(wǎng)格。其中 表示網(wǎng)格節(jié)點；表示網(wǎng)格內(nèi)點（位于開矩形中的網(wǎng)格節(jié)點）的集合；表示位于閉矩形中的網(wǎng)格節(jié)點的集合；表示-網(wǎng)格邊界點的集合。表示定義在網(wǎng)點處的待求近似解，?，F(xiàn)在對方程進行差分近似：(一) 向前差分格式 (1.2.2)， =0 (1.2.3)計算后得： (1.2.4)其中，。顯然，這是一個四點顯示格式，每一層各個節(jié)點上的值是通過一個方程組求解到的。方程組如

3、下： (1.2.5)若記，則顯格式(1.2.4)可寫成向量形式 (1.2.6)其中而對于向前差分格式，當(dāng)網(wǎng)比時穩(wěn)定，當(dāng)時不穩(wěn)定。這就意味著給定空間步長以后，時間步長必須足夠小，才能保證穩(wěn)定。1.3拋物型熱傳導(dǎo)方程數(shù)值算例對于(1.2.1)所描述的擴散方程，取已知方程的精確解為 ： (1.3.1)設(shè)空間步長，時間步長為，網(wǎng)格比。向前格式： 邊值條件：，.初值條件： 對時間和空間進行分割，令M=40，N=1600，通過Matlab計算得到該方程的解析解，數(shù)值解以及相對誤差如下：圖(1)解析解的圖像圖(2)數(shù)值解的圖像圖(3)M=40,N=1000的相對誤差的圖像我們?nèi)〔糠志_解和數(shù)值解進行比較，結(jié)

4、果如表(1)數(shù)值解精確解相對誤差0.10.10.11510.1152 0.20.20.08150.08160.30.30.04180.04190.40.40.01830.01840.50.450.01170.01180.60.350.03000.03010.70.250.06840.06860.80.150.50840.50850.90.050.36540.3654表(1)數(shù)值解與精確解的比較由表(1)我們可以看出，精確解和數(shù)值解的絕對誤差在以內(nèi)，因此可以得出，在分割M=40，N=1600下，該有限差分方法對方程(1.3.1)是收斂和穩(wěn)定的。下面，我們比較在不同的分割下對有限差分算法精度的影響

5、。在擴散系數(shù)不變的情況下，講時間和空間進行更加細(xì)密的分割，取，其中，M表示空間上的分割，N表示時間上的分割。觀察數(shù)值解與精確解在節(jié)點處的絕對誤差值，如下圖所示：圖(4)M=50,N=10000的相對誤差的圖像由圖(3)和圖(4)，兩者在節(jié)點處的誤差收斂分別是在和以內(nèi)，因此，可以得出的結(jié)論是：在收斂范圍內(nèi)，隨著時間和空間的分割越細(xì)，節(jié)點數(shù)越多，精確解和解析解之間的絕對誤差也越小，有限差分法的算法精度也越高。最后，我們比較網(wǎng)比以及時擴散方程的收斂情況。當(dāng)網(wǎng)比時，此時我們?nèi)=10，N=50，這時，方程的數(shù)值解與解析解還有相對誤差圖如下：圖(5)M=10,N=50的解析解的圖像圖(6)M=10,N=

6、50的數(shù)值解的圖像圖(7)M=10,N=50的絕對誤差的圖像此時，我們觀察絕對誤差發(fā)現(xiàn)，擴散方程(1.3.1)時不收斂不穩(wěn)定的。而前面我們已經(jīng)知道，到網(wǎng)格比為時，方程是收斂穩(wěn)定的。所以，我們可以驗證，當(dāng)網(wǎng)比時穩(wěn)定，當(dāng)時不穩(wěn)定。參考文獻1李榮華，劉播.微分方程數(shù)值解法M.北京.高等教育出版社.2009.1.2王曰朋.偏微分方程數(shù)值解OL. 3未知.偏微分方程的Matlab解法OL.4周品，何正風(fēng).MATLAB數(shù)值分析.M.北京.機械工業(yè)出版社.2009.1.附錄：L=1;M=40;N=1600;alfa=1;lambda=0.5;%網(wǎng)格比%*%h=L/M;%空間步長x=0:h:L;x=x'

7、;tao=lambda*h2/alfa;%時間步長tm=N*tao;%熱傳導(dǎo)的總時間tm%tm=0.1;t=0:tao:tm;t=t'%計算初值和邊值U=zeros(M+1,N+1);U(:,1)=sin(pi*x);U(1,:)=0;U(M+1,:)=0;%*用差分法求出溫度U，與桿長L，時間T的關(guān)系*%for k=1:N j=2; while j<=M U(j,k+1)=lambda*U(j+1,k)+(1-2*lambda)*U(j,k)+lambda*U(j-1,k); j=j+1; endend length(U);%*設(shè)置立體網(wǎng)格*%for i=1:N+1 X(:,i

8、)=x;endfor j=1:M+1 Y(j,:)=t;endmesh(X,Y,U);legend('數(shù)值解');xlabel('X');ylabel('T');zlabel('U');z=zeros(M+1,N+1);for j=1:M+1 for k=1:N+1 z(j,k)=exp(-pi*pi*t(k)*sin(pi*x(j); endend%mesh(x,t,z')legend('解析解');xlabel('X');ylabel('T');zlabel('Z');for j=1:M+1 for k=1:N+1 y(j,k)=abs