14日前公開・14日前更新・0 pv・10 min read

微分方程式と数値解析の基本②

数学Python数値解析

https://cdn.magicode.io/media/notebox/blob_XEAGpVz

微分方程式と数値解析の基本①ではとても単純な自由落下で理論解と数値解を計算したが、今回はもうちょっとだけ複雑になった単振動を対象にする。

単振動を表す微分方程式と理論解

水平バネによる単振動

質量

m

のおもりが付いているバネ（バネ定数

k

）が水平に置かれており一端が固定された状態を考える。

バネが固定されている位置を原点

x=0

、バネが自然長の時のおもりの位置を

x=x_e

とし、

x=x_0

の位置で静かに（=初速度0）おもりから手を離す。

ただし、おもりと床の間の摩擦力や空気抵抗などは無視し、運動は

x

軸方向のみ考えるものとする。

この時の運動方程式は

m \frac{d^2x}{dt^2} = -k(x - x_e)

となるので、解析対象の微分方程式は

\frac{d^2x}{dt^2} = -\frac{k}{m}(x - x_e)\dots(☆)

ここで

X(t) = x(t) - x_e

と置くと

\frac{dX}{dt}=\frac{dx}{dt},\space \frac{d^2X}{dt^2}=\frac{d^2x}{dt^2}

であるから、☆をXの微分方程式に書き換えることができる。

\frac{d^2X}{dt^2} = -\frac{k}{m}X\dots(☆☆)

☆☆のような微分方程式は、その形から解が三角関数や指数関数になると推定して解く方法などがあるが、ここでは三角関数を使う方法で解いてみる。

☆☆は

X

を

t

で2回微分すると

X

自身の

-\frac{k}{m}

倍となることを表している。

ここで三角関数

\sin, \cos

を2回微分すると

\frac{d^2\sin\omega t}{dt^2} = -\omega^2\sin\omega t, \space \frac{d^2\cos\omega t}{dt^2} = -\omega^2\cos\omega t

のように元の三角関数の

-\omega^2

倍になることを利用し、

\omega=\sqrt{\frac{k}{m}}

と置けば

X_1 = \sin\omega t

と

X_2 = \cos\omega t

は☆☆の解になる。

したがって、☆☆のような線形同次常微分方程式については解の重ね合わせが成り立つので、一般解は

X = AX_1 + BX_2 = A\sin\omega t + B\cos\omega t （A, Bは定数）

つまり☆の一般解は

\begin{aligned} x &= A\sin\omega t + B\cos\omega t + x_e \\[0.5em] \frac{dx}{dt} &= A\omega\cos\omega t - B\omega\sin\omega t \end{aligned}

後は初期値を代入して定数A, Bを求めればよく

\begin{aligned} x(0) &= B + x_e = x_0 \to B = x_0 - x_e \\[0.5em] \left. \frac{dx}{dt}\right| _{t=0} &= A\omega = 0 \to A = 0 \end{aligned}

以上より水平バネによる単振動の理論解は

\begin{aligned} x(t) &= (x_0 - x_e)\cos\omega t + x_0 \\[0.5em] v(t) &= -(x_0 - x_e)\omega\sin\omega t \end{aligned}

鉛直バネによる単振動

水平の時と同じバネとおもりを天井からぶら下げた状態を考える。

鉛直下向きを正として

y

軸を取り、バネが固定されている位置を原点

y=0

、バネが自然長の時のおもりの位置を

y=y_e

とし、

y=y_0

の位置で静かに（=初速度0）おもりから手を離す。

運動は

y

軸方向のみ考えるものとする。

なお、鉛直の場合は重力の影響がある分バネの自然長の位置よりも少し伸びた位置でおもりが静止するので、自然長からの伸びの長さを

b

とすると力の釣合いの条件から

kb = mg \to b = \frac{mg}{k}

というように

b

を計算できる。

この鉛直バネの単振動の運動方程式は

m \frac{d^2y}{dt^2} = -k(y - y_e)+mg

となるので、解析対象の微分方程式は

\frac{d^2y}{dt^2} = -\frac{k}{m}(y - y_e)+g\dots(★)

後は水平バネの時と同じように★の右辺の置き換えを行って解いていけばいいが、仮に

Y(t) = y(t)-y_e

と置くと

\frac{d^2Y}{dt^2} = -\frac{k}{m}Y+g

となり非同次項

g

（

Y

や

Y

の微分に関わらない項）が出てきてしまい後の解法に一手間加える必要があるのでちょっと面倒。

なので、★の右辺全体を

-\frac{k}{m}

で括り出す変形をし

\frac{d^2y}{dt^2} = -\frac{k}{m}(y - y_e - \frac{mg}{k})

として

Y(t) = y(t) - y_e - \frac{mg}{k}

と置いてやれば

\frac{d^2Y}{dt^2} = -\frac{k}{m}Y\dots(★★)

というように水平バネと同じ形の線形同次常微分方程式にできる。

後は同じ手順で解いてやれば★★の一般解は

Y = A\sin\omega t + B\cos\omega t （A, Bは定数）

よって、★の一般解は

\begin{aligned} y &= A\sin\omega t + B\cos\omega t + y_e + \frac{mg}{k} \\[0.5em] \frac{dy}{dt} &= A\omega\cos\omega t - B\omega\sin\omega t \end{aligned}

初期値を代入して定数A, Bを求めると

\begin{aligned} y(0) &= B + y_e + \frac{mg}{k}= y_0 \to B = y_0 - y_e -\frac {mg}{k} \\[0.5em] \left. \frac{dy}{dt}\right| _{t=0} &= A\omega = 0 よりA = 0 \end{aligned}

以上より鉛直バネによる単振動の理論解は

\begin{aligned} y(t) &= (y_0 - y_e - \frac{mg}{k})\cos\omega t + y_e + \frac{mg}{k} \\[0.5em] v(t) &= -(y_0 - y_e - \frac{mg}{k})\omega \sin\omega t \end{aligned}

数値計算する

計算するための準備など

微分方程式と数値解析の基本①でやったのと同様に4次のルンゲ・クッタ法で解いてみる。

今回も解きたい微分方程式は2階の常微分方程式なので、媒介変数

v

を使って2本の連立1階微分方程式に分ける。

\begin{aligned} \frac{dv}{dt} &= -\frac{k}{m}x = f_v(t, y, v) \\[0.5em] \frac{dy}{dt} &= v = f_y(t, y, v) \end{aligned}

速度と位置の計算手順も特に変わりなく、

k_1

〜

k_4

と

l_1

〜

l_4

を順番に計算して時間ステップ毎の値を計算していく。

前回の記事と異なるのは

\frac{dv}{dt}

の式だけで他は基本的に変わらない。

計算コード

まず水平バネの単振動を計算するためのPythonコードを載せる。

使うバネは適当に自然長0.3mでバネ定数50N/m、おもりは質量0.5kgで固定値とする。

import numpy as np
import matplotlib.pyplot as plt
import math

h = 0.01  # 時間刻み幅
xe = 0.3  # バネの自然長[m]
k = 50  # バネ定数[N/m]
m = 0.5  # おもりの質量[kg]
v0 = 0  # おもりの初速度[m/s]（静かに手を離すので必ず0）


def f_v(t, x, v):
    """ dv/dt """
    return - k / m * (x - xe)


def f_x(t, x, v):
    """ dx/dt """
    return v


def calc_numerical_solution(t, x, v):
    """ 4次のルンゲ・クッタ法により次ステップの値を計算する """
    k1 = f_v(t, x, v) * h
    l1 = f_x(t, x, v) * h
    k2 = f_v(t + h / 2, x + l1 / 2, v + k1 / 2) * h
    l2 = f_x(t + h / 2, x + l1 / 2, v + k1 / 2) * h
    k3 = f_v(t + h / 2, x + l2 / 2, v + k2 / 2) * h
    l3 = f_x(t + h / 2, x + l2 / 2, v + k2 / 2) * h
    k4 = f_v(t + h, x + l3, v + k3) * h
    l4 = f_x(t + h, x + l3, v + k3) * h
    v_next = v + (k1 + 2 * k2 + 2 * k3 + k4) / 6
    x_next = x + (l1 + 2 * l2 + 2 * l3 + l4) / 6

    return x_next, v_next


def calc_theoretical_solution(t):
    """ 理論解を計算する """
    omega = math.sqrt(k / m)
    x = (x0 - xe) * math.cos(omega * t) + xe
    v = -(x0 - xe) * omega * math.sin(omega * t)

    return x, v


def simple_harmonic_motion(x0, t_list):
    """ 単振動の計算 """
    xn = x0
    vn = v0
    # t=0の理論解を計算する
    xt, vt = calc_theoretical_solution(0)
    # 初期条件を解のリストに追加
    xn_list = [xn]  # xの数値解リスト
    vn_list = [vn]  # vの数値解リスト
    xt_list = [xt]  # xの理論解リスト
    vt_list = [vt]  # vの理論解リスト

    # t>0について解の計算を行う
    for t in t_list[1:]:
        # 数値解を計算
        xn, vn = calc_numerical_solution(t, xn, vn)
        xn_list.append(xn)
        vn_list.append(vn)
        # 理論解を計算
        xt, vt = calc_theoretical_solution(t)
        xt_list.append(xt)
        vt_list.append(vt)

    return xn_list, vn_list, xt_list, vt_list


def plot_graph(xn_list, vn_list, xt_list, vt_list):
    fig = plt.figure(figsize=(8, 4))
    fig.suptitle('Horizontal simple harmonic motion')

    # 位置のプロット
    plt.subplot(2, 1, 1)
    plt.plot(t_list, xn_list, label='numerical',
             color='r', linewidth=1)  # 数値解
    plt.plot(t_list, xt_list, label='theoretical', color='g', linestyle='None',
             marker='o', markersize=1)  # 理論解
    plt.hlines(xe, t_list[0], t_list[-1], linestyles='dotted',
               colors='k', linewidth=1)  # 自然長の線
    plt.ylabel("x [m]")
    # 凡例
    plt.legend(['numerical', 'theoretical'],
               loc='upper left', bbox_to_anchor=(1, 1))

    # 速度のプロット
    plt.subplot(2, 1, 2)
    plt.plot(t_list, vn_list, color='r', linewidth=1)  # 数値解
    plt.plot(t_list, vt_list, color='g', linestyle='None',
             marker='o', markersize=1)  # 理論解
    plt.hlines(0, t_list[0], t_list[-1], linestyles='dotted',
               colors='k', linewidth=1)  # 速度0の線
    plt.xlabel("time [s]")
    plt.ylabel("velocity [m/s]")

    fig.tight_layout()
    plt.show()


if __name__ == '__main__':
    x0 = float(input('初期位置[x]:'))
    t_end = float(input('計算終了時刻[s]:'))

    # 計算時間ステップのリストを生成する（最後の時刻t_endも含めたいので+hしている）
    t_list = np.arange(0, t_end + h, h)

    # 単振動の計算
    xn_list, vn_list, xt_list, vt_list = simple_harmonic_motion(x0, t_list)

    plot_graph(xn_list, vn_list, xt_list, vt_list)

鉛直バネの単振動のコードも水平の場合とほぼ同じ（バネ、おもりの条件も同じ）で、

x

座標ではなく

y

座標にしているのと、

f_y(t, y, v)

の式および理論解の式が異なるだけなのでその部分だけ載せる。

import numpy as np
import matplotlib.pyplot as plt
import math

h = 0.01  # 時間刻み幅
ye = 0.3  # バネの自然長[m]
k = 50  # バネ定数[N/m]
m = 0.5  # おもりの質量[kg]
v0 = 0  # おもりの初速度[m/s]（静かに手を離すので必ず0）
g = 9.80665  # 重力加速度[m/s^2]
yb = ye + m * g / k  # おもりの釣合いの位置[m]


def f_v(t, y, v):
    """ dv/dt """
    return g - k / m * (y - ye)


def calc_theoretical_solution(t):
    """ 理論解を計算する """
    omega = math.sqrt(k / m)
    y = (y0 - ye - m * g / k) * math.cos(omega * t) + m * g / k + ye
    v = -(y0 - ye - m * g / k) * omega * math.sin(omega * t)

    return y, v

計算結果のプロット

水平、鉛直ともにおもりの初期位置（手を離す位置）は0.6mとして、0〜2秒まで計算している。

振動の中心はバネの自然長の位置（黒の点線）で振幅0.3m程度の振動になる。

数値解（赤線）は理論解（緑の点）をよく再現できている。

鉛直バネの場合は重力の影響が含まれるので、振動の中心はバネの自然長の位置ではなくおもりの釣合いの位置（黒の点線）になる。

振幅も0.2m程度であり水平の単振動に比べて小さくなっている。