Python: Написание программы для имитации движения 1D волны вдоль строки

Question

Я работаю над программой, которая имитирует движение волны вдоль одномерной строки, чтобы в конечном итоге имитировать различные волновые пакеты. Я нашел программу в книге «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() 

Answer 1

ниже код должен работать:

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, вы получите функцию х:

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 - время остановки.