ディルクレの波動方程式

ディルクレの波動方程式

ディルクレの波動方程式は、数学や物理学で使用される偏微分方程式の一つで、波の伝播をモデル化します。

以下では、ディルクレの波動方程式の概要、物理的意味、数値解法について詳しく説明します。

ディルクレの波動方程式とは

波動方程式の一般形は次のように表されます：

$$
\frac{\partial^2 u}{\partial t^2} = c^2 \nabla^2 u
$$

  • $( u(x, y, t) ) $は位置 $ ( (x, y) ) $と時間 $ ( t ) $における波の振幅を示します。
  • $( c ) $は波の速度です。
  • $( \nabla^2 ) $はラプラス演算子であり、2次元の場合、次のように書かれます：

$$
\nabla^2 u = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}
$$

ディルクレ境界条件とは、領域の境界において関数$ ( u ) $がゼロになる条件を指します。
つまり、波動方程式の解$ ( u ) $が領域の端でゼロになるようにします。

物理的意味

波動方程式は、弦の振動や電磁波、音波など、さまざまな物理現象をモデル化します。
ディルクレ境界条件は、固定された端を持つ弦や膜の振動を表す場合によく使われます。
例えば、ギターの弦の両端が固定されている状態を考えると、弦の端での変位がゼロになるディルクレ条件を適用することができます。

数値解法の概要

波動方程式を数値的に解くためには、有限差分法がよく使用されます。
ここでは、空間と時間を離散化し、差分方程式を用いて波動方程式を解きます。
数値解法のステップは以下の通りです：

  1. 離散化:

    • 空間 $ ( x ) $と$ ( y ) $を格子状に分割し、時間 $ ( t ) $も小さなステップ$ ( dt ) $で分割します。

  2. 初期条件と境界条件の設定:

    • 初期条件として、時間 $( t=0 ) $における波の形状 $( u(x, y, 0) ) $とその時間微分 $( \frac{\partial u}{\partial t}(x, y, 0) ) $を設定します。
    • ディルクレ境界条件として、領域の境界上での波の振幅をゼロに設定します。

  3. 時間ステップごとの更新:

    • 中心差分法を用いて、各時間ステップにおける波の振幅を計算します。
      具体的には、次の差分方程式を使用します：

$$
u_{i,j}^{n+1} = 2u_{i,j}^n - u_{i,j}^{n-1} + \frac{c^2 dt^2}{dx^2} (u_{i+1,j}^n + u_{i-1,j}^n + u_{i,j+1}^n + u_{i,j-1}^n - 4u_{i,j}^n)
$$

ここで、$( u_{i,j}^n ) $は時間ステップ$ ( n ) $における格子点$ ( (i, j) ) $での波の振幅です。

プログラム例

ディルクレの波動方程式は、典型的には次のような形式をとる二次元の波動方程式です：

$$
\frac{\partial^2 u}{\partial t^2} = c^2 \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right)
$$

ここで、$( u(x, y, t) ) $は時間 $( t ) $における位置 $( (x, y) ) $での波の振幅、$ ( c )$ は波の速度です。

ディルクレ境界条件は、領域の境界上で$ ( u ) $がゼロになることを意味します。

この波動方程式を数値的に解いて、3次元グラフとして可視化します。

ここでは、単純な正方形領域内の初期条件と境界条件を使用します。

次のPythonコードを用いて、ディルクレの波動方程式を解き、3Dグラフをプロットします。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# パラメータの設定
L = 1.0   # 領域の長さ
T = 0.01  # 時間の長さ
c = 1.0   # 波の速度
dx = 0.01 # 空間ステップ
dt = 0.001 # 時間ステップ

# 空間格子と時間格子の生成
x = np.arange(0, L+dx, dx)
y = np.arange(0, L+dx, dx)
t = np.arange(0, T+dt, dt)
nx = len(x)
ny = len(y)
nt = len(t)

# 初期条件の設定
u = np.zeros((nx, ny, nt))
u[:, :, 0] = np.sin(np.pi * x[:, None]) * np.sin(np.pi * y[None, :])
u[:, :, 1] = u[:, :, 0]

# 波動方程式の数値解法（中心差分法）
for n in range(1, nt-1):
    for i in range(1, nx-1):
        for j in range(1, ny-1):
            u[i, j, n+1] = (2 * u[i, j, n] - u[i, j, n-1] +
                            c**2 * dt**2 / dx**2 *
                            (u[i+1, j, n] + u[i-1, j, n] + u[i, j+1, n] + u[i, j-1, n] - 4 * u[i, j, n]))

# 3Dプロットの作成
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

X, Y = np.meshgrid(x, y)
Z = u[:, :, nt//2]

ax.plot_surface(X, Y, Z, cmap='viridis')

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('U')
ax.set_title('Wave Equation at t = {:.3f}'.format(t[nt//2]))

plt.show()

このコードでは、波動方程式を数値的に解くために中心差分法を使用しています。

初期条件としては、正弦波の形状を使用しています。

また、ディルクレ境界条件を満たすために、境界上では常にゼロになるようにしています。

このスクリプトを実行すると、時間 $( t = T/2 ) $における波動の振幅 $( u(x, y, t) ) $を3Dプロットとして視覚化することができます。

[実行結果]

ソースコード解説

このPythonスクリプトは、ディルクレの波動方程式を2次元空間内で数値的に解き、時間のある時点での波動の振幅を3Dグラフとして可視化します。

以下では、スクリプトを章立てて詳しく説明します。

1. 必要なライブラリのインポート

1
2
3
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

この部分では、数値計算に必要なNumPyと、グラフ作成のためのMatplotlib、および3Dプロット用のツールキットをインポートしています。

2. パラメータの設定

1
2
3
4
5
6
# パラメータの設定
L = 1.0   # 領域の長さ
T = 0.01  # 時間の長さ
c = 1.0   # 波の速度
dx = 0.01 # 空間ステップ
dt = 0.001 # 時間ステップ

ここでは、計算領域のサイズ$ ( L )$、シミュレーションの時間長$ ( T )$、波の速度$ ( c )$、空間ステップ $( dx )$、および時間ステップ$ ( dt ) $を設定しています。

3. 空間格子と時間格子の生成

1
2
3
4
5
6
7
# 空間格子と時間格子の生成
x = np.arange(0, L+dx, dx)
y = np.arange(0, L+dx, dx)
t = np.arange(0, T+dt, dt)
nx = len(x)
ny = len(y)
nt = len(t)

ここでは、空間と時間の格子点を生成しています。

np.arangeを使って、$0$から$( L )$までの間をステップサイズ$( dx )$で分割した配列$ ( x ) $と$ ( y )$ を作成します。

時間についても同様に配列$ ( t ) $を作成します。

4. 初期条件の設定

1
2
3
4
# 初期条件の設定
u = np.zeros((nx, ny, nt))
u[:, :, 0] = np.sin(np.pi * x[:, None]) * np.sin(np.pi * y[None, :])
u[:, :, 1] = u[:, :, 0]

この部分では、波動の振幅を格納する3次元配列$ ( u ) $を初期化しています。

初期条件として、正弦波の形状を与えています。

x[:, None]とy[None, :]を用いて二次元配列を作成し、初期時刻$ ( t=0 ) $での振幅を設定します。

次の時間ステップ$ ( t=dt ) $も同じ初期条件に設定します。

5. 波動方程式の数値解法（中心差分法）

1
2
3
4
5
6
7
# 波動方程式の数値解法（中心差分法）
for n in range(1, nt-1):
    for i in range(1, nx-1):
        for j in range(1, ny-1):
            u[i, j, n+1] = (2 * u[i, j, n] - u[i, j, n-1] +
                            c**2 * dt**2 / dx**2 *
                            (u[i+1, j, n] + u[i-1, j, n] + u[i, j+1, n] + u[i, j-1, n] - 4 * u[i, j, n]))

この部分では、中心差分法を用いて波動方程式を数値的に解いています。

各時間ステップ$ ( n ) $に対して、空間格子点$ ( (i, j) ) $での波動の振幅を更新します。

中心差分法を用いることで、時間と空間の二階微分を近似的に計算しています。

6. 3Dプロットの作成

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
# 3Dプロットの作成
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

X, Y = np.meshgrid(x, y)
Z = u[:, :, nt//2]

ax.plot_surface(X, Y, Z, cmap='viridis')

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('U')
ax.set_title('Wave Equation at t = {:.3f}'.format(t[nt//2]))

plt.show()

最後に、時間$ ( t = T/2 ) $における波動の振幅を3Dプロットとして表示します。

np.meshgridを用いて2次元の格子点$ ( X ) $と$ ( Y ) $を作成し、それに対応する波動の振幅 $( Z )$ を取得します。

plot_surfaceを用いて3Dプロットを描画し、軸ラベルとタイトルを設定しています。

総括

このスクリプトは、ディルクレの波動方程式を2次元空間内で数値的に解き、特定の時点での波動の振幅を3Dグラフとして視覚化するものです。

波動方程式の数値解法として中心差分法を使用し、初期条件として正弦波を設定しています。

結果として得られた波動の振幅は、時間の中間点で3Dプロットとして表示されます。

グラフ解説

[実行結果]

以下に、グラフに表示される内容を詳細に説明します。

1. X軸とY軸:

  • グラフの$X$軸と$Y$軸は、空間内の位置を表しています。
    具体的には、正方形領域内の点$ ( (x, y) ) $を示しています。
    この領域は、$0$から$1$の範囲に設定されています。

2. Z軸:

  • $Z$軸は、波動の振幅 $( u(x, y, t) ) $を示しています。
    これは、時間のある特定の時点$ ( t ) $における空間の各点での波の高さを意味します。

3. 波動の振幅の変化:

  • グラフは、時間 $( t = T/2 ) $における波動の振幅を3次元的に表示しています。

波動の振幅は、初期条件として与えた正弦波形に基づいて、時間と共に変化します。
この特定の時間における波の形状が、3Dプロットとして描画されます。

4. 色の変化:

  • プロットの色は振幅の高さに対応しています。
    通常、色のスキーム（ここではviridisカラーマップ）が使用され、低い振幅は寒色（青や緑）で、高い振幅は暖色（黄やオレンジ）で示されます。