Код ниже должен работать:
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)
которое является одним из важнейших дифференциальных уравнений в физике. Решение этого УЧП или волны задается функцией, которая является функцией пространства и времени 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 — длина домена, в котором мы хотим решить УЧП; также константа.n — количество точек сетки в пространстве.dt - шаг по времени.tstop - время остановки.
person
anzupe
schedule
17.04.2012
