使用 Python 求解拉普拉斯方程

02-02

問題描述

我們求解的問題是：在給定邊界溫度的情況下，求出二維平面內每點的穩定溫度（即為拉普拉斯方程的解）。

本文follow實例使用 Python 解決計算物理問題

原文

Using the code

This is the Laplace equation in 2-D cartesian coordinates (for heat equation)

Where T is temperature, x is x-dimension, and y is y-dimension. x and y are function of position in cartesian coordinates. If you are interested to see the analitical solution of the equation above, you can look it here.

Here we only need to solve 2-D form of the Laplace equation. The problem to solve is shown below.

What we will do is to find the steady state temperature inside the 2-D plat (which also means the solution of Laplace equation) above with the given boundary conditions (temperature of the edge of the plat). Next, we will discretize the region of the plat, and devide it into meshgrid, and then we discretize the Laplace equation above with finite difference method. This is the discretized region of the plat.

We set Δx = Δy = 1 cm, and then make the grid as shown below.

Note that the green nodes are the nodes that we want to know the temperature there (the solution), and the white nodes are the boundary conditions (known temperature). Here is the discrete form of Laplace Equation above.

In our case, the final discrete equation is shown below.

Now, we are ready to solve the equation above. To solve this, we use "guess value" of interior grid (green nodes), here we set it 30 degree Celsius (or we can set it 35 or other value), because we dont know the value inside the grid (of course, those are the values that we want to know). Then we will iterate the equation until the difference between value before iteration and the value until iteration is "small enough", we call it convergence. In process of iterating, the temperature value in interior grid will adjust itself, its "selfcorrecting", so when we set guess value closer to its actual solution, the faster we get the "actual" solution.

We are ready for the source code. In order to use Numpy library, we need to import Numpy in our source code, and we also need to import Matplolib.Pyplot module to plot our solution. So the first step is import necessary modules.

Hide Copy Code

import numpy as npnimport matplotlib.pyplot as pltn

and then, we set the initial variables into our Python source code

Hide Copy Code

# Set maximum iterationnmaxIter = 500nn# Set Dimension and deltanlenX = lenY = 20 #we set it rectangularndelta = 1nn# Boundary conditionnTtop = 100nTbottom = 0nTleft = 0nTright = 0nn# Initial guess of interior gridnTguess = 30n

What we will do next is to set the "plot window" and meshgrid. Here is the code.

Hide Copy Code

# Set colour interpolation and colour map.n# You can try set it to 10, or 100 to see the differencen# You can also try: colourMap = plt.cm.coolwarmncolorinterpolation = 50ncolourMap = plt.cm.jetnn# Set meshgridnX, Y = np.meshgrid(np.arange(0, lenX), np.arange(0, lenY))n

np.meshgrid() creates the mesh grid for us (we use this to plot the solution), the first parameter is for the x-dimension, and the second parameter is for the y-dimension. We use np.arange() to arange 1-D array with element value starts from some value to some value, in our case its from 0 to lenX and from 0 to lenY. Then we set the region: we define 2-D array, define the size and fill the array with guess value, then we set the boundary condition, look at the syntax of filling the array element for boundary condition above here.

Hide Copy Code

# Set array size and set the interior value with TguessnT = np.empty((lenX, lenY))nT.fill(Tguess)nn# Set Boundary conditionnT[(lenY-1):, :] = TtopnT[:1, :] = TbottomnT[:, (lenX-1):] = TrightnT[:, :1] = Tleftn

Then we are ready to apply our final equation into Python code below. We iterate the equation using for loop.

Hide Copy Code

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

You should aware of the indentation of the code above, Python does not use bracket but it uses white space or indentation. Well, the main logic is finished. Next, we write code to plot the solution, using Matplotlib.

Hide Copy Code

# Configure the contournplt.title("Contour of Temperature")nplt.contourf(X, Y, T, colorinterpolation, cmap=colourMap)nn# Set Colorbarnplt.colorbar()nn# Show the result in the plot windownplt.show()nnprint("")n

Thats all, This is the complete code.

Hide Shrink

Copy Code

# Simple Numerical Laplace Equation Solution using Finite Difference Methodnimport numpy as npnimport matplotlib.pyplot as pltnn# Set maximum iterationnmaxIter = 500nn# Set Dimension and deltanlenX = lenY = 20 #we set it rectangularndelta = 1nn# Boundary conditionnTtop = 100nTbottom = 0nTleft = 0nTright = 30nn# Initial guess of interior gridnTguess = 30nn# Set colour interpolation and colour mapncolorinterpolation = 50ncolourMap = plt.cm.jet #you can try: colourMap = plt.cm.coolwarmnn# Set meshgridnX, Y = np.meshgrid(np.arange(0, lenX), np.arange(0, lenY))nn# Set array size and set the interior value with TguessnT = np.empty((lenX, lenY))nT.fill(Tguess)nn# Set Boundary conditionnT[(lenY-1):, :] = TtopnT[:1, :] = TbottomnT[:, (lenX-1):] = TrightnT[:, :1] = Tleftnn# Iteration (We assume that the iteration is convergence in maxIter = 500)nprint("Please wait for a moment")nfor iteration in range(0, maxIter):n for i in range(1, lenX-1, delta):n for j in range(1, lenY-1, delta):n T[i, j] = 0.25 * (T[i+1][j] + T[i-1][j] + T[i][j+1] + T[i][j-1])nnprint("Iteration finished")nn# Configure the contournplt.title("Contour of Temperature")nplt.contourf(X, Y, T, colorinterpolation, cmap=colourMap)nn# Set Colorbarnplt.colorbar()nn# Show the result in the plot windownplt.show()nnprint("")n

Its pretty short huh? Okay, you can copy-paste and save the source code, name it findif.py. To execute the Python source code, open your Terminal, and go to the directory where you locate the source code, type

$ python findif.py

and press enter. Then the plot window will appear.

You can try to change the boundary conditions value, for example you change the value of right edge temperature to 30 degree Celcius (Tright = 30), then the result will look like this

代碼如下

import matplotlib.pyplot as plt #繪圖用的模塊 nfrom mpl_toolkits.mplot3d import Axes3D #繪製3D坐標的函數 nimport numpy as np nnplt.figure(num=1,dpi=100,facecolor=#D3D3D3,edgecolor=r,figsize=(12,9))nn#設置最大迭代次數nmaxIter=5000nn#設置維度和變化率nlenX=lenY=10ndelta=1nn#邊界條件nTtop=100nTbottom=0nTright=0nTleft=0nn#內部網格的猜測溫度nTguess=30nn#設置格點nncolorinterpolation=50nX,Y=np.meshgrid(np.arange(0,lenX),np.arange(0,lenY))nn#設置數組尺寸以及內部溫度TguessnT=np.empty((lenX,lenY))nT.fill(Tguess)nn#設置邊界條件nT[(lenY-1):,:]=TtopnT[:1,:]=TbottomnT[:,(lenX-1):]=TrightnT[:,:1]=Tleftn nn#迭代循環nprint(Please wait for a time)nfor iteration in range(0,maxIter):n for i in range(1,lenX-1,delta):n for j in range(1,lenY-1,delta):n T[i,j]=0.25*(T[i+1][j]+T[i-i][j]+T[i][j+1]+T[i][j+1])n n nprint(Iteration finished)nnplt.title(熱傳導)nplt.contourf(X,Y,T,colorinterpolation,cmap=hot)nplt.colorbar()nnnplt.show()n

三維拉普拉斯方程的求解

import matplotlib.pyplot as plt #繪圖用的模塊 nfrom mpl_toolkits.mplot3d import Axes3D #繪製3D坐標的函數 nimport numpy as np nnfig = plt.figure()nax = Axes3D(fig)nn#設置最大迭代次數nmaxIter=1000nnn#設置維度和變化率nlenX=lenY=lenZ=10ndelta=1nn#邊界條件nTtop=500nTbottom=0nTright=100nTleft=0nTfront=300nTbehind=0n#內部網格的猜測溫度nTguess=30nn#設置格點nncolorinterpolation=50nX,Y,Z=np.meshgrid(np.arange(0,lenX),np.arange(0,lenY),np.arange(0,lenZ))nn#設置數組尺寸以及內部溫度TguessnT=np.empty((lenX,lenY,lenZ))nT.fill(Tguess)nn#設置邊界條件nT[(lenY-1):,:,:]=TtopnT[:1,:,:]=TbottomnT[:,(lenX-1):,:]=TrightnT[:,:1,:]=TleftnT[:,:,:1]=TfrontnT[:,:,-2:]=Tbehind nn#迭代循環nprint(Please wait for a time)nfor iteration in range(0,maxIter):n n for i in range(1,lenX-1,delta):n for j in range(1,lenY-1,delta):n for k in range(1,lenZ-1,delta):n T[i,j]=1/6*(T[i+1][j][k]+T[i-i][j][k]+T[i][j+1][k]+T[i][j-1][k]+T[i][j][k+1]+T[i][j][k-1])n n