利用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) 库直观地呈现了在给定边界温度的情况下，利用拉普拉斯方程求出二维平面内每点的稳定温度。