使用Python仿真波的叠加

前言：

1.波动知识
波动是一种传递能量的方式，根据物理学的知识，波动具有独立性，叠加性和相干性。独立性是指当若干个波相遇于同一区域时，可以各自保持自己的特性（频率、振幅和振动方向等），按照原来的传播方向继续前进，彼此不受影响。叠加性是指在相遇区域内，介质质点的合位移是各波分别单独传播时在该点引起位移的矢量和。这种叠加性是以独立性为基础的线性叠加。
叠加分为相干叠加和非相干叠加，这取决于各波的特性，以两个沿同一直线的且频率相同简谐波为例，
其中：
中括号括起来的便是干涉项，干涉的结果取决于各波的初始相位，只有具有固定相位差的波才能发生干涉。

2.图像与数组
图像其实就是数组，比如，彩色的图像可以表示为一个三维数组(M,N,C)，其中M和N表示分辨率，C表示颜色通道[R,G,B]，利用python的numpy库可以很方便地实现对数组的访问和修改，基于这一点，想到了使用一个二维数组（灰度图）表示一个空间，图上每一点代表空间的坐标(x,y)，当波在该空间传播时，会引起质点的振动，导致每一点的灰度值发生变化（灰度值相当于质点的位移），同时，利用numpy数组的广播机制可以轻松实现波的叠加。

正文:

需要的库：

import numpy as np
from matplotlib import pyplot as plt
from scipy import sparse

没有的话可以使用pip安装(pip install xxx)。

首先编写描述波的类，初始化函数要覆盖波的特性：振幅、频率、波数、初始相位。我们用1，2，3，4表示波源的位置，1是左上角，2是右上角，3是右下角，4是左下角。我们使用函数sparse.coo_matrix(space)获取空间上每一点的坐标（这里需要用到元组拆包的技巧），获取到坐标后，计算每一点与波源的距离，然后计算波源在该点引起的振动，振动的大小与该点灰度值的大小一一对应，也就是说，颜色越深的点，位移越大。
此外，波还随时间变化，我们定义一个函数，当调用这个函数时，时间会自加若干秒，同时质点空间振动的分布发生改变。

class Wave():
    def __init__(self,A,W,K,φ,corner):
        self.A = A
        self.W = W
        self.K = K
        self.φ= φ
        self.corner = corner
        self.t = 0.0
        self.space = None
        
    def spread(self,space):
        self.space = space.copy()
        coo = sparse.coo_matrix(space)
        if self.corner == 1:
            for (x,y) in zip(coo.col,coo.row):
                r = (x**2+y**2)**(1/2)
                self.space[y][x] = self.A*np.cos(self.W*self.t-self.K*r+self.φ) +  self.A             
        elif self.corner == 2:
            for (x,y) in zip(coo.col,coo.row):
                r = ((space.shape[1]-x)**2+y**2)**(1/2)
                self.space[y][x] = self.A*np.cos(self.W*self.t-self.K*r+self.φ) +  self.A            
        elif self.corner == 3:
            for (x,y) in zip(coo.col,coo.row):
                r = ((space.shape[1]-x)**2+(space.shape[0]-y)**2)**(1/2)
                self.space[y][x] = self.A*np.cos(self.W*self.t-self.K*r+self.φ) +  self.A 
        elif self.corner == 4:
            for (x,y) in zip(coo.col,coo.row):
                r = (x**2+(space.shape[0]-y)**2)**(1/2)
                self.space[y][x] = self.A*np.cos(self.W*self.t-self.K*r+self.φ) +  self.A 
        else:
            raise Exception('Corner must be 1,2,3,4,please check your param')
            
    def move(self):
        self.t += 5
        self.spread(self.space)

定义一个空间（100*100），初始化四个波，输入各自的参数，让各个波在空间传播，并可视化

space = np.ones(shape=(100,100),dtype='float32')

#振幅、频率、波数、初始相位、波源位置
wave_1 = Wave(100,1*np.pi/3,5,0,1) 
wave_2 = Wave(100,2*np.pi/3,5,0,2) 
wave_3 = Wave(100,3*np.pi/3,5,0,3) 
wave_4 = Wave(100,4*np.pi/3,5,0,4)

for wave in [wave_1,wave_2,wave_3,wave_4]:
	wave.spread(space)
        
while True:
	plt.subplot(2,2,1)
	plt.imshow(wave_1.space)

	plt.subplot(2,2,2)
	plt.imshow(wave_2.space)

	plt.subplot(2,2,4)
	plt.imshow(wave_3.space)

	plt.subplot(2,2,3)
	plt.imshow(wave_4.space)

	for wave in [wave_1,wave_2,wave_3,wave_4]:
		wave.move()
	plt.pause(0.1)

效果如图（实际是动态的）中间形状就像水的波纹，两侧像是山脊。:

利用numpy数组的广播机制可以轻松实现波的叠加：

wave_space = wave_1.space+wave_2.space+wave_3.space+wave_4.space
wave_space /= np.max(wave_space)

效果如图（实际是动态的）：可见相干相长的点（颜色较深的点)和相干相消的点(颜色较浅的点)在空间有规律的分布

完整代码如下所示：

import numpy as np
from matplotlib import pyplot as plt
from scipy import sparse

class Wave():
    def __init__(self,A,W,K,φ,corner):
        self.A = A
        self.W = W
        self.K = K
        self.φ= φ
        self.corner = corner
        self.t = 0.0
        self.space = None
             
    def spread(self,space):
        self.space = space.copy()
        coo = sparse.coo_matrix(space)
        if self.corner == 1:
            for (x,y) in zip(coo.col,coo.row):
                r = (x**2+y**2)**(1/2)
                self.space[y][x] = self.A*np.cos(self.W*self.t-self.K*r+self.φ) +  self.A             
        elif self.corner == 2:
            for (x,y) in zip(coo.col,coo.row):
                r = ((space.shape[1]-x)**2+y**2)**(1/2)
                self.space[y][x] = self.A*np.cos(self.W*self.t-self.K*r+self.φ) +  self.A            
        elif self.corner == 3:
            for (x,y) in zip(coo.col,coo.row):
                r = ((space.shape[1]-x)**2+(space.shape[0]-y)**2)**(1/2)
                self.space[y][x] = self.A*np.cos(self.W*self.t-self.K*r+self.φ) +  self.A 
        elif self.corner == 4:
            for (x,y) in zip(coo.col,coo.row):
                r = (x**2+(space.shape[0]-y)**2)**(1/2)
                self.space[y][x] = self.A*np.cos(self.W*self.t-self.K*r+self.φ) +  self.A 
        else:
            raise Exception('Corner must be 1,2,3,4,please check your param')     
            
    def move(self):
        self.t += 5
        self.spread(self.space)
        
if __name__=='__main__':

    space = np.ones(shape=(100,100),dtype='float32')
    
    wave_1 = Wave(100,1*np.pi/3,5,0,1) 
    wave_2 = Wave(100,2*np.pi/3,5,0,2) 
    wave_3 = Wave(100,3*np.pi/3,5,0,3) 
    wave_4 = Wave(100,4*np.pi/3,5,0,4)
    
    for wave in [wave_1,wave_2,wave_3,wave_4]:
        wave.spread(space)
        
    while True:             
        plt.subplot(2,2,1)
        plt.imshow(wave_1.space)
        
        plt.subplot(2,2,2)
        plt.imshow(wave_2.space)
        
        plt.subplot(2,2,4)
        plt.imshow(wave_3.space)
        
        plt.subplot(2,2,3)
        plt.imshow(wave_4.space)
        
        for wave in [wave_1,wave_2,wave_3,wave_4]:
            wave.move()
                    
        plt.pause(0.1)
        
        '''Superposition of waves
        wave_space = wave_1.space+wave_2.space+wave_3.space+wave_4.space
        wave_space /= np.max(wave_space)
        for wave in [wave_1,wave_2,wave_3,wave_4]:
            wave.move()
        plt.imshow(wave_space)
        plt.pause(0.1)
        '''