熱導方程的Matlab數值解方法

這是一個很久很久以前的一個故事，久到能夠讓人忘記原來這這些方程是如此的貼近自己的學習。你學或者不學，它都在這裡，不難也不簡單。過冷水今天就和大家分享一下一維熱傳導方程特別案例的具體求解方法。

熱傳導是一個很常見的現象。當物體內部的溫度分布不均勻時，熱量就會從溫度較高的地方流動，這個過程中，溫度是空間和時間的函數。熱傳導方程就是溫度所滿足的偏微分方程，它的解給出任意時刻物體內的溫度分布。

為了建立熱導方程，我們首先介紹熱導系統置於x軸，考查系統在任意x處的橫截面上的一個單位面積，設熱流沿x軸方向傳遞，x處的溫度為u(x),溫度梯度為du(x)/dx。傅立葉指出：在單位面積內流經該單位面積的熱量q與該處的溫度梯度成正比即：

k：熱導率,負號表示與溫度梯度方向相反。現在假設這個一維熱傳導系統的長度為l，橫截面面積為s杆的兩個端點處於x=0和x=l.假定杆在初始t=0時刻溫度分布為Φ(x)，在隨後的時間（t>0）,熱量在杆中流動。現在我們要確定在任意t時刻，杆中任意位置x(0<x<l)的溫度u(x,t)。

我們考查系統在x位置的一段∆x,根據傅立葉定律在∆t時間內從∆x前端流入的熱量為：

    另一方面，在該時間內從後端流出的熱量為：

在沒有其他熱源的情況下，體積元S∆x吸收的熱量使之溫度升高。而溫度升高的描述則是基於熱體比熱c的定義：

其中，m是物體的質量，c表示單位質量的物體溫度升高1K所需要的熱量。這樣體積元S∆x吸收的熱量為：

其中，ρ=m/(S∆x)是系統質量體密度。熱量守恆要求：

則：

對其進行進一步變化可得：

這就是所謂的一維系統的熱傳導方程。我們對熱傳導方程進行一個簡單的分析，若時間的微商項du(x)/dt=0這是穩態過程。則d2u(x)/dx2=0則:

有熱源的熱傳導方程為：

我們來看一個比較簡單形式的求解方法。

該條件下的熱導方程求解，採用兩種不同的形式分離變量法和差分法。我們先來看分離變量法：

則：

由邊界條件u(0,t)=0,u(l,t)=0可得X(0)=X(l)=0,求邊值問題:

解：

這裡需要解釋一下X、、(x)+λX(x)=0微分方程根據λ<0，λ=0，λ>0；表示成不同函數類型，除λ>0能夠得到符合邊界條件的函數外，其它都不符合邊界條件。

現在考慮：

將特徵值λ帶入方程的：

通解為：

於是：

再利用初值條件：u(x,0)=φ(x)可得：

於是最終解就是給出來：

我們看一道有具體條件的題：

再利用初值條件：u(x,0)=φ(x)可得：

最終結果有沒有覺得神秘複雜的熱導方程好像也不是那麼難計算，就是一個累計加和的形式，很簡單。讀者需要注意的是熱導方程的形式是和邊界條件有關係的，不同的邊界條件最終的形式差別是很大的，我們來看一下代碼：

x=0:0.1*pi:pi;y=0:0.04:1;[x,t]=meshgrid(x,y);s=0;m=length(j);%matlab可計算的最大數   for i=1:m    s=s+(200*(1-(-1)^i))/(i*pi)*(sin(i*x).*exp(-i^2*t));end;surf(x,t,s);xlabel('x'),ylabel('t'),zlabel('T');title(' 分離變量法（無窮）');axis([0 pi 0 1 0 100])
熱導方程的數值解代碼出乎意料的簡潔。我們再來看一下另外一種求解方法：有限差分方法。
有限差分：將求解域劃分為差分網格，用有限個網格節點代替連續的求解域。有限差分方法以泰勒級數展開等方法，把控制方程中的導數用網格節點上函數值的差商代替進行離散，從而建立以網格節點上的值為未知數的代數方程組.
離散化：
其代碼實現為：
%有限差分法：u=zeros(10,25);%橫坐標為x,縱坐標為t;s=(1/25)/(pi/10)^2;fprintf('穩定性係數S為:\n');disp(s);for i=2:9;    u(i,1)=100;end;for j=1:25;    u(1,j)=0;    u(10,j)=0;end;for j=1:24;    for i=2:9;        u(i,j+1)=s*u(i+1,j)+(1-2*s)*u(i,j)+s*u(i-1,j);    endenddisp(u);[x,t]=meshgrid(1:25,1:10);surf(x,t,u);xlabel('t');ylabel('x');zlabel('T');title('有限差分法解');
這就是過冷書想要和大家分享的關於一維熱傳導方程求解的方法，數值解的代碼過程很簡單，主要是數學問題，第一種方法用到了分離變量的思想使得溫度變得簡單。第二種方法就是用具體值來近似表示熱導方程。使得問題變得簡單。看完之後才有豁然開朗的感覺，數學也沒有想像中的那麼難。限於篇幅一部分人所關注的二維熱傳導方程敬請起來後期會和大家分享二維熱導方程案例，具體實現代碼。
公眾號中長期支持投稿，若是對相關方程的計算有獨到體會可以進行投稿，在公眾號回復投稿，聯繫小編即可。
要求：有必要的算法介紹、公式推導(公式需要能編輯)、完整可執行的代碼、自行選擇一個稍微複雜的微分方程作為求解示例，提交文稿一律用Word，並提供相關m文件和圖片，以上題目網上基本上都有原始碼，參考後必須註明相關參考連結，切勿直接抄襲，需要有自己的理解和分析
如需轉載，請在公眾號中回復「轉載」獲取授權，未經授權擅自搬運抄襲的，必將追究其責任！
往期回顧>>>>>>
從泰勒級數說傅立葉級數
你所不知道的Monte Carlo形式
基於Hough變換原理實現圖像直線檢測【附原始碼】
數值計算——MATLAB數值積分原理詳講