Pythonで学ぶ連成振動の解析的解法と数値的解法

数学・物理学の具体的な計算にPythonを使って、数学・物理学もPythonも同時に学んでしまいましょう。前回の記事

pianofisica.hatenablog.com

では、Pythonを使って振動を表す運動方程式（2階常微分方程式）の解を数値的に求めました。また

pianofisica.hatenablog.com

では、数値解によって記述される質点の運動の様子をアニメーションによって可視化してみました。

今回の記事はこれらの続編・応用で、複数の質点が互いにバネでつながれた、いわゆる連成振動について取り上げてみたいと思います。

バネにつながれた１質点の運動（単振動）
- 解析的解法
- 数値的解法
バネにつながれた２質点の運動（連成振動）
- 解析的解法
- 数値的解法
バネにつながれたN質点の運動

バネにつながれた１質点の運動（単振動）

解析的解法

質量 $m_1$ をもつ質点の平衡点からの変位を $x_1$ とします（図１上部）。質点の左右はそれぞれバネ定数 $k_1,\, k_2$ をもつバネでつながっているものとしましょう。この質点の運動方程式は以下で与えられます：

$\qquad \displaystyle{m_1\frac{d^2x_1}{dt^2}=-(k_1+k_2) x_1}$

運動方程式の一般解は

$\qquad \displaystyle{x_1(t)=A_1 \cos(\omega_{12} t)+B_1\sin(\omega_{12} t)}$

で与えられます。ここで $A_1,\, B_1$ は積分定数で、初期条件によって決定されます。また

$\qquad \displaystyle{\omega_{12}=\sqrt{\frac{k_1+k_2}{m_1}}}$

は振動の各速度を与えます。よって振動の周期は

$\qquad \displaystyle{T_{12}=\frac{2\pi}{\omega_{12}}}$

によって与えられます。

数値的解法

運動の様子を図示するために運動方程式を数値的に解いてアニメーションを作成してみましょう。その方法の詳細については上記の過去記事の解説を参照してください。簡単のために $m_1=1$ 、 $k_1=k_2=1/2$ とし、初期条件として $x_1(0)=1,\,\dot{x}_1(0)=0$ とします。ここでは必要なコードだけ載せておきます：

N = 3000   # ステップ数
h = 0.01   # 差分間隔
T = [0]    # 初期時刻
X1 = [1]   # 初期位置
V1 = [0]   # 初期速度

# 差分方程式
def t(k):
    return h+T[k-1]
def x1(k):
    return X1[k-1]+h*(V1[k-1])+h**2/2*(-X1[k-1])
def v1(k):
    return V1[k-1]+h*(-X1[k-1])+h**2/2*(-V1[k-1])

for k in range(N):
      T.append( t(k+1) )
      X1.append( x1(k+1) )
      V1.append( v1(k+1) )

%matplotlib inline
# JupyterNotebookでアニメーションを表示させるために必要です
import matplotlib.animation as animation

import matplotlib.pyplot as plt
fig = plt.figure()
# グラフの設定
plt.axes().set_aspect('equal')
plt.xlim(-1.5, 1.5)  # x軸の表示範囲
plt.ylim(-0.5, 0.5)  # y軸の表示範囲

# アニメーションのフレームの作成
frames = []
num_frames = int(N/10)
for i in range(num_frames):
    x = X1[10*i]    # アニメーションのためにプロットするデータ点は10個おき
    y = 0
    frame = plt.scatter(x, y, color='blue')
    frames.append([frame])

# アニメーションの描画
ani = animation.ArtistAnimation(fig, frames, interval=50)

# GIFファイルとして保存 & アニメーションの表示
from IPython.display import Image
f = r'harmonic_ani.gif'
writergif = animation.PillowWriter(fps=50)
ani.save(f,writer=writergif)
plt.close()
Image(open('harmonic_ani.gif','rb').read())

以上から次のようなアニメーションが出力されるとともに、ローカルに harmonic_ani.gif というファイル名のGIFファイルが作成されるはずです：

バネにつながれた２質点の運動（連成振動）

解析的解法

問題設定を拡大して、図１下部のような２つの質点がバネでつながっている場合を考えてみましょう。各質点の平衡点からの変位をそれぞれ $x_1,\,x_2$ とすると運動方程式は以下で与えられます：

$\qquad \displaystyle{\begin{aligned}m_1\frac{d^2x_1}{dt^2}&=-k_1 x_1 - k_2(x_1-x_2)\\&=-(k_1+k_2)x_1+k_2x_2 \\ m_2\frac{d^2x_2}{dt^2}&=-k_2 (x_2-x_1) -k_3 x_2\\&=k_2x_1-(k_2+k_3)x_2\end{aligned}}$

よって行列を用いて

$\qquad \displaystyle{\begin{aligned}\frac{d^2}{dt^2}\begin{pmatrix} x_1\\ x_2 \end{pmatrix}&=-K\begin{pmatrix} x_1\\ x_2 \end{pmatrix}\qquad K=\begin{pmatrix}\frac{k_1+k_2}{m_1} & -\frac{k_2}{m_1} \\ -\frac{k_2}{m_2} & \frac{k_2+k_3}{m_2}\end{pmatrix}\end{aligned}}$

と書き表せることがわかります。右辺の行列 $K$ の各成分がすべて定数パラメタで構成されていること、そして適当な正則行列 $\Lambda$ （具体形はあとで述べます）によって行列 $K$ が対角化できることに注意します：

$\qquad \displaystyle{\Lambda^{-1}K\Lambda=\begin{pmatrix}\Omega^2_1 & 0 \\ 0 & \Omega^2_2\end{pmatrix}}$

ここで

$\qquad \displaystyle{ {\rm Tr}\,K=\frac{k_1+k_2}{m_1}+\frac{k_2+k_3}{m_2} \qquad\quad {\rm det}\,K=\frac{k_1k_2+k_2k_3+k_3k_1}{m_1m_2} }$

がいずれも正の値であることから、行列 $K$ の固有値は正定値であり、したがってそれらの正の根号がとれることを使いました。２次の単位行列を $I_2$ として

$\qquad \displaystyle{\Lambda\Lambda^{-1}=\Lambda^{-1}\Lambda=I_2}$

であること、そして $\Lambda$ が定数行列であることから $\Lambda^{-1}$ が時間微分をすり抜けることに注意して、行列表記した運動方程式の両辺に $\Lambda^{-1}$ を掛けてみます：

$\qquad \displaystyle{\begin{aligned}\frac{d^2}{dt^2}\Lambda^{-1}\begin{pmatrix} x_1\\ x_2 \end{pmatrix}&=-\Lambda^{-1}K\begin{pmatrix} x_1\\ x_2 \end{pmatrix}\\&=-\Lambda^{-1}K\Lambda\Lambda^{-1}\begin{pmatrix} x_1\\ x_2 \end{pmatrix}\\&=-\begin{pmatrix}\Omega^2_1 & 0 \\ 0 & \Omega^2_2\end{pmatrix}\Lambda^{-1}\begin{pmatrix} x_1\\ x_2 \end{pmatrix}\end{aligned}}$

となることがわかります。さらに

$\qquad \displaystyle{\begin{pmatrix} X_1\\ X_2 \end{pmatrix}=\Lambda^{-1}\begin{pmatrix} x_1\\ x_2 \end{pmatrix}}$

として新たな変数 $X_1,\, X_2$ を導入すると、運動方程式が２つの独立な単振動の式

$\qquad \displaystyle{\frac{d^2X_i}{dt^2}=-\Omega_i^2 X_i \qquad(i=1,2)}$

に帰着することがわかります。その一般解は

$\qquad \displaystyle{X_i(t)=A_i\cos(\Omega_it)+B_i\sin(\Omega_it) \qquad(i=1,2)}$

で与えられますから、最終的に求めていた解は

$\qquad \displaystyle{\begin{pmatrix} x_1(t)\\ x_2(t) \end{pmatrix}=\Lambda\begin{pmatrix} A_1\cos(\Omega_1t)+B_1\sin(\Omega_1t)\\ A_2\cos(\Omega_2t)+B_2\sin(\Omega_2t) \end{pmatrix}}$

となります。

一応、対角化するための正則行列 $\Lambda$ を求めておきましょう：

import sympy as sp
sp.init_printing()
sp.var('k1, k2, k3, m1, m2', positive=True)
K = sp.Matrix([[(k1+k2)/m1, -k2/m1],[-k2/m2, (k2+k3)/m2]])
L, W = K.diagonalize()
L

から、行列 $\Lambda$ の具体形が出力されます。そして対角化された行列 ${\rm diag}(\Omega_1^2,\Omega_2^2)$ が

によって出力されます。ただいずれも、３つのバネ定数、２つの質量、とパラメタの数が多いため各成分はこれらの複雑な式になります。 $\Lambda^{-1}K\Lambda={\rm diag}(\Omega_1^2,\Omega_2^2)$ となることを確認してみましょう：

sp.cancel(sp.simplify( L**(-1)*K*L - W ))

から、確認できます。単純な場合として

$\qquad \displaystyle{k_1=k_2=k_3, \qquad m_1=m_2}$

としてみます：

sp.cancel(L.subs(m2, m1).subs(k2,k1).subs(k3,k1))

および

sp.cancel(W.subs(m2, m1).subs(k2,k1).subs(k3,k1))

から

$\qquad \displaystyle{\Lambda=\begin{pmatrix} 1 & -1 \\ 1 & 1 \end{pmatrix} \qquad {\rm diag}(\Omega_1^2,\Omega_2^2)=\begin{pmatrix} \frac{k_1}{m_1} & 0 \\ 0 & \frac{3k_1}{m_1} \end{pmatrix} }$

であることがわかります。

数値的解法

１質点の場合と同様に数値的に解いてみます。運動方程式を１階化します：

$\qquad \displaystyle{\begin{aligned} \frac{dx_1}{dt}&=\frac{p_1}{m_1}=v_1 \qquad \frac{dv_1}{dt}=-\frac{k_1+k_2}{m_1}x_1+\frac{k_2}{m_1}x_2 \\ \frac{dx_2}{dt}&=\frac{p_2}{m_2}=v_2 \qquad \frac{dv_2}{dt}=\frac{k_2}{m_2}x_1-\frac{k_2+k_3}{m_2}x_2\end{aligned}}$

これらから展開式

$\qquad \displaystyle{\begin{aligned} x_1(t+h)&=x_1(t)+h\frac{dx_1}{dt}(t)+\frac{h^2}{2}\frac{d^2x_1}{dt^2}(t)+\mathcal{O}(h^3)\\ &=x_1(t)+hv_1(t)+\frac{h^2}{2}\left(-\frac{k_1+k_2}{m_1}x_1(t)+\frac{k_2}{m_1}x_2(t)\right)+\mathcal{O}(h^3) \\ v_1(t+h)&=v_1(t)+h\frac{dv_1}{dt}(t)+\frac{h^2}{2}\frac{d^2v_1}{dt^2}(t)+\mathcal{O}(h^3) \\ &=v_1(t)+h\left(-\frac{k_1+k_2}{m_1}x_1(t)+\frac{k_2}{m_1}x_2(t)\right)\\ & \qquad+\frac{h^2}{2}\left(-\frac{k_1+k_2}{m_1}v_1(t)+\frac{k_2}{m_1}v_2(t)\right)+\mathcal{O}(h^3) \end{aligned}}$

を得ます。 $x_2,\, v_2$ についても同様です。ここから差分方程式が作れて、以下のコードのようにして数値計算できます：

N = 6000   # ステップ数
h = 0.01   # 差分間隔
T = [0]    # 初期時刻
X1 = [-0.8]   # 質点１の初期位置
V1 = [-0.3]   # 質点１の初期速度
X2 = [0.1]    # 質点２の初期位置
V2 = [-0.4]   # 質点２の初期速度

# パラメタ値の設定
m1 = 1.0
m2 = 6.0
k1 = 1.0
k2 = 2.0
k3 = 5.0

# 差分方程式
def t(k):
    return h+T[k-1]
def x1(k):
    return X1[k-1]+h*(V1[k-1])+h**2/2*( -(k1+k2)/m1*X1[k-1] + k2/m1*X2[k-1] )
def v1(k):
    return V1[k-1]+h*( -(k1+k2)/m1*X1[k-1] + k2/m1*X2[k-1] )\
        +h**2/2*( -(k1+k2)/m1*V1[k-1] + k2/m1*V2[k-1] )
def x2(k):
    return X2[k-1]+h*(V2[k-1])+h**2/2*( k2/m2*X1[k-1] - (k2+k3)/m2*X2[k-1] )
def v2(k):
    return V2[k-1]+h*( k2/m2*X1[k-1] - (k2+k3)/m2*X2[k-1] )\
        +h**2/2*( k2/m2*V1[k-1] - (k2+k3)/m2*V2[k-1] )

for k in range(N):
      T.append( t(k+1) )
      X1.append( x1(k+1) )
      V1.append( v1(k+1) )
      X2.append( x2(k+1) )
      V2.append( v2(k+1) )

%matplotlib inline
# JupyterNotebookでアニメーションを表示させるために必要です
import matplotlib.animation as animation

import matplotlib.pyplot as plt
fig = plt.figure()
# グラフの設定
plt.axes().set_aspect('equal')
plt.xlim(-3.0, 3.0)  # x軸の表示範囲
plt.ylim(-0.5, 0.5)  # y軸の表示範囲

# アニメーションのフレームの作成
frames = []
num_frames = int(N/15)
for i in range(num_frames):
    y = 0
    x1 = X1[15*i] -1.0   # アニメーションのためにプロットするデータ点は15個おき
    x2 = X2[15*i] +1.0
    frame1 = plt.scatter(x1, y, color='blue')
    frame2 = plt.scatter(x2, y, color='orange')
    frames.append([frame1,frame2])

# アニメーションの描画
ani = animation.ArtistAnimation(fig, frames, interval=30)

# GIFファイルとして保存 & アニメーションの表示
from IPython.display import Image
f = r'two_harmonics_ani.gif'
writergif = animation.PillowWriter(fps=30)
ani.save(f,writer=writergif)
plt.close()
Image(open('two_harmonics_ani.gif','rb').read())