本文格式為Word版，下載可任意編輯——畢業設計二維熱傳導方程有限差分法的MATLAB實現蘇佳園：二維熱傳導方程有限差分法的MATLAB實現

l所構成，其中xj?j?x?h,?x?h?J;tn?n?t?n?.

3.2插值函數的選擇

選擇不同的插值函數對偏微分方程進行估計，可得到不同的差分方程，進而穩定性和精度會有所不同。

用Taylor級數展開方法是最常用的方法，下面建立差分格式的同時引入一些基本概念及術語。

我們主要從對流方程的初值問題

?u??u?a?0,x?R,t?0,?(3.1)?x??t??u(x,0)?g(x),x?R,和擴散方程的初值問題

?u?2u?a2,?t?xu(x,0)?g(x),（其中a?0）進行探討。

x?R,t?0,(3.2)x?R.

假定偏微分方程初值問題的解u(x,y)是充分官話的，由Taylor級數展開有

u(xj,tn?1)?u(xj,tn)?un????(?t),?t?tj?u(xj,tn?1)?u(xj,tn?1)2??un???(?t),2?t?tj?u(xj?1,tn)?u(xj,tn)?(3.3)?un???(h),h?xj?u(xj,tn)?u(xj?1,tn)??un???(h),?h?xj?u(xj?1,tn)?u(xj?1,tn)2?un???(h),?2h?xj?u(xj?1,tn)?2u(xj,tn)?u(xj?1,tn)n2?2u????(h).22h?xj?

????????????n其中???n或用，表示看括號內的函數在節點(xj,tn)處取的值。利用（3.3）表達式???jj中的第1式和第3式有

u(xj,tn?1)?u(xj,tn)?

?au(xj?1,tn)?u(xj,tn)h6

?u?un?[?a]j??(??h).?t?x假使u(x,t)是滿足偏微分方程（3.1）的光滑解，則

[?u?u?a]nj?0.?t?x由此看一看出，偏微分方程（3.1）在u(xj,tn)處可以近似的用下面的方程來代替

?1nun?ujjnun?uj?1j?n?ah?0,(3.4)

j?0,?1,?2,?,n?0,1,2?.

其中uj為u(xj,tn)的近似值。（3.4）式稱為迫近微分方程（3.1）的有限差分方程或簡稱差分方程。

差分方程再加上初始條件的離散型式就可以按時間逐層推進，算出各層的值。差分格式（3.4）和初始條件的離散形式結合在一起構成了一個差分格式。

3.3方程組的建立

將離散后的差分方程轉化為方程組的形式，便于求解。

3.4方程組的求解

利用矩陣的解法求解方程組，再用MATLAB對矩陣求解方法進行程序化，以便對以后類似的方程進行求解。

隱式差分格式方程矩陣化后，得到的矩陣是嚴格的對角占優三對角矩陣，我們可以根據線性方程組的求解方法對其求解。其中這要應用的是追趕法，追趕法對于此類線性方程組的求解十分便利，用MATLAB對追趕法進行編程，就可以輕松實現矩陣的求解，進而解出差分方程的近似解。

7

蘇佳園：二維熱傳導方程有限差分法的MATLAB實現

第4章二維熱傳導方程

4.1網格剖分

在區域D:?(x,y,t)0?x?X,0?y?Y,0?t?T?中，我們設二維熱傳導方程的初始值和邊界條件如下：

2u?2u??u???t?a(2?2),0?x?X,0?y?Y,0?t?T,?x?y???u(x,y,0)?u0(x,y),0?x?X,0?y?Y??u(x0,y,t)?f(y,t),u(xX,y,t)?g(y,t),??u(x,y0,t)?h(x,t)，u(x,yY,t)?k(x,t)，其中a為正常數。

（4.1）

通過已知方程，建立一個關于時間和步長的函數，這樣就把初始區域劃分為一個網格圖。先將定義域

D:?(x,y,t)0?x?X,0?y?Y,0?t?T?

剖分為網格

Dh?{(xj,yh,tn)xj?j?x,j?0,1,?,J,J?x?X;yl?l?y,l?0,1?,J,J?y?Y;tn?n?t,n?0}.其中?t?

T為時間步長，?x?X,?y?Y分別為x軸和y軸的空間步長。NMxMy4.2穩定性分析

利用有限差分格式進行計算時是按時間層逐層推進的。那么計算第n?1層上的值時要用到第n層上計算出來的結果值，而計算第n層結果值時的舍入誤差必然會影響到第n?1層的值。從而就要分析這種誤差傳播的狀況。希望誤差不至于越來越大，以至掩蓋差分格式的解的面貌，這便是穩定性問題。

我們先考慮一維差分格式

?1nnnun?u?a?(u?ujjj?1j)(4.2)

的穩定性，其中??

?t為網格比，假設a?0差分格式從初層開始計算，當時始數據?x8

存在誤差時考察這個誤差在以后計算中的在傳播狀況。為便利起見，不考慮計算過程中的舍入誤差。及確定初始數據誤差絕對值為?，則差分格式在(xj,tn)處的誤差為

nnmm?n(?1)n?m??Cm(1?a?)(?a?)m?0nn???Cm(1?a?)m(a?)m?nm?0?(1?a?)n?.于是，對于固定網格比?及a?0的狀況，差分格式的解的誤差隨時間步長的步數

n的增加而增加。初始數據的誤差將必定掩蓋了差分格式的解的面貌，所以我們認為

差分格式（4.2）時不穩定的。

差分格式的穩定性不僅與差分格式本身有關，而且還與網格比的大小有關。差分格式的穩定性在差分方法的研究中具有特別的意義，我們再做進一步的表達。

定義4.1[8]為了度量誤差及其他應用，引入范數

nhu00n???n2???(uj)h?,?j????n12設uj有一個誤差?j，則uj就有誤差?j。假使存在一個正常數K,使得當

?t??t0,n?t?T時，一致的有

?n?K?0,

則稱差分格式是穩定的。

差分格式一旦具有穩定性，就可以用差分格式計算出偏微分方程的近似解來。一維熱傳導方程的各種類型的差分格式可以推廣到二維熱傳導方程，利用向前差分格式

?1nun?ujjnnun?2u?uj?1jj?1?t對（4.1）式進行離散，引入記號

?h2

nnn?x2un?u?2u?ujlj?1,ljlj?1,l,nnn?y2un?u?2u?ujlj,l?1jlj,l?1.（4.3）

其中unjl為差分方程在節點(j,l,n)的計算值。

差分格式

9

蘇佳園：二維熱傳導方程有限差分法的MATLAB實現

?1nun?ujljl

?t?a(?x2unjl?x2?2n?yujl?y2),（4.4）

利用Taylor級數展開易得差分格式（4.4）的截斷誤差為?(?t??x2??y2)。

t為便利穩定性的判斷，設?x??y?h，令???2，為網格比。即改寫為：

hn?1n2n2nujl?ujl?a?(?xujl??yujl),（4.5）

用Fourier方法來分析（4.4）式的穩定性。令

nik1jhik2lhun?vee,jl把此式帶入（4.5）式中有

vn?1n?1?2a?(cosk1h?1)?2a?(cosk2h?1)v,

??因此差分格式（4.5）的增長因子是

2k1h2k2hG(?t,k)?1?4a?(sin?sin),

22其中k?(k1,k2)。假使a??件是

14，則有G(?t,k)?1由此得出差分格式（4.4）的穩定性條

a??14.

簡單看出在二維狀況下采用這樣的顯式格式是不適合的，為此我們再轉向考慮隱式格式。用一維向后差分格式

nn?1uj?uj?tnnnuj?1?2uj?uj?1?2h

的直接推廣是

n?1unj