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

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

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

4、果如表(1)數值解精確解相對誤差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)數值解與精確解的比較由表(1)我們可以看出，精確解和數值解的絕對誤差在以內，因此可以得出，在分割M=40，N=1600下，該有限差分方法對方程(1.3.1)是收斂和穩定的。下面，我們比較在不同的分割下對有限差分算法精度的影響

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

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

7、;tao=lambda*h2/alfa;%時間步長tm=N*tao;%熱傳導的總時間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的關系*%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);%*設置立體網格*%for i=1:N+1 X(:,i

8、)=x;endfor j=1:M+1 Y(j,:)=t;endmesh(X,Y,U);legend('數值解');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