HTML5による物理シミュレーション 拡散・波動編の内容をなぞった.JavaScriptが書きたくなかったので,Pythonで書いた.波動方程式の導出や解析解の導出過程は理解したものを本を参考にしながら書いた.間違いもあるかもしれない.
$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$で後進しているといえることが分かる.以上から,波動方程式の表す波動は,前進する波と後進する波の合成によって構成されることが分かる.前者を前進波,後者を後進波という.
前進波と後進波によって波動が構成される様子を示す.まず描画の準備をする.
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
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$でシミュレーションする.緑線が前進波,青線が後進波,赤線が合成波を表す.
# 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))
定常波となることが分かる.次に,$A=2,B=1$で行う.
# 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=(-3, 3), 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 = 2 * 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))
前進波の影響を大きく受けて,合成波が前方に移動している.
波動方程式の基本解 $$u(x,t)=A(k)e^{ik(x-vt)}+B(k)e^{-ik(x+vt)}\hspace{3em}(A,B\in \mathbb{C})$$ の異なる$k$の重ね合わせを積分で表すと $$u(x,t)=\int^\infty_{-\infty}\left[A(k)e^{ik(x-vt)}+B(k)e^{-ik(x+vt)}\right]dk$$ $A(k),B(k)$を初期条件から決定することで,任意の時刻の波の状態を計算できる.
また,基本解を指数表記したものの実数部をとることにより, $$ \begin{align} u(x,t)&=A\cos(kx-\omega t)+iA\sin(kx-\omega t)+B\cos(kx-\omega t)+iB\sin(kx-\omega t)\\ &=(C_1\cos(kx)+C_2\sin(kx))\cos(\omega t)+(C_3\cos(kx)+C_4\sin(kx))\sin(\omega t) \end{align} $$ とできる.$A,B\in\mathbb{C}$に注意.また,$C_1,C_2,C_3,C_4\in\mathbb{R}$.このとき,一般解は $$u(x,t)=\int^\infty_{-\infty}\left[(C_1(k)\cos(kx)+C_2(k)\sin(kx))\cos(\omega t)+(C_3(k)\cos(kx)+C_4(k)\sin(kx))\sin(\omega t)\right]dk$$
代表的な境界条件には以下のものがある.
一次元弦の変位を$u(x,t)$とすると,各境界条件は
周期端条件: $u(x+L,t)=u(x,t)$
固定端条件: $u(-L/2,t)=u(L/2,t)=0$
自由端条件: $\left.\frac{\partial u(x,t)}{\partial x}\right|_{x=-L/2}=0,\left.\frac{\partial u(x,t)}{\partial x}\right|_{x=L/2}=0$
で与えられる.
周期的境界条件を課す場合,$k$は離散的な値 $$k_n=\frac{2\pi n}{L}\hspace{3em}(n\in \mathbb{N})$$ しかとれない.このことから,一般解は $$ \begin{align} u(x,t)&=\sum^\infty_{n=-\infty}[(A_n\cos(k_nx)+B_n\sin(k_nx))\cos(\omega_nt)+(C_n\cos(k_nx)+D_n\sin(k_nx))\sin(\omega_nt)]\\ &=\frac{A_0}{2}+\sum^\infty_{n=1}[(A_n\cos(k_nx)+B_n\sin(k_nx))\cos(\omega_nt)+(C_n\cos(k_nx)+D_n\sin(k_nx))\sin(\omega_nt)] \end{align} $$ 三角関数の対称性を利用した.定数は適宜置き換えているので,注意して欲しい.上の$A_n$と下の$A_n$は異なる.一般解を構成する各項は基準振動と呼ばれ,波動現象はこれらの重ねあわせにより表現できる.
${}^\forall x,{}^\forall t$に対して境界条件が成り立つようにするには $$ \begin{align} \cos: \cos\left(-k\frac{L}{2}\right)&=\cos\left(k\frac{L}{2}\right)&=0\\ \sin: \sin\left(-k\frac{L}{2}\right)&=\sin\left(k\frac{L}{2}\right)&=0 \end{align} $$ でなければならないので, $$ \begin{align} \cos: k^{(c)}_n&=\frac{2\pi}{L}\left(n-\frac{1}{2}\right)\\ \sin: k^{(s)}_n&=\frac{2\pi n}{L} \end{align} $$ つまり,$\cos(kx),\sin(kx)$を別々に考えて境界条件を課さなければならない.よってそれぞれは異なる角振動数 $$\omega^{(c)}_n=vk^{(c)}_n,\omega^{(s)}_n=vk^{(s)}_n$$ で振動する事になる.以上より一般解は $$ \begin{align} u(x,t)&=\sum^\infty_{n=-\infty}[A_n\cos(k^{(c)}_nx)\cos(\omega^{(c)}_nt)+B_n\sin(k^{(s)}_nx)\cos(\omega^{(s)}_nt)+C_n\cos(k^{(c)}_nx)\sin(\omega^{(c)}_nt)+D_n\sin(k^{(s)}_nx)\sin(\omega^{(s)}_nt)]\\ &=\sum^\infty_{n=1}[A_n\cos(k^{(c)}_nx)\cos(\omega^{(c)}_nt)+B_n\sin(k^{(s)}_nx)\cos(\omega^{(s)}_nt)+C_n\cos(k^{(c)}_nx)\sin(\omega^{(c)}_nt)+D_n\sin(k^{(s)}_nx)\sin(\omega^{(s)}_nt)] \end{align} $$
自由端条件の場合も$\cos(kx),\sin(kx)$に対する$k$の条件を別々に考えることになる. $$ \begin{align} \cos: \sin\left(-k\frac{L}{2}\right)&=\sin\left(k\frac{L}{2}\right)&=0\\ \sin: \cos\left(-k\frac{L}{2}\right)&=\cos\left(k\frac{L}{2}\right)&=0 \end{align} $$ でなければならないので,$n\in\mathbb{Z}$を用いて $$ \begin{align} \cos: k^{(c)}_n&=\frac{2\pi n}{L}\\ \sin: k^{(s)}_n&=\frac{2\pi}{L}\left(n-\frac{1}{2}\right) \end{align} $$ 固定端条件のときと条件が$\cos$と$\sin$で逆になっている.以上から,一般解は $$ \begin{align} u(x,t)&=\sum^\infty_{n=-\infty}[A_n\cos(k^{(c)}_nx)\cos(\omega^{(c)}_nt)+B_n\sin(k^{(s)}_nx)\cos(\omega^{(s)}_nt)+C_n\cos(k^{(c)}_nx)\sin(\omega^{(c)}_nt)+D_n\sin(k^{(s)}_nx)\sin(\omega^{(s)}_nt)]\\ &=\frac{A_0}{2}+\sum^\infty_{n=1}[A_n\cos(k^{(c)}_nx)\cos(\omega^{(c)}_nt)+B_n\sin(k^{(s)}_nx)\cos(\omega^{(s)}_nt)+C_n\cos(k^{(c)}_nx)\sin(\omega^{(c)}_nt)+D_n\sin(k^{(s)}_nx)\sin(\omega^{(s)}_nt)] \end{align} $$
周期端条件の一般解に$t=0$を代入し, $$u(x,0)=\frac{A_n}{2}+\sum^\infty_{n=1}[A_n\cos(k_nx)+B_n\sin(k_nx)]\hspace{2em}(k_n=\frac{2\pi n}{L})$$ ここで,フーリエ級数展開を思い出してみる.関数$f(x)$が $$f(x)=\frac{a_0}{2}+\sum^\infty_{n=1}[a_n\cos\left(\frac{n\pi}{L}x\right)+b_n\sin\left(\frac{n\pi}{L}x\right)]$$ と展開されるとき,$a_n,b_n$は三角関数の直交性から $$\begin{align} a_n&=\frac{1}{L}\int^L_{-L}f(x)\cos\left(\frac{n\pi}{L}x\right)dx\\ b_n&=\frac{1}{L}\int^L_{-L}f(x)\sin\left(\frac{n\pi}{L}x\right)dx \end{align}$$ と導けた.$f(x)$と$u(x,0)$の表式を見比べると,$x$が$2x$になっているだけなので,$a_n,b_n$の$x$を$2x$に置き換えて $$\begin{align} &A_n=\frac{1}{L}\int^{L/2}_{-L/2}f(2x)\cos\left(\frac{2\pi n}{L}x\right)\cdot2dx&=\frac{2}{L}\int^{L/2}_{-L/2}u(x,0)\cos(k_n x)dx\\ &B_n&=\frac{2}{L}\int^{L/2}_{-L/2}u(x,0)\sin(k_n x)dx \end{align}$$ これで,空間分布の初期条件を利用することで,$A_n,B_n$が求まった.しかし,$C_n,D_n$がまだ求まっていない.そこで,$u(x,t)$の時間微分を求め,そこから$C_n,D_n$を求める.つまり,速度分布の初期条件を利用して$C_n,D_n$を求めてみる. $$\begin{align} \left.\frac{\partial u(x,t)}{\partial t}\right|_{t=0}&=\sum^\infty_{n=1}\left[C_n\cos(k_nx)+D_n\sin(k_nx)\right]\omega_n\\ C_n&=\frac{2}{\omega_nL}\int^{L/2}_{-L/2}\left.\frac{\partial u(x,t)}{\partial t}\right|_{t=0}\cos(k_nx)dx\\ D_n&=\frac{2}{\omega_nL}\int^{L/2}_{-L/2}\left.\frac{\partial u(x,t)}{\partial t}\right|_{t=0}\sin(k_nx)dx \end{align}$$ つまり,波動方程式の一般解の展開系数は,初期空間分布と初期速度分布より求めることができる.
またこのとき,各基準振動に対する周期は,一般解式より $$T_n=\frac{2\pi}{\omega_n}=\frac{L}{nv}$$ この式から $$T_n=\frac{T_1}{n}$$ であるので,$n\geq1$に対して最大の周期は$T_1$であり,全ての周期は最大の周期の整数分の1.よって波動全体の周期は $$T=\frac{L}{v}$$
周期端条件の場合とほぼ同様. $$\begin{align} u(x,0)&=\sum^\infty_{n=1}\left[A_n\cos\left(\frac{2\pi}{L}x\left(n-\frac{1}{2}\right)\right)+B_n\sin\left(\frac{2\pi n}{L}x\right)\right]\\ A_n&=\frac{2}{L}\int^{L/2}_{-L/2}u(x,0)\cos\left(\frac{2\pi}{L}\left(n-\frac{1}{2}\right)x\right)dx\\ B_n&=\frac{2}{L}\int^{L/2}_{-L/2}u(x,0)\sin\left(\frac{2\pi n}{L}x\right)dx\\ \left.\frac{\partial u(x,t)}{\partial t}\right|_{t=0}&=\sum^\infty_{n=1}\left[\omega^{(c)}_nC_n\cos\left(\frac{2\pi}{L}x\left(n-\frac{1}{2}\right)\right)+\omega^{(s)}_nD_n\sin\left(\frac{2\pi n}{L}x\right)\right]\\ C_n&=\frac{2}{\omega^{(c)}_nL}\int^{L/2}_{-L/2}\left.\frac{\partial u(x,t)}{\partial t}\right|_{t=0}\cos\left(\frac{2\pi}{L}\left(n-\frac{1}{2}\right)x\right)dx\\ D_n&=\frac{2}{\omega^{(s)}_nL}\int^{L/2}_{-L/2}\left.\frac{\partial u(x,t)}{\partial t}\right|_{t=0}\sin\left(\frac{2\pi n}{L}x\right)dx \end{align}$$
またこのとき,各基準振動に対する周期は,一般解式より $$\begin{align} T^{(c)}_n&=\frac{2\pi}{\omega^{(c)}_n}=\frac{L}{(n-1/2)v}\\ T^{(s)}_n&=\frac{2\pi}{\omega^{(s)}_n}=\frac{L}{nv} \end{align}$$ 最長の周期は$T=T^{(c)}_1$であり,全ての周期は最大の周期の整数分の1なので,波動全体の周期は $$T=\frac{2L}{v}$$
固定端条件の場合と同様に,波動全体の周期は $$T=\frac{2L}{v}$$
以上をもとに,IPython上で一次元波動シミュレータを実装してみる.
from scipy.integrate import quad
def plot_wave(u0, v0, step, delta_t, max_n, wave_v, boundary, xlim, ylim):
L = xlim[1] - xlim[0]
# calculate coefficient of solution of wave equation
if boundary in ('periodic', 'p'): # 周期端条件
kc = ks = lambda n: 2 * np.pi * n / L
elif boundary in ('dirichlet', 'd'): # 固定端条件
kc = lambda n: 2 * np.pi * (n - 0.5) / L
ks = lambda n: 2 * np.pi * n / L
elif boundary in ('neumann', 'n'): # 自由端条件
kc = lambda n: 2 * np.pi * n / L
ks = lambda n: 2 * np.pi * (n - 0.5) / L
A = lambda n: 2 / L * quad(lambda x: u0(x) * np.cos(kc(n) * x), -L/2, L/2)[0]
B = lambda n: 2 / L * quad(lambda x: u0(x) * np.sin(ks(n) * x), -L/2, L/2)[0]
C = lambda n: 2 / (kc(n) * wave_v * L) * quad(lambda x: v0(x) * np.cos(kc(n) * x), -L/2, L/2)[0]
D = lambda n: 2 / (ks(n) * wave_v * L) * quad(lambda x: v0(x) * np.sin(ks(n) * x), -L/2, L/2)[0]
# calculate time evolution
range_n = np.linspace(1, max_n, max_n)
u_sum = lambda x, t: sum([A(n) * np.cos(kc(n) * x) * np.cos(kc(n) * wave_v * t) + B(n) * np.sin(ks(n) * x) * np.cos(ks(n) * wave_v * t)
+ C(n) * np.cos(kc(n) * x) * np.sin(kc(n) * wave_v * t) + D(n) * np.sin(ks(n) * x) * np.sin(ks(n) * wave_v * t)
for n in range_n])
if boundary in ('periodic', 'p', 'neumann', 'n'):
u = lambda x, t: A(0) / 2 + u_sum(x, t)
elif boundary in ('dirichlet', 'd'):
u = u_sum
# set up the figure, the axis, and the plot element we want to animate
fig = plt.figure()
fig.set_size_inches(12, 4)
ax = plt.subplot(121, xlim=xlim, ylim=ylim, xlabel=r"$x$", ylabel=r"$u(x,t)$")
line, = ax.plot([], [], '-', lw=2)
plt.subplot(122)
plt.title("Coefficients")
plt.xlabel("$n$")
plt.ylabel("Coefficient")
plt.grid()
plt.plot(range_n, [A(n) for n in range_n], label=r"$A_n$")
plt.plot(range_n, [B(n) for n in range_n], label=r"$B_n$")
plt.plot(range_n, [C(n) for n in range_n], label=r"$C_n$")
plt.plot(range_n, [D(n) for n in range_n], label=r"$D_n$")
plt.legend(loc='best').get_frame().set_alpha(0.8)
# initialization function: plot the background of each frame
def init():
line.set_data([], [])
return line,
range_x = np.linspace(xlim[0], xlim[1], 100)
# animation function. This is called sequentially
def animfunc(t):
t *= delta_t
ax.set_title("t={0:.2f}".format(t))
line.set_data(range_x, u(range_x, t))
return line,
# call the animator. blit=True means only re-draw the parts that have changed.
return display_anim(FuncAnimation(fig, animfunc, init_func=init, frames=step, interval=100, blit=True, repeat=False))
def u0(x):
mean_x, sigma2 = 0.1, 0.01
return np.exp(-(x - mean_x)**2 / (2 * sigma2)) / np.sqrt(2 * np.pi * sigma2)
v0 = lambda x:0
plot_wave(u0, v0, 100, 0.02, 10, 1, 'p', (-0.5, 0.5), (0, 4))
両端における傾きと値が一致していることが確認できる.周期は $$T=\frac{L}{v}=1$$ なので,$t=0.5$で位相が逆になる.
plot_wave(u0, v0, 100, 0.02, 10, 1, 'd', (-0.5, 0.5), (-4, 4))
任意の時刻で$u\left(-\frac{L}{2},t\right)=u\left(\frac{L}{2},t\right)=0$であることが確認できる.周期は2.
plot_wave(u0, v0, 100, 0.02, 10, 1, 'n', (-0.5, 0.5), (0, 4))
u0 = lambda x: 1 if x > 0 else 0
plot_wave(u0, v0, 100, 0.02, 20, 1, 'd', (-0.5, 0.5), (-1, 1))
u0 = lambda x:0
def v0(x):
mean_x, sigma2 = 0.1, 0.01
return -np.exp(-(x - mean_x)**2 / (2 * sigma2)) / np.sqrt(2 * np.pi * sigma2)
plot_wave(u0, v0, 101, 0.02, 10, 1, 'p', (-0.5, 0.5), (-0.6, 0.6))
plot_wave(u0, v0, 101, 0.02, 10, 1, 'd', (-0.5, 0.5), (-0.6, 0.6))
plot_wave(u0, v0, 101, 0.02, 10, 1, 'n', (-0.5, 0.5), (-0.6, 0.6))
ちなみに,初期速度分布を $$ \left.\frac{\partial u(x,t)}{\partial t}\right|_{t=0} =\begin{cases} 1&(0<x<0.2)\\ 0&\mathrm{otherwise} \end{cases} $$ とすると,
u0 = lambda x: 0
v0 = lambda x: -1 if np.abs(x - 0.1) < 0.1 else 0
plot_wave(u0, v0, 101, 0.02, 30, 1, 'n', (-0.5, 0.5), (-0.3, 0.3))
/home/yuntan/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:16: IntegrationWarning: The integral is probably divergent, or slowly convergent. app.launch_new_instance() /home/yuntan/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:16: IntegrationWarning: The integral is probably divergent, or slowly convergent. app.launch_new_instance()