#!/usr/bin/env python # coding: utf-8 # 波動方程式の解析解シミュレーション # ==== # [HTML5による物理シミュレーション 拡散・波動編](http://www.amazon.co.jp/HTML5%E3%81%AB%E3%82%88%E3%82%8B%E7%89%A9%E7%90%86%E3%82%B7%E3%83%9F%E3%83%A5%E3%83%AC%E3%83%BC%E3%82%B7%E3%83%A7%E3%83%B3-%E6%8B%A1%E6%95%A3%E3%83%BB%E6%B3%A2%E5%8B%95%E7%B7%A8%E2%80%95JavaScript%E3%83%A9%E3%82%A4%E3%83%96%E3%83%A9%E3%83%AA%E3%81%A8Canvas-Context-Web-Workers%E3%82%92%E4%BD%BF%E3%81%86/dp/4877833129/ref=sr_1_2?ie=UTF8&qid=1412059790&sr=8-2&keywords=html5%E3%81%AB%E3%82%88%E3%82%8B%E7%89%A9%E7%90%86%E3%82%B7%E3%83%9F%E3%83%A5%E3%83%AC%E3%83%BC%E3%82%B7%E3%83%A7%E3%83%B3)の内容をなぞった．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 = """""" 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)) # 定常波となることが分かる．次に，$A=2,B=1$で行う． # In[4]: # 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}$$ # # ### 自由端条件 # $$\begin{align} # A_n&=\frac{2}{L}\int^{L/2}_{-L/2}u(x,0)\cos\left(\frac{2\pi n}{L}x\right)dx\\ # B_n&=\frac{2}{L}\int^{L/2}_{-L/2}u(x,0)\sin\left(\frac{2\pi}{L}\left(n-\frac{1}{2}\right)x\right)dx\\ # 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 n}{L}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}{L}\left(n-\frac{1}{2}\right)x\right)dx # \end{align}$$ # # 固定端条件の場合と同様に，波動全体の周期は # $$T=\frac{2L}{v}$$ # # 波動方程式の解析解シミュレータの作成 # ---- # 以上をもとに，IPython上で一次元波動シミュレータを実装してみる． # In[5]: 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)) # 波動方程式の解析解シミュレーション # ---- # ### 初期空間分布のみ与えた場合 # シミュレータを用いて，それぞれの境界条件でシミュレーションを行う．初期空間分布は，ガウス分布 # $$ # u(x,0)=\frac{1}{\sqrt{2\pi \sigma^2}}\exp\left(-\frac{(x-\bar{x})^2}{2\sigma^2}\right) # \hspace{3em}(\bar{x}=0.1,\sigma^2=0.01) # $$ # を利用する．また，波動方程式の$v$(波の伝播速度)を1，波長$L=1$とした． # # #### 周期端条件 # In[6]: 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$で位相が逆になる． # # #### 固定端条件 # In[7]: 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． # # #### 自由端条件 # In[8]: plot_wave(u0, v0, 100, 0.02, 10, 1, 'n', (-0.5, 0.5), (0, 4)) # 自由端条件を満たしていることが分かる．周期は2． # # #### 階段関数 # 少し変わった波を初期空間分布として設定してみる． # In[9]: 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)) # ### 初期速度分布のみ与えた場合 # 初期条件として， # $$ # \left.\frac{\partial u(x,t)}{\partial t}\right|_{t=0}=-\frac{1}{\sqrt{2\pi \sigma^2}}\exp\left(-\frac{(x-\bar{x})^2}{2\sigma^2}\right) # \hspace{3em}(\bar{x}=0.1,\sigma^2=0.01) # $$ # を与える． # # #### 周期端条件 # In[10]: 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)) # #### 固定端条件 # In[11]: plot_wave(u0, v0, 101, 0.02, 10, 1, 'd', (-0.5, 0.5), (-0.6, 0.6)) # #### 自由端条件 # In[12]: 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