利用Python求解拉普拉斯橢圓方程中熱傳導問題

丁小勇，鄭濱紅

（上饒幼兒師范高等專科學校，江西 上饒 334000）

拉普拉斯方程[1]是以Pierre-imon Laplace命名的二階偏微分方程。拉普拉斯方程是橢圓偏微分方程的最簡單例子。求解拉普拉斯方程是電磁學、天文學和流體力學等領域經常遇到的一類重要的數學問題，因為該方程以勢函數的形式描寫了電場、引力場和流場等物理對象的性質，求解拉普拉斯方程就成了解決所有這些問題的關鍵所在。在本文中，我們將用數值求解[2]代替積分求解，來求解拉普拉斯方程中熱傳導問題。

當我們提到數值方法時，就意味著離散化。離散化就是將連續形式的微分方程轉化為離散形式。同時它也意味著我們將積分問題轉化為線性代數問題，這給我們使用Python解決問題提供了思路。

在熱傳導中，我們如何在給定邊界溫度的情況下，利用拉普拉斯方程求出二維平面內每個點的穩定溫度，即拉普拉斯方程的解呢？本文中，我們借助計算機語言Python以及數值求解和有限差分法來求解熱傳導問題。

1 微分方程數值求解——有限差分法

1.1 有限差分法概述

有限差分法[3]（Finite Difference Method，FDM）是一種求解微分方程數值解的近似方法，其主要原理是對微分方程中的微分項進行直接差分近似，從而將微分方程轉化為代數方程組求解。有限差分法其面對的對象是微分方程，包括常微分方程和偏微分方程。此外，有限差分法需要對微分進行近似，這里的近似采取的是離散近似，使用某一點周圍點的函數值近似表示該點的微分。

1.2 有限差分法的基本原理

以一維函數f(x)為例。現在已知離散點對應，假設任何一個函數都展開為多項式（泰勒展開）

然后可以組成一個矩陣方程：

這里要求出矢量a，需要做一個求逆操作Y=X-1，并且

然后可以發現在x=0處0～n階導數已經可以求出。如果想求出其他點的n階導數，我們展開該函數即可。

2 用Python解拉普拉斯橢圓方程的優勢

Python 是一種簡單易用的動態計算機編程語言，在數學領域、計算物理領域擁有非常廣泛的應用，因為它有很多有用的庫。我們可以使用較少的Python代碼來完成更多的工作。這樣一來，我們不必糾結于如何編程，只需將精力放在想要解決的問題上面。

在解拉普拉斯橢圓方程中，使用NumPy（numerical library for Python）[4]庫和Scipy庫（numeric and scientific library for Python）可以幫我們解決很多復雜的問題，因為這些庫為我們提供了矩陣運算、線性代數操作與信號處理、傅里葉變化、統計分析優化等現成的函數。除了在數學領域和計算物理領域，Python在機器學習領域也有大量的應用。很多機器學習框架都是基于Python開發的，包括谷歌的 Tensorflow。這說明Python不僅在學術界具有廣泛的應用，在工業界也擁有很強的實用性。

3 拉普拉斯橢圓方程中熱傳導問題

拉普拉斯方程是一個二階偏微分方程,其描述了空間中無源場的分布情況，因此，我們在二維笛卡兒坐標系中設定拉普拉斯方程（用于熱方程）如圖1所示。其中T是溫度，x是x維度，y是y維度。x和y是笛卡兒坐標系中位置的函數。

圖1 二維笛卡兒坐標系

我們要做的是在給定的邊界條件（平板邊緣的溫度）下，找到上面二維平板內部的穩態溫度（拉普拉斯方程的解）。

接下來，我們將對平面的區域進行離散化，并將其劃分為網格，然后用有限差分法對上面的拉普拉斯方程進行離散。圖2展示了離散化后的平面區域，設Δx=Δy=1 cm，然后制作如圖2所示的網格圖。

圖2 平面離散區域網格圖

在拉普拉斯方程的平面離散形式中，我們使用內部網格（綠色節點）的“估計值”，這里我們將其設置為30℃（或者我們可以將其設置成35℃或其他值），因為我們不知道網格內部的值（當然，這些是我們想知道的值）。然后我們將用迭代方程進行運算，直到迭代前的值與迭代后的數值之間的差足夠小，我們稱之為收斂。在迭代過程中，內部網格中的溫度值會自行調整，這是自校正的，所以當我們將估計值設置為更接近其實際解時，我們得到實際解的速度就更快。值得注意的是，綠色節點（內部網格溫度見圖3）是我們想知道那里的溫度（解決方案）的節點，而白色節點是邊界條件（已知溫度）。

圖3 內部網格溫度圖

4 Python下拉普拉斯橢圓方程熱傳導問題解決方案

本文使用的編程語言是Python，應用的庫有：Numpy(numerical library for Python)和 Matplotlib(library for plotting and visualizing data using Python)。下面將看到，本文可以使用Python 代碼來實現要求解的值。為了解決這個問題，本文要使用Numpy庫。需要在源代碼中導入Numpy，還需要導入Matplolib.Pyplot模塊來繪制我們的解決方案。

4.1 導入Numpy庫和Matplolib.Pyplot模塊并初始化變量

在編寫代碼之前，第一步是導入必要的模塊。然后，我們將初始變量設置到Python源代碼中。

4.2 設置“繪圖窗口”和網格

np.meshgrid（）為我們創建的網格（我們使用它來繪制解決方案），第一個參數用于x維，第二個參數用于y維。我們使用np.arange（）來排列元素值從某個值開始到某個值的一維數組。在我們的例子中，它從0到lenX，從0到lenY。然后我們設置區域：定義一個二維數組，定義其形狀，并使用猜測值對其進行初始化，然后設置邊界條件，使用前面給出的邊界條件進行初始化。

4.3 使用for循環實現迭代方程

將我們的最終公式應用到下面的Python代碼中。我們使用for循環迭代方程。

我們應該注意到上面代碼的縮進，Python代碼使用的是空格或者縮進，而不使用大括號。至此我們的主程序部分已經完成了，接下來，我們使用 Matplotlib來畫出最后的結果。

4.4 使用Matplotlib編寫代碼來繪制解決方案

接下來，我們使用Matplotlib編寫代碼來繪制解決方案。

以上代碼運行結果如圖4所示。

圖4 邊緣溫度值35℃結果圖

嘗試更改邊界條件的值，例如，將右邊緣溫度的值更改為30攝氏度（Tright=30），運行結果如圖5所示。

圖5 邊緣溫度值30℃結果圖

5 結束語

本文利用Python中Numpy(numerical library for Python) 庫和 Matplotlib(library for plotting and visualizing data using Python) 庫直觀地呈現了在給定邊界溫度的情況下，利用拉普拉斯方程求出二維平面內每點的穩定溫度。