波动方程矩阵运算代码的有效性

Question

我写了一个波叠加程序，它与多个波源的波方程重叠，然后给出一个包含所有相长和相消干扰的单波，并绘制叠加的强度图。此代码有效，但效率非常低。

• 如果我把这里的点作为1000x1000的网格，整个程序需要一段时间才能运行。

• 有没有什么方法可以使用以下一个或全部(函数、lambda函数、可映射、直接或类似地定义2D numpy数组(使此代码更加高效和干净。

• 如果是，有没有办法衡量运行操作所需的时间。这不是家庭作业，我正试图为我的光学研究建立一些自己的东西。非常感谢您提前提供的帮助，我真的很感激

  import numpy as np
  from matplotlib import pyplot as plt
  from sklearn import mixture
  import matplotlib as mpl
  from mpl_toolkits.mplot3d import Axes3D
  from matplotlib import cm
  import scipy
  import scipy.ndimage
  import scipy.misc 
  xmin, xmax = 0,25
  ymin, ymax = -12.500,12.500
  xpoints, ypoints = 500,500
  points,amp,distance,wave=[],[],[],[]
  numsource=11
  x=np.linspace(xmin, xmax, xpoints)
  y=np.linspace(ymin, ymax, ypoints)
  xx,yy=np.meshgrid(x,y, sparse=False)
  pointer = np.concatenate([xx.reshape(-1,1),yy.reshape(-1, 1)], axis=-1)
  counter=len(pointer)
  A=[0]*counter #for storing amplitudes later
  # Arrays of point source locations
  source=tuple([0,(((numsource-1)/2)-i)*2.5] for i in range(0,numsource-1))
  # Arrays of Subtraction of Coordinates of Sources from Point sources (For Distance from source)   
  points=[(pointer-source[p]) for p in range(0,numsource-1)]
  # distance of each point in map from its source (Sqrt of Coordinate difference sum)
  distance=[(points[i][:,0]**2 + points[i][:,1]**2)**0.5 for i in range(0,numsource-1)]
  # Amplitudes of each wave defined arbitrarily
  amp= np.array([4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4])
  k=20
  # wave equation for each wave based on defined amplitude and distance from source
  wave=([amp[i] * np.sin (k*distance[i]) for i in range(0,numsource-1)])
  #superposition
  for i in range(0,numsource-1):
  A=A+wave[i]
  A=np.asarray(A)
  print(A)
  intensity = A**2
  #constructive, destructive superposition plot
  plt.figure(figsize=(10,10))
  plt.xlim(xmin,xmax)
  plt.ylim(ymin,ymax)
  plt.scatter(pointer[:,0], pointer[:,1],c=intensity, cmap='viridis')
  plt.colorbar()
  

Answer 1

我设法将机器上的计算时间从166毫秒减少到146毫秒。

您希望去掉A初始化为25000x1列表然后转换为数组的部分。您可以直接将其初始化为数组。我禁用的部分将作为注释留在您的代码中。

#import matplotlib as mpl
from matplotlib import pyplot as plt
#from mpl_toolkits.mplot3d import Axes3D
#from matplotlib import cm
import numpy as np
#import scipy
#import scipy.ndimage
#import scipy.misc
#from sklearn import mixture
import time
class time_this_scope():
"A handy context manager to measure timings."""
def __init__(self):
self.t0 = None
self.dt = None
def __enter__(self):
self.t0 = time.perf_counter()
return self
def __exit__(self, *args):
self.dt = time.perf_counter() - self.t0
print(f"This scope took {int(self.dt*1000)} ms.")

def main():
xmin, xmax = 0, 25
ymin, ymax = -12.500, 12.500
xpoints, ypoints = 500, 500
points, amp, distance, wave = [], [], [], []
numsource = 11
x = np.linspace(xmin, xmax, xpoints)
y = np.linspace(ymin, ymax, ypoints)
xx, yy = np.meshgrid(x, y, sparse=False)
pointer = np.concatenate([xx.reshape(-1, 1), yy.reshape(-1, 1)], axis=-1)

#    counter = len(pointer)
#    A = [0]*counter #for storing amplitudes later
A = np.zeros(len(pointer))
# Arrays of point source locations
source = tuple((0, (((numsource-1)/2)-i)*2.5) for i in range(0, numsource-1))
# Arrays of Subtraction of Coordinates of Sources from Point sources (For Distance from source)
points = [(pointer-source[i]) for i in range(0, numsource-1)]
# distance of each point in map from its source (Sqrt of Coordinate difference sum)
distance = [(points[i][:,0]**2 + points[i][:,1]**2)**0.5 for i in range(0, numsource-1)]
# Amplitudes of each wave defined arbitrarily
#    amp = np.array([4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4])
amp = np.array([4]*numsource)
k = 20
# wave equation for each wave based on defined amplitude and distance from source
wave = [amp[i] * np.sin(k*distance[i]) for i in range(0, numsource-1)]
#superposition
for i in range(0, numsource-1):
A = A + wave[i]
#        A = np.asarray(A)

print(A)
intensity = A**2
#constructive, destructive superposition plot
#    plt.figure(figsize=(10, 10))
#    plt.xlim(xmin, xmax)
#    plt.ylim(ymin, ymax)
#    plt.scatter(pointer[:,0], pointer[:,1], c=intensity, cmap='viridis')
#    plt.colorbar()
# Timing.
total_time = 0
N = 20
for _ in range(N):
with time_this_scope() as timer:
main()
total_time += timer.dt
print(total_time/N)

现在绘图本身仍然需要500毫秒，所以我认为imshow将是比scatter更好的选择。我怀疑matplotlib是否很乐意拥有25万个标记来渲染。