15日前公開・15日前更新・0 pv・20 min read

微分方程式と数値解析の基本③ [強制振動・減衰振動]

数学Python数値解析SciPy

https://cdn.magicode.io/media/notebox/7c75f033-565a-400c-98f3-8a48c20dfa07.jpeg

微分方程式の数値解析の基本②の単振動から発展して減衰を考慮した強制振動を対象にする。

強制振動の微分方程式と理論解

対象の微分方程式

質量

m

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

k

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

バネの左端を原点

x=0

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

x=x_e

とし、

x=x_0

の位置でおもりから手を離す。

その後、バネの左端には

a\cos\omega t

の振動（

a

：振幅、

\omega

：角振動数）を強制的に与え続ける。

任意の時刻

t

で自然長の位置が元の

x_e

から

a\cos\omega t

だけずれていると考えれば、この時のバネの伸びは

x - (x_e +a\cos \omega t)

になる。

さらに、おもりには速度に比例する抵抗力が速度と逆向きに働くものとし、その抵抗係数を

b

（正の値）とする。

それ以外の摩擦力などは全無視。この時の運動方程式は

m\ddot{x} = -k(x - x_e - a\cos\omega t) - b\dot{x}

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

\ddot{x} + \frac{b}{m} \dot{x} + \frac{k}{m} (x-x_e) = \frac{ka}{m} \cos\omega t

ただし、

x

の時間微分を

\dot{x},\ddot{x}

で表している。

これをもう少し扱いやすくするために

X=x-x_e

、

ω_0=\sqrt{\frac{k}{m}}

、

\gamma = \frac{b}{2m}

、

ka=F

とそれぞれ置き換えると

\ddot{X} + 2\gamma \dot{X} + \omega _0^2 X = \frac{F}{m} \cos\omega t \qquad \cdots (1)

となるのでこの(1)を解くことになる。

しかし、(1)は

\frac{F}{m} \cos\omega t

という非同次（非斉次）項、つまり

X

やその微分に関わらない項が入っている非同次（非斉次）線形微分方程式と言われるもので、単振動のような同次線形微分方程式とは異なる解き方が必要になる。

非同次線形微分方程式の一般解は、その方程式の非同次項を無視した同次線形微分方程式の一般解と、非同次線形微分方程式を満たす特殊解を何かしら1つ見つけてそれらを重ね合わせたものになる。

つまり、今回の場合で言えば(1)の右辺を無視した同次線形微分方程式

\ddot{X} + 2\gamma \dot{X} + \omega _0^2 X = 0 \qquad \cdots (2)

の一般解

X_h

と、(1)を満たす特殊解

X_s

を見つけ、それらを合わせた

X_h+X_s

が微分方程式(1)の一般解になる。

同次線形微分方程式の一般解

まずは同次線形微分方程式(2)の一般解

X_h

を求めよう。

この形の微分方程式を解く常套手段で、

X_h=e^{\lambda t}

になるだろうと仮定して(2)に代入すると

\lambda ^2 e^{\lambda t} + 2\gamma \lambda e^{\lambda t} + \omega_0^2 e^{\lambda t} = 0

e^{\lambda t} \not = 0

なので

\lambda ^2 + 2\gamma \lambda + \omega_0^2 = 0

これは

\lambda

に関する二次方程式であるから解の公式で

\lambda

を求められる。

最初に

\gamma = \frac{b}{2m}

とわざわざ1/2を付けて置き換えたのは

\lambda

をちょっとだけ見やすくするためです。

\lambda = -\gamma \pm \sqrt{\gamma^2 - \omega_0^2}

ここで注意しなければならないのがルートの中身が負、0、正のどれになるかで場合分けが必要になることである。

中身が負、つまり

\gamma^2 < \omega_0^2

の場合は

\lambda

は虚数解となり、(2)の解は減衰振動を表す。

中身が0、つまり

\gamma^2 = \omega_0^2

の場合は

\lambda

は実数の重解となり、(2)の解は臨界減衰を表す。

中身が正、つまり

\gamma^2 > \omega_0^2

の場合は

\lambda

は2つの実数解となり、(2)の解は過減衰を表す。

減衰振動はおもりが振動しながら徐々に速度が0に収束するが、臨界減衰と過減衰はいずれも振動することなく指数関数的に速度が0に収束する。

各パターンで得られる解の形が変わってくるので本当はそれぞれ求めないといけないのだが、今回は減衰振動になる場合に限定する。

この時

\lambda = -\gamma \pm i\sqrt{\omega_0^2 - \gamma^2}

となり、

\alpha=\sqrt{\omega_0^2 - \gamma^2}

と置いてやれば

\lambda = -\gamma \pm i\alpha =\lambda_{\pm}

と

\lambda

が2つ求まるので、(2)の解を

X_1

X_2

とすると以下のようになる。

\begin{aligned} X_1 &=e^{\lambda_+ t} = e^{(-\gamma +i\alpha )t} = e^{-\gamma t} e^{i\alpha t} = e^{-\gamma t} (\cos \alpha t + i \sin \alpha t) \\[0.5em] X_2 &=e^{\lambda_- t} = e^{(-\gamma -i\alpha )t} = e^{-\gamma t} e^{-i\alpha t} = e^{-\gamma t} (\cos \alpha t - i \sin \alpha t) \end{aligned}

ただし、途中でオイラーの公式

e^{i\theta} = \cos \theta + i\sin \theta

を使っている。

したがって、任意定数A, Bを用いると(2)の一般解

X_h

は以下のように表せる。（ここの式展開の詳細は補足に記載）

\color{red} X_h = A\frac{X_1+X_2}{2} + B\frac{X_1-X_2}{2i} = e^{-\gamma t} (A\cos \alpha t + B\sin \alpha t) \qquad \cdots (3)

e^{-\gamma t}

は減衰を表し、

A\cos \alpha t + B\sin \alpha t

は単振動と同じ形である。

非同次線形微分方程式の特殊解と一般解

次に非同次線形微分方程式(1)を満たす特殊解

X_s

をなんとか1個でいいので求めたい。

とは言ってもこのままの形で

X_s

を直接予想して確かめるのは非常に困難なので、これまたトリッキーな常套手段で求めていくことになる。

まず、(1)の

X

を

Y

に置き換え（単に

X

と区別するため）右辺の外力を

\frac{F}{m} \sin \omega t

に変えた式

\ddot{Y} + 2\gamma \dot{Y} + \omega _0^2 Y = \frac{F}{m} \sin\omega t \qquad \cdots (1)^{\prime}

を考える。

(1)’の両辺に虚数単位

i

をかけると

i\ddot{Y} + 2i\gamma \dot{Y} + i\omega _0^2 Y = \frac{F}{m} i\sin\omega t \qquad \cdots (1)^{\prime \prime}

(1)’ + (1)’’として整理すると

(\ddot{X} + i\ddot{Y}) + 2\gamma (\dot{X} + i\dot{Y}) + \omega _0^2 (X+iY) = \frac{F}{m} (\cos\omega t + i\sin\omega t)

ここで

Z=X+iY

として右辺はオイラーの公式を使えば上の式は以下のようにまとめられる。

\ddot{Z} + 2\gamma \dot{Z} + \omega_0^2 Z = \frac{F}{m} e^{i\omega t} \qquad \cdots (1)^{\prime \prime \prime}

で、この微分方程式(1)’’’を満たす複素数解

Z

を求めてその実部を取ると、それが元の微分方程式(1)の解とみなせるというのが定石らしい。

ちなみに、(1) で外力として与える振動が

\sin \omega t

なら

Z

の虚部を取ればいい。

それでは(1)’’’の解を求めるため、その解を

Z=Ce^{i\omega t}

と仮定して(1)’’’に代入すると（Cは任意定数）

-\omega^2 Ce^{i\omega t} + 2i\gamma C\omega e^{i\omega t} + C \omega_0^2 e^{i\omega t} = \frac{F}{m} e^{i\omega t}

e^{i\omega t} \not =0

なので、辺々割ってCを求めると

C=\frac{F}{m} \frac{1}{\omega_0^2 - \omega^2 + 2\gamma \omega i} \qquad \cdots (4)

(4)の分母分子に

\omega_0^2 - \omega^2 - 2\gamma \omega i

（分母の共役複素数）をかけると

C=\frac{F}{m} \frac{\omega_0^2 - \omega^2 - 2\gamma \omega i}{(\omega_0^2 - \omega^2)^2 + 4\gamma^2 \omega^2}

さらにもう一工夫として分子にある

\omega_0^2 - \omega^2 - 2\gamma \omega i

について複素平面上で考えれば

\begin{aligned} \cos \phi &= \frac{\omega_0^2 - \omega^2}{\sqrt{(\omega_0^2 - \omega^2)^2 + 4\gamma^2 \omega^2}} \\[0.5em] \sin \phi &= \frac{-2\gamma \omega}{\sqrt{(\omega_0^2 - \omega^2)^2 + 4\gamma^2 \omega^2}} \end{aligned}

と表せるので、結果的に

\begin{aligned} \omega_0^2 - \omega^2 - 2\gamma \omega i &= \sqrt{(\omega_0^2 - \omega^2)^2 + 4\gamma^2 \omega^2} (\cos \phi + i\sin \phi) \\[0.5em] &= \sqrt{(\omega_0^2 - \omega^2)^2 + 4\gamma^2 \omega^2} e^{i\phi} \end{aligned}

とまとめられる。

\phi

は複素数

\omega_0^2 - \omega^2 - 2\gamma \omega i

の偏角である。

よって、最終的なCの形は

C = \frac{F}{m} \frac{e^{i\phi}}{\sqrt{(\omega_0^2 - \omega^2)^2 + 4\gamma^2 \omega^2}}

であるから、(1)’’’の解

Z

は

Z = \frac{F}{m} \frac{e^{i(\omega t +\phi)}}{\sqrt{(\omega_0^2 - \omega^2)^2 + 4\gamma^2 \omega^2}} \qquad \cdots (5)

そして、前述したように(1)の特殊解

X_s

は(5)の実部を取ればいいので

\color{red} X_s = \frac{F}{m} \frac{\cos (\omega t + \phi)}{\sqrt{(\omega_0^2 - \omega^2)^2 + 4\gamma^2 \omega^2}} \qquad \cdots (6)

となる。これで(3)(6)の通り

X_h

X_s

が求まったので、(1)の一般解は

X=X_h+X_s=e^{-\gamma t} (A\cos \alpha t + B\sin \alpha t) + \frac{F}{m} \frac{\cos (\omega t + \phi)}{\sqrt{(\omega_0^2 - \omega^2)^2 + 4\gamma^2 \omega^2}}

であり、

X=x-x_e

と置き換えていたので最終的な

x

の一般解は以下のようになる。

\color{red} x=e^{-\gamma t} (A\cos \alpha t + B\sin \alpha t) + \frac{F}{m} \frac{\cos (\omega t + \phi)}{\sqrt{(\omega_0^2 - \omega^2)^2 + 4\gamma^2 \omega^2}} + x_e \qquad \cdots (7)

初期値問題を解く

(7)は一般解なので今回想定する初期条件を代入して定数A, Bを求める。

(7)を時間微分して速度

\dot{x}

を計算すると

\dot{x} = -e^{\gamma t} \lbrace (A\gamma - B\alpha)\cos \alpha t + (A\alpha + B\gamma)\sin \alpha t \rbrace -\frac{F}{m} \frac{\omega \sin (\omega t + \phi)}{\sqrt{(\omega_0^2 - \omega^2)^2 + 4\gamma^2 \omega^2}}

なので、初期条件

t=0

で

x(0)=x_0

、

\dot{x}(0)=0

を代入して頑張って計算すればA, Bは以下のように求められる。

\begin{aligned} A &= -\frac{F}{m} \frac{\cos \phi}{\sqrt{(\omega_0^2 - \omega^2)^2 + 4\gamma^2 \omega^2}} + x_0 - x_e \\[0.5em] B &= \frac{F}{m} \frac{\omega \sin \phi - \gamma \cos \phi}{\alpha \sqrt{(\omega_0^2 - \omega^2)^2 + 4\gamma^2 \omega^2}} + \frac{(x_0 - x_e)\gamma}{\alpha} \end{aligned}

最終的にまとめると、今回の初期条件の下で強制振動の位置を表す

x

の解は

x=e^{-\gamma t} (A\cos \alpha t + B\sin \alpha t) + \frac{F}{m} \frac{\cos (\omega t + \phi)}{\sqrt{(\omega_0^2 - \omega^2)^2 + 4\gamma^2 \omega^2}} + x_e

ただし、AとBは上記の形で

\begin{aligned} \omega_0 = \sqrt{\frac{k}{m}} \qquad \gamma = \frac{b}{2m} \qquad F=ka \\[0.5em] \alpha=\sqrt{\omega_0^2 - \gamma^2} \qquad \phi=\tan^{-1} \Big( \frac{-2\gamma \omega}{\omega_0^2 - \omega^2} \Big) \end{aligned}

はい、できましたー…

理論解のプロット

求めた

x

の理論解を実際にプロットして振動の様子を見る。

各種係数は強制振動と減衰の様子がいい感じに見やすくなるように適当に設定している。

係数の計算がこの話の肝かつ煩わしいところなので、最初に計算して使いまわせるようにinit_constant_params()という初期化用の関数にまとめている。

import numpy as np
import matplotlib.pyplot as plt
from numpy import pi, sqrt, sin, cos, atan2, exp


def init_constant_params(x0, amplitude):
    """ 計算で使う各種定数を初期化する
        外力の振幅amplitudeは引数として受け取り色々な強制振動パターンで計算できるようにする
    """
    # 物理量の設定
    xe = 0.3  # バネの自然長[m]
    k = 0.5  # バネ定数[N/m]
    m = 0.5  # おもりの質量[kg]
    b = 0.1  # 減衰係数[kg/s]
    omega = pi * 2 / 3  # 外力の角振動数[rad/s]
    # 計算用の各種係数
    force = k * amplitude  # 外力[N]
    gamma = b / (2 * m)
    omega_0 = sqrt(k / m)  # バネ系の固有振動数
    alpha = sqrt(omega_0**2 - gamma**2)
    phi = atan2(-2 * gamma * omega, omega_0**2 - omega**2)
    denominator = sqrt((omega_0**2 - omega**2) **
                       2 + 4 * gamma**2 * omega**2)
    A = - force / m * cos(phi) / denominator + x0 - xe
    B = force / m * (omega * sin(phi) - gamma * cos(phi)) / \
        (alpha * denominator) + (x0 - xe) * gamma / alpha

    return {
        'equilibrium_length': xe,
        'mass': m,
        'gamma': gamma,
        'omega': omega,
        'alpha': alpha,
        'force': force,
        'phi': phi,
        'denominator': denominator,
        'A': A,
        'B': B,
    }


def calc_theoretical_solution(t, params):
    """ 理論解を計算する
    """
    # 計算に使う定数を設定
    xe = params['equilibrium_length']
    m = params['mass']
    gamma = params['gamma']
    omega = params['omega']
    alpha = params['alpha']
    force = params['force']
    phi = params['phi']
    denominator = params['denominator']
    A = params['A']
    B = params['B']

    # 減衰振動だけを想定した計算であり臨界減衰や過減衰ではないので注意！！！
    x = exp(-gamma * t) * (A * cos(alpha * t) + B * sin(alpha * t)) \
        + force / m * cos(omega * t + phi) / denominator + xe

    return x


def forced_vibration(t_list, params):
    """ 全時刻の理論解を計算する """
    x_list = [calc_theoretical_solution(t, params) for t in t_list]
    return x_list


def plot_graph(t_list, xu_list, xf_list):
    fig = plt.figure(figsize=(8, 4))
    fig.suptitle('Forced vibration (theoretical solution)')
    plt.plot(t_list, xu_list, label='unforced', color='b', linewidth=1)
    plt.plot(t_list, xf_list, label='forced', color='r', linewidth=1)
    plt.xlabel("time [s]")
    plt.ylabel("x [m]")
    plt.legend(['unforced', 'forced'],
               loc='upper right', bbox_to_anchor=(1, 1))
    plt.grid()
    plt.show()


if __name__ == '__main__':
    t_end = float(input('計算終了時刻[s]:'))
    x0 = float(input('おもりの初期位置 x0[m]: '))
    v0 = 0  # 初速度は0で固定
    h = 0.01  # 時間刻み幅

    # 初期化した各種定数を取得
    params_u = init_constant_params(x0, amplitude=0)  # 強制力なし unforced
    params_f = init_constant_params(x0, amplitude=0.1)  # 強制力あり forced

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

    # 振動の計算
    xu_list = forced_vibration(t_list, params_u)  # 強制力なし
    xf_list = forced_vibration(t_list, params_f)  # 強制力あり

    plot_graph(t_list, xu_list, xf_list)

強制振動があるパターンを赤線、ないパターン（強制振動の振幅が0）を青線でプロット。

odeintを使った数値解

久々に強制振動の理論解の導出をやってみたが結構しんどい。

しかしこれでも微分方程式の初歩に過ぎない…

で、↑みたいなめんどくさい導出なんかやりたくねんだ！という場合は数値的に解こう。

「微分方程式の数値解析の基本①」「微分方程式の数値解析の基本②」でやったのと同様に自分で解析部分のプログラムを作成するのでもいいが、今回は常微分方程式の数値解を計算できるSciPyのodeintを使ってみる。

odeintとは

SciPyは様々な数学的な処理を行う関数やアルゴリズムを提供するPythonのライブラリで、その中に常微分方程式を数値的に解くことができるodeintという関数がある。

odeint — SciPy v1.14.1 Manual

odeint (func, y0, t, args = (),...)

公式ドキュメントを見ると指定可能な引数は他にもたくさんあるが基本はこれだけでOK。

func

計算対象の微分方程式を定義しておく関数。

連立微分方程式の場合はそれらをすべてこの中で定義しておく。

関数定義のところでは引数を指定しておき、前ステップの計算結果yのリスト、時刻、その他のパラメータ（argsが指定されていれば）などを指定する。

計算対象となる状態変数の初期値のリスト。

連立微分方程式の場合はそれぞれの状態変数の初期値をリスト形式で渡す。

計算の各時刻（時間ステップ）のリスト。

あらかじめnumpyのarange()やlinspace()などを使って生成した時刻のリストを渡すとよい。

args

y0、t以外で何かfuncに渡して使いたいパラメータがあればタプル形式で指定する。

特に何もなければargsは不要。

odeintを使う準備

解きたい微分方程式は2階の微分方程式なので、odeint()で計算するためには媒介変数を使って2本の連立1階微分方程式に分ける必要がある。

\ddot{x} + \frac{b}{m} \dot{x} + \frac{k}{m} (x-x_e) = \frac{ka}{m} \cos\omega t

この微分方程式について媒介変数

v

（速度）を使えば、

\dot{x}=v

、

\ddot{x}=\dot{v}

であるから以下のように分けることができる。

\begin{cases} \dot{x} = v \\ \dot{v} = \frac{1}{m} \lbrace ka\cos \omega t -bv -k(x-x_e) \rbrace \end{cases}

ベクトルの成分表記っぽくすれば

\dfrac{d}{dt} \begin{pmatrix} x \\ v \end{pmatrix} = \begin{pmatrix} v \\ \frac{1}{m} \lbrace ka\cos \omega t -bv -k(x-x_e) \rbrace \end{pmatrix}

これらの式をfuncの中で定義してodeint()に渡してあげればよい。

そして、計算結果であるodeint()の戻り値の中身はNumpy配列で、各時刻の

x

と

v

の結果の配列になっている。

これをそのままmatplotlibのプロット関数に渡すと

x

と

v

両方の結果がプロットされるが、

# 各時刻の結果が[x, v]の順で格納されている
result[:, 0] # x
result[:, 1] # v

のようにすれば片方だけの結果を抽出できる。

数値解のプロット

理論解の時と同じ条件で数値計算を行う。

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from numpy import pi, cos


def init_constant_params(amplitude):
    """ 計算で使う各種定数を初期化する
        外力の振幅amplitudeは引数として受け取り色々な強制振動パターンで計算できるようにする
    """
    k = 0.5  # バネ定数[N/m]
    force = k * amplitude  # 外力[N]

    return {
        'equilibrium_length': 0.3,  # バネの自然長[m]
        'spring_constant': k,
        'mass': 0.5,  # おもりの質量[kg]
        'damping_coefficient': 0.1,  # 減衰係数[kg/s]
        'force': force,  # 外力[N]
        'omega_f': pi * 2 / 3,  # 外力の角振動数[rad/s]
    }


def func(xv, t, params):
    # 計算に使う定数を設定
    xe = params['equilibrium_length']
    k = params['spring_constant']
    m = params['mass']
    b = params['damping_coefficient']
    force = params['force']
    omega_f = params['omega_f']

    x, v = xv
    dx_dt = v
    dv_dt = (force * cos(omega_f * t) - b * v - k * (x - xe)) / m

    return [dx_dt, dv_dt]


def plot_graph(t_list, xu_list, xf_list):
    fig = plt.figure(figsize=(8, 4))
    fig.suptitle('Forced vibration (numerical solution)')
    plt.plot(t_list, xu_list, label='unforced', color='b', linewidth=1)
    plt.plot(t_list, xf_list, label='forced', color='r', linewidth=1)
    plt.xlabel("time [s]")
    plt.ylabel("x [m]")
    plt.legend(['unforced', 'forced'],
               loc='upper right', bbox_to_anchor=(1, 1))
    plt.grid()
    plt.show()


if __name__ == '__main__':
    t_end = float(input('計算終了時刻[s]:'))
    x0 = float(input('初期位置[m]:'))
    v0 = 0  # おもりの初速度[m/s]
    h = 0.01  # 時間刻み幅

    # 初期化した各種定数を取得
    params_u = init_constant_params(amplitude=0)  # 強制力なし unforced
    params_f = init_constant_params(amplitude=0.1)  # 強制力あり forced

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

    # 振動の計算
    result_u = odeint(func, [x0, v0], t_list, args=(params_u,))  # 強制力なし
    result_f = odeint(func, [x0, v0], t_list, args=(params_f,))  # 強制力あり

    # xのみプロットしたいので抽出
    xu_list = result_u[:, 0]
    xf_list = result_f[:, 0]

    plot_graph(t_list, xu_list, xf_list)

本当はちゃんと誤差などを出して定量的に評価しないといけないけれど、結果は理論解の時とほぼ同じなのでヨシ。

補足

同次線形微分方程式の解の変形

\begin{aligned} X_1 &= e^{-\gamma t} (\cos \alpha t + i \sin \alpha t) \\[0.5em] X_2 &=e^{-\gamma t} (\cos \alpha t - i \sin \alpha t) \end{aligned}

X_1,X_2

は共役複素数なので

\begin{aligned} \frac{X_1 + X_2}{2} &= e^{-\gamma t} \cos \alpha t \\[0.5em] \frac{X_1 - X_2}{2i} &=e^{-\gamma t} \sin \alpha t \end{aligned}

という計算でそれぞれの実部と虚部を取り出すことができる。

X_1, X_2

は線形微分方程式(2)の解なので、それらを定数倍して足し合わせた（線形結合した）

\frac{X_1 + X_2}{2},\frac{X_1 - X_2}{2i}

もそれぞれ(2)の解になる。さらにそれらを線型結合した

A \frac{X_1 + X_2}{2} + B\frac{X_1 - X_2}{2i} = e^{-\gamma t} (A\cos \alpha t + B\sin \alpha t)

も(2)の解になるので、結局

X_h = e^{-\gamma t} (A\cos \alpha t + B\sin \alpha t)

のように式変形できる。

変位

X_h

は物理量なので実数でなければいけないから、この変形により虚数単位

i

を使わずに解を表現することができる。

数値計算での減衰の場合分けについて

理論計算だと

\gamma

と

\omega_0

の大小関係に応じて減衰振動、過減衰、臨界減衰の式を切り替えられるプログラムにする必要があるが、数値計算はそんなことを考えなくても条件に応じた答えを出してくれる。

例えば

減衰振動：

\gamma < \omega_0

つまり

b^2<4mk

臨界減衰：

\gamma = \omega_0

つまり

b^2=4mk

過減衰：

\gamma > \omega_0

つまり

b^2>4mk

となるように係数を定めると、場合分けしなくてもそれぞれの条件に応じた運動を計算して結果を出してくれる。（下の図は強制振動なしで計算した結果）

微分方程式と数値解析の基本③ [強制振動・減衰振動]

強制振動の微分方程式と理論解

対象の微分方程式

同次線形微分方程式の一般解

非同次線形微分方程式の特殊解と一般解

初期値問題を解く

理論解のプロット

odeintを使った数値解

odeintとは

odeintを使う準備

数値解のプロット

補足

同次線形微分方程式の解の変形

数値計算での減衰の場合分けについて

Discussion

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