1.1 波動方程式の導出
弦の振動から1次元波動方程式を導出します。
📐 理論
1次元波動方程式:
\[\frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2}\]ここで \(u(x,t)\) は変位、\(c\) は波の伝播速度
導出: 張力 \(T\) の弦、線密度 \(\rho\) のとき \(c = \sqrt{T/\rho}\)
💻 コード例 1: 波動方程式の数値シミュレーション（有限差分法）
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
# パラメータ
L = 10.0 # 弦の長さ
c = 1.0 # 波の伝播速度
T = 20.0 # 時間範囲
Nx = 200 # 空間格子点数
Nt = 1000 # 時間ステップ数
# 格子
x = np.linspace(0, L, Nx)
t = np.linspace(0, T, Nt)
dx = x[1] - x[0]
dt = t[1] - t[0]
# CFL条件の確認
r = c * dt / dx # クーラント数
print(f"クーラント数 r = {r:.3f} (安定条件: r ≤ 1)")
# 初期条件: 三角波
def initial_displacement(x):
u = np.zeros_like(x)
peak = L / 2
width = L / 5
mask = np.abs(x - peak) < width
u[mask] = 1 - np.abs(x[mask] - peak) / width
return u
# 初期化
u = np.zeros((Nt, Nx))
u[0] = initial_displacement(x)
u[1] = u[0].copy() # 初速度=0
# 有限差分法（陽的スキーム）
for n in range(1, Nt-1):
for i in range(1, Nx-1):
u[n+1, i] = (2*(1-r**2)*u[n, i] - u[n-1, i] +
r**2*(u[n, i+1] + u[n, i-1]))
# 境界条件: u(0,t) = u(L,t) = 0
u[n+1, 0] = 0
u[n+1, -1] = 0
# 可視化
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# スナップショット
times = [0, Nt//4, Nt//2, 3*Nt//4]
for idx, n in enumerate(times):
ax = axes[idx//2, idx%2]
ax.plot(x, u[n], 'b-', linewidth=2)
ax.axhline(0, color='gray', linewidth=0.5)
ax.grid(True, alpha=0.3)
ax.set_xlabel('位置 x', fontsize=12)
ax.set_ylabel('変位 u(x,t)', fontsize=12)
ax.set_title(f't = {t[n]:.2f}', fontsize=12)
ax.set_ylim(-1.2, 1.2)
plt.suptitle('波動方程式の数値解（有限差分法）', fontsize=14)
plt.tight_layout()
plt.show()
# 時空間プロット
plt.figure(figsize=(12, 6))
plt.contourf(x, t, u, levels=50, cmap='RdBu_r')
plt.colorbar(label='変位 u(x,t)')
plt.xlabel('位置 x', fontsize=12)
plt.ylabel('時間 t', fontsize=12)
plt.title('波動の伝播（時空間図）', fontsize=14)
plt.tight_layout()
plt.show()
print("\n波動方程式の性質:")
print("- 波は左右に伝播し、境界で反射")
print("- エネルギーは保存される")
print(f"- 伝播速度: c = {c} m/s")