連立1階常微分方程式をPythonで数値解析

タイトルの通り、SciPyを用いて連立1階常微分方程式をPythonで数値解析してみました。
この手の例としてローレンツ方程式を対象とした記事は他にも調べたら出てくるので、今回はUpper Convected Maxwell Model(以下UCMモデル)の構成方程式を対象として、Pythonで数値解析できればと思います。

そもそもですが、1階常微分方程式をPythonで数値解析するのに参考になったサイトを紹介させていただきます。
https://home.hirosaki-u.ac.jp/jupyter/python-sk/
とても丁寧に解説させているので、Pythonに慣れていない方や常微分方程式をPythonで数値解析するのに手こずっている方はまずこのサイトを参照してみるのが良いのではないでしょうか。

また、UCMモデルとはなんぞやという方は以下のWikipedia等を参照してください。
https://en.wikipedia.org/wiki/Upper-convected_Maxwell_model

以下のコードは

1. 構成方程式で, λ=G=1
2. 剪断速度関数を$\dot\gamma(t) = \gamma\omega cos(\omega t), \gamma=2$

として実装しています。

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# UCMモデルの右辺
def upper_convected_maxwell_stress(t, y, shear_rate):
    G = 1.0
    λ = 1.0

    # 剪断速度関数
    gamma_dot = shear_rate(t)

    #UCMモデルの微分方程式
    dtau_xx_dt = 2 * G * gamma_dot * y[1]  - (1 / λ) * y[0]
    dtau_xy_dt = G * gamma_dot - (1 / λ) * y[1]
    dtau_yy_dt = - (1 / λ) * y[2]

    return [dtau_xx_dt, dtau_xy_dt, dtau_yy_dt]

# 剪断速度関数
def shear_rate(t, omega=1):
    return 2 * omega * np.cos(omega * t)

# 初期条件
tau0_xx = 0.0
tau0_xy = 0.0
tau0_yy = 0.0
initial_conditions = [tau0_xx, tau0_xy, tau0_yy]

# 時間の範囲
t_span = (0, 100)

# 解く
solution = solve_ivp(
    upper_convected_maxwell_stress,
    t_span,
    initial_conditions,
    args=(shear_rate,),
    method='RK45',
    t_eval=np.linspace(0, 100, 1000)
)

# 結果のプロット
plt.plot(solution.t, shear_rate(solution.t), label='Shear Rate', color='black')
plt.plot(solution.t, solution.y[0], label='Shear Stress (τxx)', color='red')
plt.plot(solution.t, solution.y[1], label='Shear Stress (τxy)', color='green')
plt.plot(solution.t, solution.y[2], label='Shear Stress (τyy)', color='blue')
plt.xlabel('Time')
plt.ylabel('Shear Stress')
plt.legend()
plt.show()

この実行結果は以下になります。

今回はUCMモデルを対象としましたが、他にも様々な微分方程式を立式して数値解析することが可能です。
そこまで複雑なことをしていないので細かい解説は省きましたが、これから授業の課題や先行論文の再現等で同様の作業が必要になった方に参考になれば幸いです。