微分方程式を解く

2021年11月30日

はじめに

SciPy の odeint モジュールを用いて微分方程式を解く。

自由落下

自由落下問題を解く。

モジュールの用意。

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

パラメータ。

g = 9.81

関数の用意。

def func(s, t, m):
    x, v = s
    f = -m*g
    dsdt = [v, f/m]
    return dsdt 

解く方程式は、力 f、質量 m、重力加速度 g として、f = -mg である。関数の引数は、s に変数リスト、t に時刻、その他はパラメータとする。s はここでは位置 x、x の一回微分 (v) のリストである。返り値は、それらの一回微分のリストである。この場合、速度 v と加速度 f/m である。

実行。

m = 1
x0 = [5, 0]
t = np.linspace(0, 1, 101)
sol = odeint(func, x0, t, args=(m,))

odeint に関数を渡して方程式を解いている。x0 は初期値、t は時刻のリスト、args はパラメータである。sol に結果が入る。

結果を図示すると次のようになる。

plt.plot(t, sol[:, 0])
plt.xlabel("t [s]")
plt.ylabel("x [m]")
plt.show()

毛管現象

毛管現象の方程式を解く。パラメータは以下の通り。

rho = 1000.
mu = 1e-3
d = 1e-5
r = d/2.
sigma = 70e-3
theta = 0.

rho は密度、mu は粘性係数、d は毛管径、r は毛管半径、sigma は表面張力、theta は接触角である。

重力とバランスしたときの高さ (Jurin の高さ) は、次のように計算できる。

h = 2.*sigma*np.cos(theta)/(rho*r*g)
print("h = ", h*1000., "[mm]")

結果。

h =  2854.2303771661573 [mm]

関数を用意。

def func(s, t):
    z, v = s
    m = rho*np.pi*r**2*z
    f = 2.*np.pi*r*sigma*np.cos(theta) - 8.*np.pi*mu*z*v - m*g
    dsdt = [v, f/m]
    return dsdt

z が液柱の高さ、v が速度である。運動方程式を解くのは自由落下の場合と同じだが、この場合は質量が液柱の高さに依存する。かかる力としては、吸い上げる表面張力と重力の他に、粘性抵抗があり、これはポアズイユの法則による。

実行。

z0 = [5e-3, 0]
t = np.linspace(0, 500000, 500001)
sol = odeint(func, z0, t)

描画。

plt.plot(t, sol[:, 0]*1000.)
plt.plot([0, 500000], [h*1000., h*1000.])
plt.xlim([0, 500000])
plt.ylim(0)
plt.xlabel("t [s]")
plt.ylabel("z [mm]")
plt.grid()
plt.show()