Ccmmutty logo
Commutty IT
0 pv10 min read

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

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

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

水平バネによる単振動

質量mmのおもりが付いているバネ(バネ定数kk)が水平に置かれており一端が固定された状態を考える。
バネが固定されている位置を原点x=0x=0、バネが自然長の時のおもりの位置をx=xex=x_eとし、x=x0x=x_0の位置で静かに(=初速度0)おもりから手を離す。
ただし、おもりと床の間の摩擦力や空気抵抗などは無視し、運動はxx軸方向のみ考えるものとする。
この時の運動方程式は
md2xdt2=k(xxe)m \frac{d^2x}{dt^2} = -k(x - x_e)
となるので、解析対象の微分方程式は
d2xdt2=km(xxe)()\frac{d^2x}{dt^2} = -\frac{k}{m}(x - x_e)\dots(☆)
ここでX(t)=x(t)xeX(t) = x(t) - x_eと置くと
dXdt=dxdt, d2Xdt2=d2xdt2\frac{dX}{dt}=\frac{dx}{dt},\space \frac{d^2X}{dt^2}=\frac{d^2x}{dt^2}
であるから、☆をXの微分方程式に書き換えることができる。
d2Xdt2=kmX(☆☆)\frac{d^2X}{dt^2} = -\frac{k}{m}X\dots(☆☆)
☆☆のような微分方程式は、その形から解が三角関数や指数関数になると推定して解く方法などがあるが、ここでは三角関数を使う方法で解いてみる。
☆☆はXXttで2回微分するとXX自身のkm-\frac{k}{m}倍となることを表している。
ここで三角関数sin,cos\sin, \cosを2回微分すると
d2sinωtdt2=ω2sinωt, d2cosωtdt2=ω2cosωt\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
のように元の三角関数のω2-\omega^2倍になることを利用し、ω=km\omega=\sqrt{\frac{k}{m}}と置けば X1=sinωtX_1 = \sin\omega tX2=cosωtX_2 = \cos\omega tは☆☆の解になる。
したがって、☆☆のような線形同次常微分方程式については解の重ね合わせが成り立つので、一般解は
X=AX1+BX2=Asinωt+BcosωtA,Bは定数)X = AX_1 + BX_2 = A\sin\omega t + B\cos\omega t (A, Bは定数)
つまり☆の一般解は
x=Asinωt+Bcosωt+xedxdt=AωcosωtBωsinωt\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を求めればよく
x(0)=B+xe=x0B=x0xedxdtt=0=Aω=0A=0\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}
以上より水平バネによる単振動の理論解は
x(t)=(x0xe)cosωt+x0v(t)=(x0xe)ωsinωt\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}

鉛直バネによる単振動

水平の時と同じバネとおもりを天井からぶら下げた状態を考える。
鉛直下向きを正としてyy軸を取り、バネが固定されている位置を原点y=0y=0、バネが自然長の時のおもりの位置をy=yey=y_eとし、y=y0y=y_0の位置で静かに(=初速度0)おもりから手を離す。
運動はyy軸方向のみ考えるものとする。
なお、鉛直の場合は重力の影響がある分バネの自然長の位置よりも少し伸びた位置でおもりが静止するので、自然長からの伸びの長さをbbとすると力の釣合いの条件から
kb=mgb=mgkkb = mg \to b = \frac{mg}{k}
というようにbbを計算できる。
この鉛直バネの単振動の運動方程式は
md2ydt2=k(yye)+mgm \frac{d^2y}{dt^2} = -k(y - y_e)+mg
となるので、解析対象の微分方程式は
d2ydt2=km(yye)+g()\frac{d^2y}{dt^2} = -\frac{k}{m}(y - y_e)+g\dots(★)
後は水平バネの時と同じように★の右辺の置き換えを行って解いていけばいいが、仮にY(t)=y(t)yeY(t) = y(t)-y_eと置くと
d2Ydt2=kmY+g\frac{d^2Y}{dt^2} = -\frac{k}{m}Y+g
となり非同次項ggYYYYの微分に関わらない項)が出てきてしまい後の解法に一手間加える必要があるのでちょっと面倒。
なので、★の右辺全体をkm-\frac{k}{m}で括り出す変形をし
d2ydt2=km(yyemgk)\frac{d^2y}{dt^2} = -\frac{k}{m}(y - y_e - \frac{mg}{k})
としてY(t)=y(t)yemgkY(t) = y(t) - y_e - \frac{mg}{k}と置いてやれば
d2Ydt2=kmY(★★)\frac{d^2Y}{dt^2} = -\frac{k}{m}Y\dots(★★)
というように水平バネと同じ形の線形同次常微分方程式にできる。
後は同じ手順で解いてやれば★★の一般解は
Y=Asinωt+BcosωtA,Bは定数)Y = A\sin\omega t + B\cos\omega t (A, Bは定数)
よって、★の一般解は
y=Asinωt+Bcosωt+ye+mgkdydt=AωcosωtBωsinωt\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を求めると
y(0)=B+ye+mgk=y0B=y0yemgkdydtt=0=Aω=0よりA=0\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}
以上より鉛直バネによる単振動の理論解は
y(t)=(y0yemgk)cosωt+ye+mgkv(t)=(y0yemgk)ωsinωt\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階の常微分方程式なので、媒介変数vvを使って2本の連立1階微分方程式に分ける。
dvdt=kmx=fv(t,y,v)dydt=v=fy(t,y,v)\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}
速度と位置の計算手順も特に変わりなく、k1k_1k4k_4l1l_1l4l_4を順番に計算して時間ステップ毎の値を計算していく。
前回の記事と異なるのはdvdt\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)
鉛直バネの単振動のコードも水平の場合とほぼ同じ(バネ、おもりの条件も同じ)で、xx座標ではなくyy座標にしているのと、fy(t,y,v)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程度であり水平の単振動に比べて小さくなっている。

Discussion

コメントにはログインが必要です。