SageMathで波動方程式の一般解を求める

# はじめに
　波の性質を表す波動方程式は偏微分方程式なので、変数分離を行って、常微分方程式として一般解を求めます。一般解を求めるためには、微分積分の公式集を見る必要があります。しかし、pythonの文法で記述可能なSegaMathとよばれるプログラムを使うことで、微分方程式の一般解をコンピュータに求めさせることができます。
　この記事では、最初に、高校でなじみのある正弦波の式を用いた古典的な波動方程式の立式を示し、その後、波動方程式の一般解をSegaMathで求める例を示します。

# 正弦波の式で波動方程式をつくる
　波動方程式の説明では、一般的な関数の
``` math
\begin{align}
y=f(x-ct)	
\end{align}
```
を用いた導出が知られていますが、具体的な波の式を使って導出するほうがピンとくると思うので、ここでは、弦の振動等を表す正弦波を例に波動方程式を導出します。
　高校の物理で取り上げられる正弦波の波の要素には、以下のものがあります。
``` math
\begin{align}
v &= \frac{\lambda}{T} \qquad    (1) \\
f &=\frac{1}{T}	\qquad    (2)\\
v &=f{\lambda} \qquad    (3)\\
\end{align}
```
ここで、$v$ [m/s]は波の速さ、$f$[Hz]は振動数、$\lambda$[m]は波長、$T$ [s]は周期を表しています。
　単振動する波源から、$x$軸の正の向きに速度$v$[m/s]で伝わる正弦波の原点での変位は、以下の式(4)になります。
``` math
\begin{align}
y=A\sin\frac{2\pi}{T}t \qquad    (4)
\end{align}
```
ここで、$y$は原点（$x=0$）にある波源の時刻$t$ [s]での変位、$A$[m]は振幅 、$T$[s]は周期を表します。
　さらに、時刻$t$ [s]における、位置$x$ [m]の媒質Pの変位$y$ [m]は、以下の式(5)のようになります。
``` math
\begin{align}
y=A\sin\frac{2\pi}{T}  (t-\frac{x}{v}) \qquad    (5)
\end{align}
```

ここで、$v$ [m/s]は波の進む速度、$\frac{x}{v}$は原点から位置$x$に振動が伝わるまでの時間の遅れを表します。
　式(5)を変形した式(6)を使って波動方程式を立式します。
　まず、式(6)を$x$で2回偏微分します。
``` math
\begin{align}
y &=A\sin(\frac{2\pi}{T}t-\frac{2\pi}{Tv}x) \qquad    (6)\\
\frac{\partial y}{\partial x} &=A(-\frac{2\pi}{Tv})\cos(\frac{2\pi}{T}t-\frac{2\pi}{Tv}x)\\
\frac{\partial^2 y}{\partial x^2}&=-A(-\frac{2\pi}{Tv})^2\sin(\frac{2\pi}{T}t-\frac{2\pi}{Tv}x)\\
\frac{\partial^2 y}{\partial x^2}&=-(-\frac{2\pi}{Tv})^2y  \qquad    (7) \\
\end{align}
```
　次に、式(6)を$t$で2回偏微分します。
``` math
\begin{align}
\frac{\partial y}{\partial t} &=A(-\frac{2\pi}{T})\cos(\frac{2\pi}{T}t-\frac{2\pi}{Tv}x)\\
\frac{\partial^2 y}{\partial t^2} &=-A(-\frac{2\pi}{T})^2\sin(\frac{2\pi}{T}t-\frac{2\pi}{Tv}x)\\
\frac{\partial^2 y}{\partial t^2} &=-(-\frac{2\pi}{T})^2y \qquad    (8)\\
\end{align}
```
式(7)を以下のように書き直します。
``` math
\begin{align}
\frac{\partial^2 y}{\partial x^2} &=-\frac{1}{v^2}(-\frac{2\pi}{T})^2y \qquad    (9)\\
\end{align}
```
式(8)と式(9)を比較すると、

$\frac{\partial^2 y}{\partial t^2}
$に$\frac{1}{v^2}$をかけると、$\frac{\partial^2 y}{\partial x^2}$に一致することがわかり、以下のようになります。

``` math
\begin{align}
\frac{\partial^2 y}{\partial x^2}=\frac{1}{v^2}\frac{\partial^2 y}{\partial t^2}
\end{align}
```
同様の導出を
``` math
\begin{align}
y=f(x-ct)	
\end{align}
```
の関数で行います。この式と正弦波の式との違いは、正弦波の式の変数の単位が時間の単位であったのに対し、変数の単位が距離（位置）の単位であることです。
``` math
\begin{align}
\frac{\partial y}{\partial x} &=f'(x-ct)\\
\frac{\partial^2 y}{\partial x^2}&=f''(x-ct)\\
\frac{\partial y}{\partial t}&=-cf'(x-ct)\\
\frac{\partial^2 y}{\partial t^2}&=c^2f''(x-ct)\\
\frac{\partial^2 y}{\partial t^2}&=c^2\frac{\partial^2 y}{\partial x^2}\\
\frac{\partial^2 y}{\partial x^2}&=\frac{1}{c^2}\frac{\partial^2 y}{\partial t^2}
\end{align}
```

　$y=f(x-ct)$は波動方程式の解となります。
　次に、$y=f(x-ct)$は$x$軸の正の方向に進む波を表しているので、$x$軸の負の方向に進む波$y=g(x+ct)$が波動方程式を満たすかを確認します。
``` math
\begin{align}
\frac{\partial y}{\partial x}&=f'(x+ct)\\
\frac{\partial^2 y}{\partial x^2}&=f''(x+ct)\\
\frac{\partial y}{\partial t}&=cf'(x+ct)\\
\frac{\partial^2 y}{\partial t^2}&=c^2f''(x+ct)\\
\frac{\partial^2 y}{\partial t^2}&=c^2\frac{\partial^2 y}{\partial x^2}\\
\frac{\partial^2 y}{\partial x^2}&=\frac{1}{c^2}\frac{\partial^2 y}{\partial t^2}
\end{align}
```
 　最後に、2つの波が線形結合された$y=f(x-ct)+g(x+ct)$が波動方程式を満たすかを確認します。
``` math
\begin{align}
\frac{\partial y}{\partial x}&=f'(x-ct)+g'(x+ct)\\
\frac{\partial^2 y}{\partial x^2}&=f''(x-ct)+g''(x+ct)\\
\frac{\partial y}{\partial t}&=-cf'(x-ct)+cg''(x+ct)\\
\frac{\partial^2 y}{\partial t^2}&=c^2f''(x-ct)+c^2g''(x+ct)\\
\frac{\partial^2 y}{\partial t^2}&=c^2\frac{\partial^2 y}{\partial x^2}\\
\frac{\partial^2 y}{\partial x^2}&=\frac{1}{c^2}\frac{\partial^2 y}{\partial t^2}
\end{align}
```
$y=f(x-ct)+g(x+ct)$も波動方程式の解であることがわかりました。この解のことを「ダランベールの解」とよびます。

# 1次元の波動方程式の解
　解を具体的に与えてみます。
　まず、2つの変数のある偏微分方程式を変数分離して、常微分方程式とします。
　解となる関数は$y=u(x,t)=F(x)G(t)$のように$F(x)$と$G(t)$の積とします。ここでは、弦の振動を考え、波の伝わる速度を$v$とすると、波動方程式
``` math
\begin{align}
\frac{\partial^2 y}{\partial x^2}=\frac{1}{v^2}\frac{\partial^2 y}{\partial t^2}
\end{align}
```
は、
``` math
\frac{d^2 F(x)G(t)}{d x^2}=\frac{1}{v^2}\frac{d^2 F(x)G(t)}{d t^2}
```
となるので、

``` math
G(t)\frac{d^2 F(x)}{d x^2}=\frac{1}{v^2}F(x)\frac{d^2 G(t)}{d t^2}
```
と変形できます。変数$x$の関数と変数$t$の関数でまとめたいので、両辺を$F(x)G(t)$で割ります。
``` math
\frac{1}{F(x)}\frac{d^2 F(x)}{d x^2}=\frac{1}{v^2G(t)}\frac{d^2 G(t)}{d t^2}
```
この式は、変数$x$だけが変化しても、変数$t$だけが変化しても両辺は等しくなることを意味します。この条件を満たすのは、両辺が定数（例えば$K$）となる場合です。
``` math
\frac{1}{F(x)}\frac{d^2 F(x)}{d x^2}=\frac{1}{v^2G(t)}\frac{d^2 G(t)}{d t^2}=K
```
つまり
``` math
\begin{align}
\frac{\partial^2 y}{\partial x^2}=\frac{1}{v^2}\frac{\partial^2 y}{\partial t^2}
\end{align}
```
の偏微分方程式を解くことは以下の、
``` math
\begin{align}
\frac{1}{F(x)}\frac{d^2 F(x)}{dl x^2}=K\\
\frac{1}{v^2G(t)}\frac{d^2 G(t)}{d t^2}=K
\end{align}
```
の常微分方程式を解くことになります。
　$K$は任意の数値ですので、正、0、負の3つの場合に分けて計算します。
　公式集を用いて手計算するのは大変なので、コンピュータのプログラムで計算します。用いるのは、SageMathです。
　https://sagecell.sagemath.org/　
にアクセスして、計算します。
「Type some Sage code below and press Evaluate.」の下の枠に入力し、「Evaluate」を押すと、入力した方程式を解いてくれます。入力の文法はPythonの文法に従います。
以下に、
``` math
\frac{d^2 F(x)}{d x^2}=KF(x)
```
を解く際の入力と得られる出力を示します。
$K=0$の場合
入力
``` 
F=function('F')(x)
desolve( diff(F,x,2)==0,F)
``` 
出力
```
_K2*x + _K1
``` 
この出力は、以下の式を意味します。
``` math
F(x) = k_2x + k_1
```
 $K>0$の場合
$K=k^2>0$として、
``` math
\frac{d^2 F(x)}{d x^2}=k^2F(x)
```
つまり
``` math
\frac{d^2 F(x)}{d x^2}-k^2F(x)=0
```
の方程式を解きます。
入力
``` 
var('k')
F=function('F')(x)
assume(k>0)
desolve( diff(F,x,2)-k*k*F==0, F, ivar=x)
``` 
出力
``` 
_K1*e^(k*x) + _K2*e^(-k*x)
``` 
この出力は、以下の式を意味します。
``` math
F(x) = k_1e^{kx} + k_2e^{-kx}
```

$K<0$の場合
$K=-k^2<0$として、
``` math
\frac{d^2 F(x)}{d x^2}=-k^2F(x)
```
つまり
``` math
\frac{d^2 F(x)}{d x^2}+k^2F(x)=0
```
の方程式を解きます。
入力
``` 
var('k')
F=function('F')(x)
assume(k>0)
desolve( diff(F,x,2)+k*k*F==0, F, ivar=x)
``` 
出力
``` 
_K2*cos(k*x) + _K1*sin(k*x)
``` 
この出力は、以下の式を意味します。
``` math
F(x) = k_2\cos(kx) + k_1\sin(kx)
```

SageMathで波動方程式の一般解が求まりました。
　$G(t)$についても、同様にして解き、この後は、物理的意味のある解かどうかを判断した後、境界条件を使って、特殊解を求めます。

# おわりに
　今回、はじめに、古典的一次元波動方程式を立式する方法を説明しました。次に、SageMathというプログラムを使って、一般解を求めました。境界条件も含めた入力にすると、一気に特殊解が求まると思います。うまく特殊解が得られたら、記事にする予定です。