使用 Python 解决计算物理学问题

引言

拉普拉斯方程是一个简单的二阶偏微分方程。它也是椭圆型偏微分方程中最简单的例子。这个方程在科学，尤其是在物理学中非常重要，因为它描述了电势和引力势的行为，以及热传导。在热力学（热传导）中，我们将拉普拉斯方程称为稳态热方程或热传导方程。

在本文中，我们将使用数值方法而不是解析/微积分方法来求解拉普拉斯方程。当我们说数值方法时，我们指的是“离散化”。离散化是将微分方程的连续形式“转换”为微分方程的离散形式的过程；这也意味着通过离散化，我们可以将微积分问题转化为矩阵代数问题，这是编程所偏爱的。

在这里，我们将使用有限差分法来解决一个简单的热传导问题。我们将使用 Python 编程语言、Numpy（Python 的数值库）和 Matplotlib（用于在 Python 中绘制和可视化数据的库）作为工具。我们还将看到，使用 Python 可以写更少的代码，做更多的事情。

背景

在计算物理学中，我们“总是”使用编程来解决问题，因为计算机程序可以“快速”地进行大规模和复杂的计算。计算物理学可以用下图表示。

如今，有许多编程语言被用来解决许多数值问题，例如 Matlab。但在这里，我们将使用 Python，这是一种“易于学习”的编程语言，当然，它是免费的。它还拥有强大的数值、科学和数据可视化库，如 Numpy、Scipy 和 Matplotlib。Python 还提供并行执行，我们可以在计算机集群上运行它。

回到拉普拉斯方程，我们将在下一节中使用 Python 求解一个简单的二维热传导问题。这里，我假设读者对有限差分法有基本了解，因此我不详细介绍有限差分法的细节、离散化误差、稳定性、一致性、收敛性以及最快/最优迭代算法。我们将在此处跳过许多计算公式的步骤。

我们不进行数值-解析验证来求解问题，而是仅演示如何使用 Python、Numpy 和 Matplotlib 求解问题，当然，还会带有一些简化的计算物理学概念，以便不专门研究计算物理学的普通读者也能理解这里的源代码。

准备工作

为了生成下面的结果，我使用了这个环境

• 操作系统：Linux Ubuntu 14.04 LTS
• Python：Python 2.7
• Numpy：Numpy 1.10.4
• Matplotlib：Matplotlib 1.5.1

如果您正在运行 Ubuntu，您可以使用 pip 安装 Numpy 和 Matplotlib，或者您可以在终端中运行此命令。

$ sudo apt-get install python-numpy

并使用此命令安装 Matplotlib

$ sudo apt-get install python-matplotlib

请注意，Python 已预装在 Ubuntu 14.04 中。要尝试 Python，只需在终端中输入 Python 并按 Enter 键。

您也可以在 Windows 操作系统中使用 Python、Numpy 和 Matplotlib，但我更喜欢使用 Ubuntu。

Using the Code

这是二维笛卡尔坐标系下的拉普拉斯方程（对于热方程）

其中 T 是温度，x 是 x 维度，y 是 y 维度。x 和 y 是笛卡尔坐标系中位置的函数。如果您有兴趣查看上述方程的解析解，可以从这里查看。

这里，我们只需要求解拉普拉斯方程的二维形式。要解决的问题如下所示

我们将要做的，是找出上述二维平板内的稳态温度（这也意味着拉普拉斯方程的解），并给定边界条件（平板边缘的温度）。接下来，我们将离散化平板区域，并将其划分为网格，然后我们使用有限差分法离散化上述拉普拉斯方程。这是平板的离散化区域。

我们设 Δx = Δy = 1 cm，然后作成如下网格

请注意，绿色节点是我们想要知道温度（解）的节点，白色节点是边界条件（已知温度）。这里是上述拉普拉斯方程的离散形式。

在我们的例子中，最终的离散方程如下所示。

现在，我们准备好求解上述方程了。为了求解它，我们使用内部网格（绿色节点）的“猜测值”，这里我们将其设为 30 摄氏度（或者我们可以将其设为 35 度或其他值），因为我们不知道网格内部的值（当然，这些是我们想要知道的值）。然后，我们将迭代方程，直到迭代前的值与迭代后的值之差“足够小”，我们称之为收敛。在迭代过程中，内部网格的温度值会自我调整，它是“自纠正”的，所以当我们设置的猜测值越接近实际解，我们就能越快得到“实际”解。

我们已准备好编写源代码。为了使用 Numpy 库，我们需要在源代码中导入 Numpy，我们还需要导入 Matplotlib.Pyplot 模块来绘制我们的解。因此，第一步是导入必要的模块。

    import numpy as np
    import matplotlib.pyplot as plt

然后，我们将初始变量设置到我们的 Python 源代码中

    # Set maximum iteration
    maxIter = 500

    # Set Dimension and delta
    lenX = lenY = 20 #we set it rectangular
    delta = 1

    # Boundary condition
    Ttop = 100
    Tbottom = 0
    Tleft = 0
    Tright = 0

    # Initial guess of interior grid
    Tguess = 30

接下来，我们将设置“绘图窗口”和 meshgrid。这是代码

    # Set colour interpolation and colour map.
    # You can try set it to 10, or 100 to see the difference
    # You can also try: colourMap = plt.cm.coolwarm
    colorinterpolation = 50
    colourMap = plt.cm.jet

    # Set meshgrid
    X, Y = np.meshgrid(np.arange(0, lenX), np.arange(0, lenY))

np.meshgrid() 为我们创建网格（我们用它来绘制解），第一个参数是 x 维度，第二个参数是 y 维度。我们使用 np.arange() 来安排一个从某个值到某个值的元素值的一维数组，在我们的例子中，是从 0 到 lenX 和从 0 到 lenY。然后我们设置区域：我们定义一个二维数组，定义大小并将数组填充为猜测值，然后我们设置边界条件，请注意上面这里填充数组元素的语法。

    # Set array size and set the interior value with Tguess
    T = np.empty((lenX, lenY))
    T.fill(Tguess)

    # Set Boundary condition
    T[(lenY-1):, :] = Ttop
    T[:1, :] = Tbottom
    T[:, (lenX-1):] = Tright
    T[:, :1] = Tleft

然后我们就可以将最终方程应用到下面的 Python 代码中。我们使用 for 循环来迭代方程。

    # Iteration (We assume that the iteration is convergence in maxIter = 500)
    print("Please wait for a moment")
    for iteration in range(0, maxIter):
    	for i in range(1, lenX-1, delta):
	    	for j in range(1, lenY-1, delta):
		    	T[i, j] = 0.25 * (T[i+1][j] + T[i-1][j] + T[i][j+1] + T[i][j-1])

    print("Iteration finished")

您应该注意上面的代码缩进，Python 不使用括号，而是使用空格或缩进。好了，主要逻辑已经完成了。接下来，我们编写代码来绘制解，使用 Matplotlib。

    # Configure the contour 
    plt.title("Contour of Temperature")
    plt.contourf(X, Y, T, colorinterpolation, cmap=colourMap)

    # Set Colorbar
    plt.colorbar()

    # Show the result in the plot window
    plt.show()

    print("")

就这样，这是 **完整代码**。

    # Simple Numerical Laplace Equation Solution using Finite Difference Method
    import numpy as np
    import matplotlib.pyplot as plt

    # Set maximum iteration
    maxIter = 500

    # Set Dimension and delta
    lenX = lenY = 20 #we set it rectangular
    delta = 1

    # Boundary condition
    Ttop = 100
    Tbottom = 0
    Tleft = 0
    Tright = 30

    # Initial guess of interior grid
    Tguess = 30

    # Set colour interpolation and colour map
    colorinterpolation = 50
    colourMap = plt.cm.jet #you can try: colourMap = plt.cm.coolwarm

    # Set meshgrid
    X, Y = np.meshgrid(np.arange(0, lenX), np.arange(0, lenY))

    # Set array size and set the interior value with Tguess
    T = np.empty((lenX, lenY))
    T.fill(Tguess)

    # Set Boundary condition
    T[(lenY-1):, :] = Ttop
    T[:1, :] = Tbottom
    T[:, (lenX-1):] = Tright
    T[:, :1] = Tleft

    # Iteration (We assume that the iteration is convergence in maxIter = 500)
    print("Please wait for a moment")
    for iteration in range(0, maxIter):
	    for i in range(1, lenX-1, delta):
		    for j in range(1, lenY-1, delta):
			    T[i, j] = 0.25 * (T[i+1][j] + T[i-1][j] + T[i][j+1] + T[i][j-1])

    print("Iteration finished")

    # Configure the contour 
    plt.title("Contour of Temperature")
    plt.contourf(X, Y, T, colorinterpolation, cmap=colourMap)

    # Set Colorbar
    plt.colorbar()

    # Show the result in the plot window
    plt.show()

    print("")

它相当简短，*不是吗*？好的，您可以复制粘贴并保存源代码，将其命名为 findif.py。要执行 Python 源代码，请打开您的终端，然后进入您存放源代码的目录，键入

$ python findif.py

然后按 **Enter**。然后将出现绘图窗口。

您可以尝试更改边界条件的值，例如，将右侧边缘温度的值更改为 30 摄氏度（Tright = 30），结果将如下所示

关注点

Python 是一种“易于学习”的动态类型编程语言，它为计算物理学或其他科学学科提供了（开源的）强大库。由于 Python 是一种解释型语言，与 C 或 C++ 等编译型语言相比，它的速度较慢，但同样，它易于学习。我们也可以用更少的代码做更多的事情，这样我们就不会在编程上挣扎，而是可以专注于我们想要解决的问题。

在计算物理学中，使用 Numpy 和 Scipy（Python 的数值和科学库），我们可以解决许多复杂问题，因为它提供了矩阵求解器（特征值和特征向量求解器）、线性代数运算，以及信号处理、傅里叶变换、统计、优化等。

除了在计算物理学中的应用，Python 也被用于机器学习，甚至 Google 的 TensorFlow 也使用 Python。

历史

• 2016 年 3 月 21 日：初始版本