波动现象在生活中非常常见，比如你随便扔一颗石子到平静的湖面上，一圈圈的波纹图案就会出现。波动现象的控制方程为波动方程，下面不要眨眼，请欣赏美丽的波纹：正方形域内波反射图案矩形区域波反射图案三角形区域(一条边为无反射边界)波反射图案

只要我们求解出波动方程我们就可以得到上面美丽的图案，那么什么是波动方程呢，二维的波动方程它大概长这个样子：

这个方程可以描述薄膜的振动，可以理解为波在一个绷紧的弹性膜中的传播。学过数理方程的小伙伴可能可以手算出这个方程在矩形区域里面的解析解，但是我们今天不求解析解，我们利用有限差分法数值求解波动方程，这里采用的语言是Python，Python这个语言很友好，可以说通俗易懂。既然要用有限差分法对方程求解，第一步就是要把上述的偏微分方程离散为差分方程(这里采用的是三点差分格式)：

上述差分方程U的上标为时间标，下标为空间标。将上述的方程做一些化简：

其中：

这样k+1时间步的结果就可以通过k和k-1时间步的结果推进得到，这里会有一个很关键的地方，比如当前时间步是k=0，k-1=-1这个是不存在的，我们没有负的时间步，那该怎么办呢？这里利用我们的初始条件我们可以外插得到我们-1时间步的结果：

v0为初始状态下弹性膜中每一点的速度。为了保证求解的收敛性，三点差分格式要求步长比小于等于1：

这里有一个注意的地方，除了初始条件(定义t=0时弹性膜上每一点位移与速度)之外我们还需要边界条件，用来定义弹性膜边界在以后的每个时刻的状态，如下图红线圈出的点是位于边界处的点，边界上的的值是不需要求解的，可以用边界条件计算得出。

下面我们求解一个2X2方形区域内的波动方程，边界条件与初始条件如下：

下面我要放出代码了，说实话如果没有上面的推导，直接去看代码很多人估计看不懂，如果真的对偏微分方程的数值解有兴趣可以看看《偏微分方程的数值解》李荣华版，如果大家想要电子书的话可以私信我。

import numpy as np

import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

from matplotlib import cm

from PIL import Image

class WaveEquationFD:

def __init__(self, N, D, Mx, My):

self.N = N

self.D = D

self.Mx = Mx

self.My = My

self.tend = 6

self.xmin = 0

self.xmax = 2

self.ymin = 0

self.ymax = 2

self.initialization()

self.eqnApprox()

def initialization(self):

self.dx = (self.xmax - self.xmin)/self.Mx

self.dy = (self.ymax - self.ymin)/self.My

self.x = np.arange(self.xmin, self.xmax+self.dx, self.dx)

self.y = np.arange(self.ymin, self.ymax+self.dy, self.dy)

#----- Initial condition -----#

self.u0 = lambda r, s: 0.2*np.exp(-((r-(self.xmax + self.xmin)/2)**2+(s-(self.ymax + self.ymin)/2)**2)/0.01)

#----- Initial velocity -----#

self.v0 = lambda a, b: 0

#----- Boundary conditions -----#

self.bxyt = lambda left, right, time: 0

self.dt = (self.tend - 0)/self.N

self.t = np.arange(0, self.tend+self.dt/2, self.dt)

# Assertion for the condition of r < 1, for stability

r = 4*self.D*self.dt**2/(self.dx**2+self.dy**2);

assert r < 1, "r is bigger than 1!"

def eqnApprox(self):

#----- Approximation equation properties -----#

self.rx = self.D*self.dt**2/self.dx**2

self.ry = self.D*self.dt**2/self.dy**2

self.rxy1 = 1 - self.rx - self.ry

self.rxy2 = self.rxy1*2

#----- Initialization matrix u for solution -----#

self.u = np.zeros((self.Mx+1, self.My+1))

self.ut = np.zeros((self.Mx+1, self.My+1))

self.u_1 = self.u.copy()

#----- Fills initial condition and initial velocity -----#

for j in range(1, self.Mx):

for i in range(1, self.My):

self.u[i,j] = self.u0(self.x[i], self.y[j])

self.ut[i,j] = self.v0(self.x[i], self.y[j])

def solve_and_animate(self):

u_2 = np.zeros((self.Mx+1, self.My+1))

xx, yy = np.meshgrid(self.x, self.y)

fig = plt.figure()

ax = fig.add_subplot(111, projection='3d')

wframe = None

count=1

k = 0

nsteps = self.N

while k < nsteps:

if wframe:

ax.collections.remove(wframe)

self.t = k*self.dt

#----- Fills in boundary condition along y-axis (vertical, columns 0 and Mx) -----#

for i in range(self.My+1):

self.u[i, 0] = self.bxyt(self.x[0], self.y[i], self.t)

self.u[i, self.Mx] = self.bxyt(self.x[self.Mx], self.y[i], self.t)

for j in range(self.Mx+1):

self.u[0, j] = self.bxyt(self.x[j], self.y[0], self.t)

self.u[self.My, j] = self.bxyt(self.x[j], self.y[self.My], self.t)

if k == 0:

for j in range(1, self.My):

for i in range(1, self.Mx):

self.u[i,j] = 0.5*(self.rx*(self.u_1[i-1,j] + self.u_1[i+1,j])) \

+ 0.5*(self.ry*(self.u_1[i,j-1] + self.u_1[i,j+1])) \

+ self.rxy1*self.u[i,j] + self.dt*self.ut[i,j]

else:

for j in range(1, self.My):

for i in range(1, self.Mx):

self.u[i,j] = self.rx*(self.u_1[i-1,j] + self.u_1[i+1,j]) \

+ self.ry*(self.u_1[i,j-1] + self.u_1[i,j+1]) \

+ self.rxy2*self.u[i,j] - u_2[i,j]

u_2 = self.u_1.copy()

self.u_1 = self.u.copy()

wframe = ax.plot_surface(xx, yy, self.u, cmap='rainbow', linewidth=2,

antialiased=False)

ax.set_xlim3d(0, 2.0)

ax.set_ylim3d(0, 2.0)

ax.set_zlim3d(-1.5, 1.5)

ax.set_xticks([0, 0.5, 1.0, 1.5, 2.0])

ax.set_yticks([0, 0.5, 1.0, 1.5, 2.0])

ax.set_xlabel("x")

ax.set_ylabel("y")

ax.set_zlabel("U")

plt.savefig(str(count))

plt.pause(0.01)

k += 0.5

count+=1

def main():

simulator = WaveEquationFD(200, 0.1, 100, 100)

simulator.solve_and_animate()

plt.show()

if __name__ == "__main__":

main()

#N = 200

#D = 0.25

#Mx = 50

#My = 50

imgs = []

下面就是见证奇迹的时刻：

差分方程的推导很简单，但是如果没有推导明白的很难理解上面的程序，如果有疑问的小伙伴可以私信我，想要《偏微分方程数值解》的小伙伴也可以私信我。