1

私は、最終的にさまざまな波束をシミュレートするために、1 次元の文字列に沿って波動をシミュレートするプログラムに取り組んでいます。「計算科学のための Python スクリプティング」という本の中で、波動を記述すると主張するプログラムを見つけましたが、それを実装する方法はよくわかりません (この本は Google ブックスにあり、その前後のテキストは表示されません)。コード）。

たとえば、「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)

これは、物理学で最も重要な微分方程式の1つです。この偏微分方程式または波の解は、空間と時間の関数である関数によって与えられu(x,t)ます。

波の概念を視覚化するために、空間と時間の2つの次元を考えてみましょう。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停止時間です。
于 2012-04-17T00:50:46.557 に答える