波動方程式の解析解シミュレーション

HTML5による物理シミュレーション 拡散・波動編の内容をなぞった.JavaScriptが書きたくなかったので,Pythonで書いた.波動方程式の導出や解析解の導出過程は理解したものを本を参考にしながら書いた.間違いもあるかもしれない.

Lagrangianを用いて波動方程式を導出する

$f$個のおもり(質量$m$)を自然長$a$の$f+1$本のばねでつないだ系を考える.それぞれのおもりの$x_j=ja$からの変位を$\eta_j\hspace{1em}(1\leq j\leq f)$とおく.系の全長を$L=a(f+1)$とおく.Lagrangianは $$\begin{align} \mathcal{L} &=K-U\\ &=\sum^f_{j=1}\frac{m}{2}\dot{\eta}^2_j-\left(\sum^{f-1}_{j=1}\frac{k}{2}(\eta_{j-1}-\eta_j)^2+\frac{k}{2}\eta^2_1+\frac{k}{2}\eta^2_f\right)\\ &=\frac{1}{2}\boldsymbol{\dot{\eta}}^T\hat{M}\boldsymbol{\dot{\eta}}-\frac{1}{2}\boldsymbol{\dot{\eta}}^T\hat{k}\boldsymbol{\dot{\eta}}\\ \end{align}$$ ただし, $$\boldsymbol{\dot{\eta}}=\begin{pmatrix}\eta_1\\\eta_2\\\vdots\\\eta_f\end{pmatrix}, \hat{M}=m\hat{1}=\begin{pmatrix} m&& && &&\boldsymbol{o}\\ &&m&& && \\ && &&\ddots&& \\ \boldsymbol{o}&& && &&m \end{pmatrix}, \hat{k}=k\begin{pmatrix} 2&&-1&& 0&&\cdots&&\cdots&&\cdots\\ -1&& 2&&-1&& 0&&\cdots&&\cdots\\ 0&&-1&& 2&&-1&& && \\ \vdots&& 0&&-1&&\ddots&& &&\\ \vdots&&\vdots&& && &&\ddots&&-1\\ \vdots&&\vdots&& && &&-1&&2 \end{pmatrix}$$ Lagrange方程式は $$\hat{M}\boldsymbol{\ddot{\eta}}+\hat{k}\boldsymbol{\eta}=0$$

ここで,全系の長さ$L$を保ちつつ,$f\rightarrow+\infty (a\rightarrow0)$の極限をとることを考える.このとき,$\eta_j(t)$を$u(x_j,t)$と置き換える.また,単位長さあたりの質量(質量密度)は保存するとしたいので,密度を$\rho$とおき, $$m=\rho a$$ とする.さらに,一様な材質のばねでは,バネ定数はばねの長さに反比例するため,系全体でのバネ定数$\kappa$というものを考え, $$\kappa L=ka$$ 以上を用いてLagrange方程式を書き直す.まず,Lagrange方程式を成分表示に直して, $$\begin{align} m&\frac{\partial^2}{\partial t^2}u(x_j,t)+k(2u(x_j,t)-u(x_{j+1},t)-u(x_{j-1},t))\hspace{2em}(j=1,2,\cdots,f)\\ \rho a&\frac{\partial^2}{\partial t^2}u(x_j,t)=-\kappa\frac{L}{a}(2u(x_j,t)-u(x_{j+1},t)-u(x_{j-1},t)) \end{align}$$ となる.次に,$u(x,t)$がなめらかであると考え, $$\begin{align} u(x_{j\pm1},t) &=u(x_j\pm a,t)\\ &=u(x_j,t)\pm\left.\frac{\partial u}{\partial x}\right|_{x=x_j}a+\frac{1}{2}\left.\frac{\partial^2 u}{\partial x^2}\right|_{x=x_j}a^2+O(a^3) \end{align}$$ と評価し,代入すると $$\rho a\frac{\partial^2}{\partial t^2}u(x_j,t)=\kappa\frac{L}{a}\left(\left.\frac{\partial^2 u}{\partial x^2}\right|_{x=x_j}a^2\right)$$ $a$を消去し,$x_j=x$と置き換えて $$\rho\frac{\partial^2}{\partial t^2}u(x,t)=\kappa L\frac{\partial^2}{\partial x^2}u(x,t)$$ $v=\sqrt{\frac{\kappa L}{\rho}}$とおくと, $$\frac{1}{v^2}\frac{\partial^2}{\partial t^2}u(x,t)=\frac{\partial^2}{\partial x^2}u(x,t)$$ この偏微分方程式を波動方程式とよぶ.

波動方程式の基本解

変数分離法を用いる.空間依存部分と時間依存部分を分離し, $$u(x,t)=X(x)T(t)$$ とする.これを波動方程式に代入すると, $$\begin{align} \frac{1}{v^2}X(x)\frac{d^2T(t)}{dt^2}&=T(t)\frac{d^2X(x)}{dx^2}\\ \frac{1}{v^2T(t)}\frac{d^2T(t)}{dt^2}=\frac{1}{X(x)}\frac{d^2X(x)}{dx^2} \end{align}$$ と分離できる.任意の$x$と$t$に対し上式が成立するためには,両辺の値が必ず定数となる必要がある.なので,分離定数を導入し, $$\frac{1}{v^2T(t)}\frac{d^2T(t)}{dt^2}=\frac{1}{X(x)}\frac{d^2X(x)}{dx^2}=(\mathrm{分離定数})$$ とできる.両辺は簡単に解け,解はそれぞれ $$\begin{align} X(x)&=A\mathrm{e}^{ikx}+B\mathrm{e}^{-ikx}\\ T(t)&=C\mathrm{e}^{i\omega t}+D\mathrm{e}^{-i\omega t} \end{align}$$ ただし,$A,B,C,D$は積分定数で,複素数.$k=\pm\sqrt{(\mathrm{分離定数})},\omega=vk$.$k$は先のバネ定数とは異なる.これらを$u(x,t)=X(x)T(t)$に代入し, $$u(x,t)=A\mathrm{e}^{ikx-i\omega t}+B\mathrm{e}^{-ikx-i\omega t}+C\mathrm{e}^{ikx+i\omega t}+D\mathrm{e}^{-ikx+i\omega t}$$ ただし,$A,B,C,D$は適宜置き換えている.ここで, $$\begin{align} A\mathrm{e}^{ikx-i\omega t}&=A(\cos(kx-\omega t)+i\sin(kx-\omega t))\\ D\mathrm{e}^{-ikx+i\omega t}&=D(\cos(kx-\omega t)-i\sin(kx-\omega t)) \end{align}$$ 定数$A,D$をうまく置き換えると二項を合わせることができる.B,Cについても同様.つまり $$u(x,t)=A\mathrm{e}^{ikx-i\omega t}+B\mathrm{e}^{-ikx-i\omega t}$$ とできる.

基本解から分かること

$\omega=vk$から $$u(x,t)=A\mathrm{e}^{ik(x-vt)}+B\mathrm{e}^{-ik(x+vt)}$$ 第一項,第二項をそれぞれ $$\begin{align} u_+(x-vt)&\equiv A\mathrm{e}^{ik(x-vt)}=A(\cos(k(x-vt))+i\sin(k(x-vt)))\\ u_-(x+vt)&\equiv B\mathrm{e}^{ik(x+vt)}=A(\cos(k(x+vt))-i\sin(k(x+vt))) \end{align}$$ と定義する.まず,$\cos(k(x-vt))$に着目する.この関数は時間発展しても,それに応じて位置が前進すると式の値が変わらない.つまり,この関数で表される波は,速度$v$で進行しているといえる.また,今度は$u_-$の$\cos(k(x+vt))$に着目すると,時間発展してもそれに応じて位置が後進すると式の値が変わらないことから,この関数で表される波が速度$v$で後進しているといえることが分かる.以上から,波動方程式の表す波動は,前進する波と後進する波の合成によって構成されることが分かる.前者を前進波,後者を後進波という.

前進波と後進波のシミュレーション

前進波と後進波によって波動が構成される様子を示す.まず描画の準備をする.

In [1]:
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
In [2]:
from tempfile import NamedTemporaryFile
import base64
from IPython.display import HTML

GIF_TAG = """<img src="data:image/gif;base64,{0}" type="image/gif">"""

def anim_to_html(anim):
    with NamedTemporaryFile(suffix='.gif') as f:
        anim.save(f.name, writer='imagemagick', fps=20)
        gif = open(f.name, "rb").read()
    return GIF_TAG.format(base64.b64encode(gif).decode())

def display_anim(anim):
    plt.close(anim._fig)
    return HTML(anim_to_html(anim))

前進波,後進波,合成波を $$\begin{align} u_+&=A\cos(kx-\omega t)\\ u_-&=B\cos(kx+\omega t)\\ u&=u_++u_- \end{align}$$ とする.また,$k=2\pi/10,\omega=2\pi/100$とした.

まず,$A=1.5,B=1.5$でシミュレーションする.緑線が前進波,青線が後進波,赤線が合成波を表す.

In [3]:
# First set up the figure, the axis, and the plot element we want to animate
fig = plt.figure()
ax = plt.axes(xlim=(0, 20), ylim=(-2, 2), xlabel=r"$x$", ylabel=r"$y$")
rline, = ax.plot([], [], '-', lw=2)  # 合成波
gline, = ax.plot([], [], '--', lw=2)  # 前進波
bline, = ax.plot([], [], '--', lw=2)  # 後進波

# initialization function: plot the background of each frame
def init():
    rline.set_data([], [])
    gline.set_data([], [])
    bline.set_data([], [])
    return rline, gline, bline

# animation function.  This is called sequentially
def animfunc(t):
    ax.set_title(r"$k=2\pi/10,\omega=2\pi/100,t={0}$".format(t))
    x = np.linspace(0, 20, 1000)
    uplus  = np.cos(2 * np.pi / 10 * x - 2 * np.pi / 100 * t)
    uminus = np.cos(2 * np.pi / 10 * x + 2 * np.pi / 100 * t)
    gline.set_data(x, uplus)
    bline.set_data(x, uminus)
    rline.set_data(x, uplus + uminus)
    return rline, gline, bline

# call the animator.  blit=True means only re-draw the parts that have changed.
display_anim(FuncAnimation(fig, animfunc, init_func=init, frames=101, interval=100, blit=True, repeat=False))
Out[3]: