1

我正在开发一个程序，该程序模拟沿一维弦的波动以最终模拟不同的波包。我在“Python Scripting for Computational Science”一书中找到了一个程序，该程序声称可以描述波动，但我不确定如何实现它（这本书在 Google Books 上，不会显示之前/之后的文本代码）。

例如，我知道“f”是 x 和 t 的函数，“I”是 x 的函数，但实际需要哪些函数来产生波？

I= 
f= 
c= 
L= 
n=
dt=
tstop=


x = linespace(0,L,n+1) #grid points in x dir
dx = L/float(n)
if dt <= 0: dt = dx/float(c) #max step time
C2 = (c*dt/dx)**2  #help variable in the scheme
dt2 = dt*dt

up = zeros(n+1) #NumPy solution array
u = up.copy()   #solution at t-dt
um = up.copy()    #solution at t-2*dt

t = 0.0
for i in iseq(0,n):
u[i] +0.5*C2*(u[i-1] - 2*u[i] +u[i+1]) + \
    dt2*f(x[i], t)
um[0] = 0;   um[n] = 0

while t<= tstop:
t_old = t; t+=dt
#update all inner points:
for i in iseq(start=1, stop= n-1):
up[i] = -um[i] +2*u[i] + \
    C2*(u[i-1] - 2*u[i] + u[i+1]) + \
    dt2*f(x[i], t_old)

#insert boundary conditions
up[0] = 0; up[n] = 0
#updata data structures for next step
um = u.copy(); u = up.copy()
4

1 回答 1

2

下面的代码应该可以工作：

from math import sin, pi
from numpy import zeros, linspace
from scitools.numpyutils import iseq

def I(x):   
    return sin(2*x*pi/L)  

def f(x,t): 
    return 0

def solver0(I, f, c, L, n, dt, tstop):
    # f is a function of x and t, I is a function of x
    x = linspace(0, L, n+1) # grid points in x dir
    dx = L/float(n)
    if dt <= 0: dt = dx/float(c) # max time step
    C2 = (c*dt/dx)**2 # help variable in the scheme
    dt2 = dt*dt

    up = zeros(n+1) # NumPy solution array
    u = up.copy() # solution at t-dt
    um = up.copy() # solution at t-2*dt

    t = 0.0
    for i in iseq(0,n):
        u[i] = I(x[i])
    for i in iseq(1,n-1):
        um[i] = u[i] + 0.5*C2*(u[i-1] - 2*u[i] + u[i+1]) + \
                dt2*f(x[i], t)

    um[0] = 0; um[n] = 0

    while t <= tstop:
          t_old = t; t += dt
          # update all inner points:
          for i in iseq(start=1, stop=n-1):
              up[i] = - um[i] + 2*u[i] + \
                      C2*(u[i-1] - 2*u[i] + u[i+1]) + \
                      dt2*f(x[i], t_old)

          # insert boundary conditions:
          up[0] = 0; up[n] = 0
          # update data structures for next step
          um = u.copy(); u = up.copy()
    return u

if __name__ == '__main__':

# When choosing the parameters you should also check that the units are correct

   c = 5100
   L = 1
   n = 10
   dt = 0.1
   tstop = 1
   a = solver0(I, f, c, L, n, dt, tstop)

它返回一个数组，其中包含时间 tstop 以及我们解决方案网格中所有点的波浪值。

在将其应用于实际情况之前，您应该阅读波动方程和有限元方法以了解代码的作用。它可用于求波动方程的数值解：

Utt + beta*Ut = c^2*Uxx + f(x,t)

这是物理学中最重要的微分方程之一。此 PDE 或波的解由一个函数给出，该函数是空间和时间的函数u(x,t)。

为了形象化波浪的概念，考虑两个维度，空间和时间。如果你固定时间，例如 t1，你会得到一个 x 的函数：

U(x) = U(x,t=t1)

然而，在空间的特定点 x1 处，波是时间的函数：

U(t) = U(x=x1, t)

这应该可以帮助您了解波的传播方式。为了找到解决方案，您需要施加一些初始条件和边界条件，以将所有可能的波限制为您感兴趣的波。对于这种特殊情况：

  • I = I(xi)是我们将施加以使扰动/波前进的初始力。
  • 该术语f = f(x,t)解释了任何产生波浪的外力。
  • c是波速；它是一个常数（假设介质是均匀的）。
  • L是我们要求解 PDE 的域的长度；也是一个常数。
  • n是空间中的网格点数。
  • dt是一个时间步长。
  • tstop是停止时间。
于 2012-04-17T00:50:46.557 回答