馮立偉

(沈陽化工大學(xué)數(shù)理系，遼寧沈陽110142)

許多工程問題需要研究熱量在物體內(nèi)部的傳導(dǎo)情況或某種物質(zhì)在液體中的擴散情況，因此研究熱傳導(dǎo)問題特別是非穩(wěn)態(tài)熱傳導(dǎo)問題十分重要.目前熱傳導(dǎo)方程已有多種求解格式［1－2］.MATLAB是目前最流行、應(yīng)用最廣泛的科學(xué)和工程計算軟件.MATLAB基于矩陣運算，具有強大的數(shù)值運算能力和圖形可視化能力，是方便實用、功能強大的數(shù)學(xué)軟件［3］.用MATLAB求解常微分方程已有大量的研究［4－6］.王飛等介紹了如何使用MATLAB實現(xiàn)有限差分法求解微分方程［7］.高理平等給出了對兩點邊值問題有限元方法的程序［8］.李燦等對熱傳導(dǎo)問題的MATLAB數(shù)值計算進行了討論［9］.本文討論求解一維熱傳導(dǎo)方程幾種不同差分格式的MATLAB編程方法，并使用算例進行檢驗和對結(jié)果進行分析.

u(x，t)表示在t時刻物體內(nèi)部坐標(biāo)為x處的溫度，a是熱傳導(dǎo)系數(shù)，f(x，t)為熱源.

1 熱傳導(dǎo)方程差分格式

1.1 差分格式的建立

首先對x-t平面進行網(wǎng)格剖分.分別取h，τ為空間步長與時間步長，用2族平行直線xj= jh，tk=kτ將矩形區(qū)域［0，T］×［0，L］分割成矩形網(wǎng)格.

顯式格式:

隱式格式:

Crank-Nicolson格式:

Du Fort Frankel格式:

由Taylor公式容易得出:它們都與一維熱傳導(dǎo)方程相容，其截斷誤差分別為O(τ+h2)，O(τ+ h2)，O(τ2+h2)和O(τ2+h2)［1］.

1.2 初、邊值條件的處理

對定解條件進行離散化.由初始條件及第一類邊界條件，可直接得到:

1.3 穩(wěn)定性分析

使用Fourier方法可知，當(dāng)r≤1/2時顯式格式穩(wěn)定.

隱式格式變?yōu)?

由Fourier方法可得隱式格式恒穩(wěn)定.

Crank-Nicolson格式變?yōu)?

對于r＞0恒有增長因子|G|≤1，恒穩(wěn)定.

Du Fort Frankel格式變?yōu)?

由Fourier方法可得Du Fort Frankel格式恒穩(wěn)定.

2 差分格式的求解

顯式格式:將(6)式與離散化的初邊值條件結(jié)合，得到求解此問題的差分方程組:

由于初始時間層上的u值為已知，由(10)式即可算出u在第一層各個節(jié)點處的近似值u1j.重復(fù)使用此式，可以逐層計算出所有的

隱式格式:將(7)式與離散化的初邊值條件聯(lián)立，得差分方程組:

將上述方程組改寫成矩陣形式

此方程組是三對角方程組，且系數(shù)矩陣嚴(yán)格對角占優(yōu)，故解存在唯一.通過在每一時間層上求解一個這樣的線性方程組得到在各個時刻各個網(wǎng)格點上的函數(shù)值.

Crank-Nicolson格式:將(8)式與離散化的初邊值條件聯(lián)立并整理，得差分方程組:

此方程組的系數(shù)矩陣嚴(yán)格對角占優(yōu)，差分方程組解存在唯一.

程序如下:

Du Fort Frankel格式:將(7)式與初始條件及第一類邊界條件式聯(lián)立并整理將格式改寫為:

由于此格式是3層格式，需要事先知道前面2個時間層上的解，第1層上的解通過對初始條件離散可得，第2層上的解使用前面的任意1種2層差分格式得到，再用此格式求解其余時間層上的解.程序核心代碼如下:

3 數(shù)值實驗

使用前面4種不同差分格式求解下面熱傳導(dǎo)方程的初邊值問題.

此問題的真解為u(x，t)=e－π2tsin(πx).

定義誤差e=u－U，

取空間步長h=0.1，求解結(jié)果如圖1～圖6所示.

圖1 r=1時各種差分方法的誤差‖ek‖2Fig.1 r=1 errors of several difference schemes‖ek‖2

圖2 r=時各種差分方法的誤差‖ek‖2Fig.2 r=errors of several difference schemes‖ek‖2

圖3 r=時各種差分方法的誤差‖ek‖2Fig.3 r=errors of several difference schemes‖ek‖2

圖4 r=時各種差分方法的誤差‖ek‖2Fig.4 r=errors of several difference schemes‖ek‖2

圖5 r=時各種差分方法的誤差‖ek‖2Fig.5 r=errors of several difference schemes‖ek‖2

圖6 r=0.67時各種差分方法的誤差‖ek‖2Fig.6 r=0.67 errors of several difference schemes‖ek‖2

4 結(jié)論

給出了使用MATLAB求解熱傳導(dǎo)方程幾種差分格式的方法和部分主要程序，通過數(shù)值實驗看到Du Fort Frankel格式和 Crank-Nicolson格式是誤差較小且實用的方法.

［1］ 胡健偉，湯懷民.微分方程數(shù)值方法［M］.2版.北京:科學(xué)出版社，2007:131－138.

［2］ Morton K W，Mayers D F.Numerical Solution of Partial Differential Equations［M］.2thed.Cambridge UK:Cambridge University Press，2005:19－26.

［3］ 鄭阿奇.MATLAB實用教程［M］.北京:電子工業(yè)出版社，2005:175－182.

［4］ 單毅.常微分方程的MATLAB解法［J］.武漢大學(xué)學(xué)報(工學(xué)版)，2003，36(z2):150－152.

［5］ 何雙.MATLAB在常微分方程初值問題的應(yīng)用［J］.長春師范學(xué)院學(xué)報(自然科學(xué)版)，2005，24 (3):17－19.［6］ 唐洪浪，桂現(xiàn)才.用MATLAB符號工具箱編程求常微分方程的通解［J］.洛陽師范學(xué)院學(xué)報，2005，24(2):81－84.

［7］ 王飛，裴永祥.有限差分方法的MATLAB編程［J］.新疆師范大學(xué)學(xué)報(自然科學(xué)版)，2003，22 (4):22－27.

［8］ 高理平，楊光.兩點邊值問題有限元方法的程序?qū)崿F(xiàn)與計算分析［J］.山東師范大學(xué)學(xué)報(自然科學(xué)版)，2009，24(1):1－5.

［9］ 李燦，高彥棟，黃素逸.熱傳導(dǎo)問題的matlab數(shù)值計算［J］.華中科技大學(xué)學(xué)報(自然科學(xué)版)，2002，30(9):91－93.