蘇佳園：二維熱傳導方程有限差分法的MATLAB實現(xiàn)PAGE28PAGE29第1章前言1.1問題背景在史策教授的《一維熱傳導方程有限差分法的MATLAB實現(xiàn)》和曹剛教授的《一維偏微分方程的基本解》中，對偏微分方程的解得MATLAB實現(xiàn)問題進行過研究，但只停留在一維中，而實際中二維和三維的應用更加廣泛。諸如粒子擴散或神經(jīng)細胞的動作電位。也可以作為某些金融現(xiàn)象的模型，諸如布萊克-斯科爾斯模型與Ornstein-uhlenbeck過程。熱方程及其非線性的推廣形式也被應用與影響分析。在科學和技術(shù)發(fā)展過程中，科學的理論和科學的實驗一直是兩種重要的科學方法和手段。雖然這兩種科學方法都有十分重要的作用，但是一些研究對象往往由于他們的特性（例如太大或太小，太快或太慢）不能精確的用理論描述或用實驗手段來實現(xiàn)。自從計算機出現(xiàn)和發(fā)展以來，模擬那些不容易觀察到的現(xiàn)象，得到實際應用所需要的數(shù)值結(jié)果，解釋各種現(xiàn)象的規(guī)律和基本性質(zhì)。科學計算在各門自然科學和技術(shù)科學與工程科學中其越來越大的作用，在很多重要領(lǐng)域中成為不可缺少的重要工具。而科學與工程計算中最重要的內(nèi)容就是求解科學研究和工程技術(shù)中出現(xiàn)的各種各樣的偏微分方程或方程組。解偏微分方程已經(jīng)成為科學與工程計算的核心內(nèi)容，包括一些大型的計算和很多已經(jīng)成為常規(guī)的計算。為什么它在當代能發(fā)揮這樣大的作用呢？第一是計算機本身有了很大的發(fā)展；第二是數(shù)值求解方程的計算法有了很大的發(fā)展，這兩者對人們計算能力的發(fā)展都是十分重要的。1.2問題現(xiàn)狀近三十年來，解偏微分方程的理論和方法有了很大的發(fā)展，而且在各個學科技術(shù)的領(lǐng)域中應用也愈來愈廣泛，在我國，偏微分方程數(shù)值解法作為一門課程，不但在計算數(shù)學專業(yè)，而且也在其他理工科專業(yè)的研究生的大學生中開設(shè)。同時，求解熱傳導方程的數(shù)值算法也取得巨大進展，特別是有限差分法方面，此算法的特點是在內(nèi)邊界處設(shè)計不同于整體的格式，將全局的隱式計算化為局部的分段隱式計算。而且精度上更好。目前，在歐美各國MATLAB的使用十分普及。在大學的數(shù)學、工程和科學系科，MATLAB被用作許多課程的輔助教學手段，MATLAB也成為大學生們必不可少的計算工具，甚至是一項必須掌握的基本技能。在我國，MATLAB在各大專院校的應用日益普遍，許多專業(yè)已把MATLAB作為基本計算工具。在科研機構(gòu)和工業(yè)界，MATLAB正得到越來越廣泛的應用。MATLAB具有強大的圖形繪制功能，為科學計算和圖形處理提供了很大的方便。我們只需制定的繪圖方式，再提供繪圖數(shù)據(jù)，有程序指令就可以得到形象、直觀的圖形結(jié)果。因此，近些年越來越多的人開始使用MATLAB來求解數(shù)值計算和圖形處理技術(shù)，我們也可以繪制出熱傳導方程數(shù)值解的二維、三維圖形，從而可以更好的理解熱傳導方程的意義。1.3問題解決目前，對于求解偏微分方程有很多方法，但差分法和有限元離散法式主要解決問題的兩種方法。一般來說，用差分法來接偏微分方程，解得得結(jié)果就是方程的準確解函數(shù)再借點上的近似值。而用變分近似的方法求解，是將近似解表示成有限維子空間中基函數(shù)的線性組合。有限元法也是基于變分原理，由于選擇了特殊的基函數(shù)，使它能適用于一般的區(qū)域。這種基函數(shù)是與區(qū)域的剖分有關(guān)的，近似解表示為基函數(shù)的線性組合，二線性組合中的系數(shù)，又是剖分節(jié)點上或其導數(shù)的近似值。有關(guān)一維熱傳導方程的有限差分法求解的MATLAB實現(xiàn)，西安建筑科技大學的史策教授已經(jīng)解決，本文借鑒史老師的求解思想，對二維熱傳導方程進行轉(zhuǎn)換，再對解法編程實現(xiàn)，從而進一步對熱傳導方程進行探討。二維熱傳導方程求解在現(xiàn)實生活中的應用也更加廣泛，所以有很好的現(xiàn)實意義。第2章預備知識定義2.1[8]含有未知函數(shù)的偏導數(shù)的方程稱為偏微分方程。定義2.2[8]方程 稱熱傳導方程（或擴散方程）。其中，是固體的傳熱過程中在處、時刻的溫度。系數(shù)稱為熱傳導系數(shù)，當時，方程為其中，為維數(shù)。定義2.3[8]在特定條件下求解方程的解。這樣的條件成為定解條件。給出了方程和定結(jié)條件，就構(gòu)成了定解問題。定義2.4[1]一般說，邊界條件有下列形式其中為邊界的外法向?qū)?shù)。有如下幾種特殊形式(1)Dirichlet（或第一類）條件：即值給定；(2)Neumann（或第二類）條件：.即的外法向?qū)?shù)給定；(3)Robbins（或第三類）條件：。定義2.5[8]只有出事條件而沒有邊界條件的定解問題。定義2.6[8]只有邊界條件而沒有初值條件的定解問題。定義2.7[8]既有邊值條件又有初值條件的定解問題。定義2.7[8]定義在上的函數(shù)的一個關(guān)系式，設(shè)，有關(guān)系式以上變換稱為Fourier變換。其中是虛數(shù)單位。定義2.9[8]由第個時間層推進到第個時間層時差分方程提供了逐點直接計算的表達式，我們稱次差分方程為顯式格式。定義2.10[8]有限差分格式在新的時間層上包含有多于一個的節(jié)點，這種有限差分格式稱為隱式格式。定義2.11[11]稱為向前差分。定義2.12[11]稱為向后差分。定義2.13[11]稱為中心差分。定義2.14[11]用微分方程的解代替差分方程的全部近似解，這樣得到的方程兩邊的差就是截斷誤差。定理2.1[8]給定一個適定的線性初值問題以及與其相容的差分格式，則差分格式的穩(wěn)定性是差分格式收斂性的充要條件。第3章求解二維熱傳導方程的基本思想基本思想是把連續(xù)的定解區(qū)域用有限個離散點構(gòu)成網(wǎng)格來代替，這些離散點稱作網(wǎng)格的節(jié)點；把連續(xù)定解區(qū)域上的連續(xù)變量的函數(shù)用在網(wǎng)格上定義的離散變量函數(shù)的近似；把原方程和定解條件中的微商用差商來近似，積分用積分來近似，于是原微分方程和定解條件近似的代之以代數(shù)方程，即優(yōu)先差分方程組，解此方程組就可以得到原問題在離散點上的近似解。下面是有限差分法數(shù)值計算的基本步驟：3.1區(qū)域的離散用有限差分方法求解偏微分方程問題必須把連續(xù)問題進行離散化。為此首先要對求解區(qū)域給出網(wǎng)格剖分，由于求解的問題不同，因此求解區(qū)域也不盡相同。下面用例子來說明不同區(qū)域的剖分離散。并引入一些常用術(shù)語。例3.1雙曲型和拋物型方程的初值問題，求解區(qū)域是我們在的上半平面畫出兩族平行于坐標軸的直線，把上班平面分成矩形網(wǎng)格。其交點稱為節(jié)點（或網(wǎng)格點）?？稍O(shè)距離，稱其為空間步長，平行線的距離按具體問題而定?？稍O(shè)距離，稱其為時間步長。這樣兩族網(wǎng)格線可以寫作網(wǎng)格節(jié)點有時記為。例3.2雙曲型和拋物型方程的初邊值問題，設(shè)求解區(qū)域是這個區(qū)域的網(wǎng)格由平行于軸的直線族與平行于軸的直線族所構(gòu)成，其中3.2插值函數(shù)的選擇選擇不同的插值函數(shù)對偏微分方程進行估計，可得到不同的差分方程，進而穩(wěn)定性和精度會有所不同。用Taylor級數(shù)展開方法是最常用的方法，下面建立差分格式的同時引入一些基本概念及術(shù)語。我們主要從對流方程的初值問題(3.1)和擴散方程的初值問題(3.2)（其中）進行討論。假定偏微分方程初值問題的解是充分官話的，由Taylor級數(shù)展開有(3.3)其中或用，表示看括號內(nèi)的函數(shù)在節(jié)點處取的值。利用（3.3）表達式中的第1式和第3式有.如果是滿足偏微分方程（3.1）的光滑解，則由此看一看出，偏微分方程（3.1）在處可以近似的用下面的方程來代替(3.4)其中為的近似值。（3.4）式稱為逼近微分方程（3.1）的有限差分方程或簡稱差分方程。差分方程再加上初始條件的離散型式就可以按時間逐層推進，算出各層的值。差分格式（3.4）和初始條件的離散形式結(jié)合在一起構(gòu)成了一個差分格式。3.3方程組的建立將離散后的差分方程轉(zhuǎn)化為方程組的形式，便于求解。3.4方程組的求解利用矩陣的解法求解方程組，再用MATLAB對矩陣求解方法進行程序化，以便對以后類似的方程進行求解。隱式差分格式方程矩陣化后，得到的矩陣是嚴格的對角占優(yōu)三對角矩陣，我們可以根據(jù)線性方程組的求解方法對其求解。其中這要應用的是追趕法，追趕法對于此類線性方程組的求解非常方便，用MATLAB對追趕法進行編程，就可以輕松實現(xiàn)矩陣的求解，進而解出差分方程的近似解。第4章二維熱傳導方程4.1網(wǎng)格剖分在區(qū)域中，我們設(shè)二維熱傳導方程的初始值和邊界條件如下：其中為正常數(shù)。通過已知方程，建立一個關(guān)于時間和步長的函數(shù)，這樣就把初始區(qū)域劃分為一個網(wǎng)格圖。先將定義域剖分為網(wǎng)格其中為時間步長，分別為軸和軸的空間步長。4.2穩(wěn)定性分析利用有限差分格式進行計算時是按時間層逐層推進的。那么計算第層上的值時要用到第層上計算出來的結(jié)果值，而計算第層結(jié)果值時的舍入誤差必然會影響到第層的值。從而就要分析這種誤差傳播的情況。希望誤差不至于越來越大，以至掩蓋差分格式的解的面貌，這便是穩(wěn)定性問題。我們先考慮一維差分格式(4.2)的穩(wěn)定性，其中為網(wǎng)格比，假設(shè)差分格式從初層開始計算，當初始數(shù)據(jù)存在誤差時考察這個誤差在以后計算中的在傳播情況。為方便起見，不考慮計算過程中的舍入誤差。及確定初始數(shù)據(jù)誤差絕對值為，則差分格式在處的誤差為于是，對于固定網(wǎng)格比及的情況，差分格式的解的誤差隨時間步長的步數(shù)的增加而增加。初始數(shù)據(jù)的誤差將必定掩蓋了差分格式的解的面貌，所以我們認為差分格式（4.2）時不穩(wěn)定的。差分格式的穩(wěn)定性不僅與差分格式本身有關(guān)，而且還與網(wǎng)格比的大小有關(guān)。差分格式的穩(wěn)定性在差分方法的研究中具有特別的意義，我們再做進一步的敘述。定義4.1[8]為了度量誤差及其他應用，引入范數(shù)設(shè)有一個誤差，則就有誤差。如果存在一個正常數(shù),使得當時，一致的有則稱差分格式是穩(wěn)定的。差分格式一旦具有穩(wěn)定性，就可以用差分格式計算出偏微分方程的近似解來。一維熱傳導方程的各種類型的差分格式可以推廣到二維熱傳導方程，利用向前差分格式對（4.1）式進行離散，引入記號（4.3）其中為差分方程在節(jié)點的計算值。差分格式（4.4）利用Taylor級數(shù)展開易得差分格式（4.4）的截斷誤差為。為方便穩(wěn)定性的判斷，設(shè)，令，為網(wǎng)格比。即改寫為：（4.5）用Fourier方法來分析（4.4）式的穩(wěn)定性。令把此式帶入（4.5）式中有因此差分格式（4.5）的增長因子是其中。如果，則有由此得出差分格式（4.4）的穩(wěn)定性條件是容易看出在二維情況下采用這樣的顯式格式是不適合的，為此我們再轉(zhuǎn)向考慮隱式格式。用一維向后差分格式的直接推廣是（4.6）此格式的截斷誤差仍為，仍用Fourier方法來分析這個格式的穩(wěn)定性，仿前可以得出其增長因子是：因此有，即差分格式（4.6）是絕對穩(wěn)定的。為了提高精度，對微分方程（4.1）也可以用Crank-Nicolson型差分格式，這也是一維問題的直接推廣。其格式可寫為（4.7）這也是二階精度格式，其增長因子是因此，對任何都有，所以(4.7)式也是絕對穩(wěn)定的?，F(xiàn)在考慮一下隱式格式（4.6）式和（4.7）式的求解方法。我們知道，在一位格式形成的方程組是系數(shù)矩陣為三對角矩陣的線性代數(shù)方程組，因此用追趕法很容易求解。而對于（4.6）式和（4.7）式導出的系數(shù)矩陣不是三對角矩陣，因此求解就不容易了。我們對于顯式格式和隱式格式的分析知道，在實際使用上都受到限制，因此構(gòu)造每層計算量不大的絕對穩(wěn)定的格式就成為一個具有現(xiàn)實意義并很有興趣的問題。在一維中，隱式格式是絕對穩(wěn)定并可用追趕法很容易求解。由此產(chǎn)生了下面將使用的交替方程隱式格式。它具有絕對穩(wěn)定、容易求解和有相當精度的特點。4.3建立方程組我們在構(gòu)造微分方程（4.1）的隱式格式和顯式格式中，對和做了同樣的處理，即同時在第層或第層取值。為了構(gòu)造一維形式的隱式格式，對二階導數(shù)用在第層上用未知的二階中心差分來代替，而則用在第層上用已知的二階中心差分來代替，這樣得到的方程組在僅方向是隱式的。比較容易求解，用追趕法就可以了。同理，為對稱起見，在下一時間層上重復上述步驟，即又僅在方向是隱式的，對方向是顯式的。這樣相鄰的兩個時間層合并起來構(gòu)成一個差分格式。故用多次追趕法就可以解出了。我們解向后差分格式的方程。令，則（4.6）式變?yōu)?（4.8）（4.8）在x方向是隱式時，代入（4.3），變形為（4.9）（4.9）化為矩陣形式：（4.10）（4.8）在y方向是隱式時，代入（4.3），可變形為（4.11）（4.11）化為矩陣形式：（4.12）第5章MATLAB編程5.1追趕法對于前面求出的（4.10）和（4.12）矩陣，我們稱之為三對角矩陣。在一些實際問題中，例如解常微分方程的邊值問題，接熱傳導方程一級船體數(shù)學放樣中建立三次樣條函數(shù)等，都會要求解系數(shù)矩陣為對角占優(yōu)的三對角方程組，解三對角矩陣我們有一種特殊的方法叫追趕法。例如：求解線性方程組，其中為三對角矩陣。首先對矩陣做如下三角分解：其中為下三角矩陣，為上三角矩陣，且主對角線元素全為1.求解等價于求解兩個三角形方程組。先求解的解，再求解的解。追趕法公式實際上就是把高斯消去法用到求解三對角方程組上去的結(jié)果。這是由于特別簡單，因此使得求解的計算公式非常簡單，而且計算量很小。另外追趕法的計算公式中不會出現(xiàn)中間結(jié)果數(shù)量級的巨大增長和舍入誤差的嚴重積累。最后在MATLAB中編程實現(xiàn)。追趕法程序：見附錄。例5.1用追趕法求解線性方程組解：在MATLAB命令窗口中輸入求解程序：>>A=[2-100;-13-20;0-24-3;00-35]A=2-100-13-200-24-300-35 >>f=[61-21]’;>>x=chase(A,f)輸出計算結(jié)果為：x=91210.即追趕法求出了方程組的解為：5.2二維熱傳導方程有限差分解的MATLAB編程程序見附錄。第6章數(shù)值舉例分析6.1算例為了便于觀察和表示微分方程的解，我們將最終解用圖像表示出來。例6.1求滿足下列條件的二維熱傳導方程的解：在MATLAB命令窗口中輸入求解程序：a=0.5;u_xy0=inline('0','x','y');u_xyt=inline('x^2*sin(y)-y^2*sin(x)','x','y','t');D=[0,5,0,5];T=1000;Mx=50;My=50;N=50;[u,x,y,t]=sjy(a,D,T,u_xy0,u_xyt,Mx,My,N);mesh(x,y,u)xlabel('x')ylabel('y')zlabel('U')輸出計算結(jié)果為：rx=1.0000e+003ry=1.0000e+003此解即為差分方程的近似解。例6.2求解如下二維熱傳導方程在MATLAB命令窗口中輸入求解程序：a=10;u_xy0=inline('0','x','y');u_xyt=inline('x*sin(pi*y)-y*sin(pi*x)','x','y','t');D=[0,2,0,4];T=500;Mx=80;My=80;N=40;[u,x,y,t]=sjy(a,D,T,u_xy0,u_xyt,Mx,My,N);mesh(x,y,u)xlabel('x')ylabel('y')zlabel('t')輸出計算結(jié)果為：rx=2.0000e+005ry=5.0000e+0046.2運行結(jié)果分析利用MATLAB的工具箱PDETOOL求解方程結(jié)果為：例6.1精確解例6.2精確解對上述數(shù)值解和精確解進行比較，將兩解得圖畫在一個坐標面內(nèi)，可以更加直觀的觀察圖形的區(qū)別，若兩圖形相差比較大則說明此算法還有不對或不完善的地方；若兩圖形重合度比較高，則說明此算法和編程過程比較好，能真實的反應偏微分方程的解的情況。如下圖：從上圖可以看出結(jié)果基本一致，重合度比較高，即說明計算差分方程的的結(jié)果比較真實，能夠很好的反應真實解的情況。若數(shù)值解的誤差在一定的范數(shù)下滿足不等式（為初始層的誤差），則熱傳導方程的數(shù)值解是穩(wěn)定的，可從上圖看出，兩圖幾乎重合，說明數(shù)值解和精確解的誤差非常小。第7章結(jié)束語本篇文章利用MATLAB數(shù)學軟件來求解二維熱傳導方程，通過對區(qū)域的剖分，把微分方程進行離散化，在選擇隱式差分格式對微分方程進行轉(zhuǎn)化，得到隱式差分方程，再把差分方程轉(zhuǎn)換成矩陣形式求解，由于矩陣形式是嚴格對角占優(yōu)的三對角矩陣，通過追趕法輕松的求解。利用MATLAB對追趕法的思想進行編程實現(xiàn)。因此，可以快速的求解類似的二維熱傳導方程。而且誤差非常小，可以控制在有效的范圍內(nèi)。以上算例表明了此方法的可行性和穩(wěn)定性。相比一維熱傳導方程有限差分法的MATLAB實現(xiàn)，二維熱傳導方程的求解更加深入，實際應用性更強。但求解的基本思想不變。致謝本論文是在導師王其林老師和張學富老師的悉心指導下完成的。導師淵博的專業(yè)知識，嚴謹?shù)闹螌W態(tài)度，精益求精的工作作風，誨人不倦的高尚師德，嚴于律己，寬于待人的崇高風范和樸實無華、平易近人的人格魅力對我影響深遠。不僅使我在學術(shù)上有了更進一步的深入、專業(yè)知識更加鞏固，最重要的是在我以后的工作和學習中明白了對于待人接物和為人處事的道理，著實受益匪淺。這將會是我人生經(jīng)歷中一筆重要的財富。文章從選題到完成，每一步都有導師的津津指導，傾注了導師的大量心血。在此，謹向?qū)煴硎境绺叩木匆夂陀芍缘母兄x。本文的順利完成，離不開各位老師、同學、親人和朋友的關(guān)心與幫助。非常感謝鄒昌文老師的指導和幫助；感謝負責管理實驗室的老師和同學的幫助；感謝室友的幫助和勉勵；感謝哥哥蘇佳斌的關(guān)心和編程方面的指導；感謝爸爸媽媽的多年教育及對我學業(yè)的大力支持；感謝好朋友孫斌的幫助和支持。在此，對以上幾位以及其他在背后默默支持和幫助我的老師、同學和朋友表示深深地感謝。沒有他們的幫助和支持是不可能完成我的學士論文的。愿他們好人一生平安、幸福。重慶交通大學理學院數(shù)學與應用數(shù)學06級2班蘇佳園 2010年6月參考文獻[1]陸金甫，關(guān)治.偏微分方程的數(shù)值解法[M]，北京：清華大學出版社，2006。[2]史策.熱傳導方程有限差分的MATLAB實現(xiàn)[J].西安建筑科技大學學報，2009，24(4)：27-29.[3]馬明書，馬菊意，谷淑敏.高維熱傳導方程的高精度分支顯格式[J]，河南師范大學學報(數(shù)學與信息科學學院)，2008，23(3)：446-452.[4]王正林，何倩，龔純.精通MATLAB科學計算[M],北京：電子工業(yè)出版社，2009。[5]王飛，裴永祥.有限差分法的MATLAB編程[J].新疆師范大學學報（自然科學版），2003，22(4)：21-27.[6]李先枝.熱傳導方程差分解法的最佳網(wǎng)絡[J],河南大學學報（自然科學版）,2004，34(3)：16-18.[7]萬正蘇，方春華，張再云.關(guān)于熱傳導方程有限差分區(qū)域分解并行算法精度的注記[J].湖南理工學院學報（自然科學版），2007，20(3)：12-14.[8]曹鋼，王桂珍，任曉榮.一維熱傳導方程的基本解[J].山東輕工業(yè)學院學報，2005，19(4)：76-80.[9]李傅山，孔翠翠.熱傳導方程的基本解與正態(tài)分布密度函數(shù)[J]，曲阜師范大學數(shù)學科學學院，2006，77(2)：42-43.[10]陽衡.一類擬線性拋物線型偏微分方程的解[J]，湖南環(huán)境生物職業(yè)技術(shù)學院學報，2009，3(17)：17-19.[11]蔡國梁，李丹.二維熱傳導方程的非古典對稱和相容性[J]，江蘇大學理學院計算科學系學報，2008，29(3)：273-276.[12]張正林，王遠弟.非局部邊界條件下的拋物型偏微分方程組[J]。上海大學理學院數(shù)學系學報，2008，22(1)：21-30.[13]劉相國，周宏宇，張海燕，閔濤.二維變系數(shù)熱傳導方程初邊值問題的交替方向隱格式[J],西安理工大學理學院，2007，27(2)：199-204.[14]谷建濤，佟玉霞，付景紅.一階偏微分方程的Mathematica和Matlab解法比較[J],河北理工大學學報，2008，30(2)：40-42.[15]LILi-kangetal.NumericalSolutionofDierentialEquations[M].Shanghai:FudanUniversity[16]SUNHong-lie.Aclassofhighaccuracyexplicitdifferenceschemeforsolvingheat-conductionequationsofmulti-dimension[J].AppliedMathematicsJournalofUniversities,Edition,1999,14(4):427-432.附錄追趕法程序functionx=chase(A,f)n=rank(A);b=zeros(n,1);a=zeros(n-1,1);c=zeros(n-1,1);fori=1:n-1b(i,1)=A(i,i);a(i,1)=A(i+1,1);c(i,1)=A(i,i+1);endb(n,1)=A(i,i);fori=2:nb(i,1)=f(i,1)-(a(i-1,1)/b(i-1,1))*c(i-1,1);f(i,1)=f(i,1)-(a(i-1,1)/b(i-1,1))*f(i-1,1);endx(n)=f(n,1)-b(n,1);fori=(n-1):-1:1x(i)=(f(i,1)-c(i,1)*x(i+1))/b(i,1);end求解熱傳導方程的程序function[u,x,y,t]=sjy(a,D,T,u_xy0,u_xyt,Mx,My,N)%a為方程系數(shù)%D為在x軸和y軸上的邊界值,D(1)<=x<=D(2),D(3)<=y<=D(4)%T為時間上限%u_xy0為t=0時的初值%u_xyt為邊界取值函數(shù)/r