100本ノック / 数値計算 / 数値計算100本ノック

製造業の常微分方程式入門｜加熱・反応・設備振動をPythonでシミュレーション

2020-06-29

加熱・反応・搬送の時間変化を予測する：製造業の意思決定に効く常微分方程式10本ノック（No.071〜No.080）

本記事では、常微分方程式（ODE）の代表的な数値解法を、架空の精密部品工場における熱処理条件の設計、反応槽の品質予測、設備振動、感染・欠勤リスク、加熱炉の運転計画へ結び付けます。目的は計算公式を覚えることではなく、計算精度、計算時間、モデルの妥当性を区別し、操業条件の意思決定に使えるシミュレーションを設計することです。

扱う範囲は No.071〜No.080（オイラー法、改良オイラー法、Runge–Kutta法、適応刻み幅法、剛性方程式、Adams法、BDF法、振り子、SIRモデル、物理シミュレーション）です。掲載データとパラメータはすべて説明用の架空値です。

[!NOTE] 本資料は、数理工房 (もしくは代表である和山個人) が過去に企業研修において使用した notebook を企業様の許可を得て再構成・編集のうえ公開しています。
掲載データはすべて架空のものであり、実在する企業・工場・数値とは一切関係ありません。

はじめに：この記事で扱う製造業の実務課題

製造現場では温度、濃度、圧力、振動、在庫などが時間とともに変化します。「現在値」だけでなく、投入量や外気温を変えたときにいつ規格内へ入り、どこまで行き過ぎ、異常時に何分の余裕があるかを知る必要があります。本記事では、状態 $y(t)$ の変化率を

$\frac{dy}{dt}=f(t,y;\theta),\qquad y(t_0)=y_0$

と表し、数値積分によって将来軌道を求めます。 $\theta$ は熱容量や反応速度などのモデルパラメータです。

現場でよくある状況

熱処理槽が目標温度へ入る時刻を見積もり、生産開始時刻を決めたい
速い反応と遅い反応が混在し、通常の計算では刻み幅を極端に小さくする必要がある
既存設備の時系列データからモデルを校正し、条件変更の影響を事前評価したい
計算結果が滑らかでも、数値誤差なのかモデル誤差なのか判断できない

なぜこの問題は判断が難しいのか

離散化した計算値には、刻み幅による打切り誤差と丸め誤差があります。一方、実設備との差には、パラメータ推定誤差、未計測外乱、モデル構造の不足も含まれます。高次の解法に変えるだけで実設備予測が正しくなるとは限りません。刻み幅を半分にした収束確認と、独立した実測データによる妥当性確認の両方が必要です。

今回扱うノックの全体像

前半では同じ冷却モデルを複数解法で解き、精度と計算量を比較します。中盤では剛性を持つ反応系と多段法を扱います。後半では振り子、SIRモデル、二領域加熱炉を題材に、モデル化から意思決定までを確認します。

Python 環境の準備

NumPy で配列計算、pandas で表、SciPy で検証済みODEソルバー、matplotlib で可視化を行います。乱数シードを固定し、同じ結果を再現できるようにします。

%matplotlib inline
import platform
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

SEED = 20260712
rng = np.random.default_rng(SEED)
pd.set_option("display.precision", 4)
print(f"Python: {platform.python_version()}")
print(f"NumPy: {np.__version__}, pandas: {pd.__version__}")
print(f"matplotlib: {matplotlib.__version__}")

Python: 3.13.1
NumPy: 2.5.1, pandas: 3.0.3
matplotlib: 3.11.0

架空データの作成

最初の比較には、炉から取り出した部品の冷却を表すNewtonの冷却則を使います。部品温度 $T$ 、周囲温度 $T_a$ 、冷却係数 $k$ に対して

$\frac{dT}{dt}=-k(T-T_a),\qquad T(t)=T_a+(T_0-T_a)e^{-kt}$

です。後半では反応濃度、振り子、欠勤リスク、炉内二領域の架空パラメータを別途定義します。解析解は数値解法の誤差を測る基準としてのみ使います。

T0, T_ambient, k_cool = 180.0, 25.0, 0.18

def cooling_rhs(t, y):
    return np.asarray([-k_cool * (y[0] - T_ambient)])

def cooling_exact(t):
    return T_ambient + (T0 - T_ambient) * np.exp(-k_cool * np.asarray(t))

base_times = np.arange(0, 21, 2)
fictional_measurements = cooling_exact(base_times) + rng.normal(0, 1.2, len(base_times))
cooling_data = pd.DataFrame({"time_min": base_times, "measured_temp_C": fictional_measurements})
display(cooling_data.head())

plt.figure(figsize=(8, 4))
plt.scatter(base_times, fictional_measurements, label="Fictional measurements")
tt = np.linspace(0, 20, 200)
plt.plot(tt, cooling_exact(tt), label="Model curve")
plt.title("Fictional workpiece cooling data")
plt.xlabel("Time [min]"); plt.ylabel("Temperature [°C]")
plt.grid(True, alpha=0.3); plt.legend(); plt.tight_layout(); plt.show()

	time_min	measured_temp_C
0	0	180.8863
1	2	133.9578
2	4	100.9478
3	6	78.1300
4	8	62.2210

png

No.071：オイラー法

実務での意味

オイラー法は、現在の変化率が短時間続くと仮定して次の状態を予測します。PLCの簡易予測やモデルの概念検証に使いやすく、刻み幅が判断結果へ与える影響を理解する出発点です。

分析・モデル化の考え方

刻み幅を $h$ とすると、 $y_{n+1}=y_n+h f(t_n,y_n)$ です。1ステップの局所誤差は $O(h^2)$ 、一定区間を進んだ大域誤差は $O(h)$ です。したがって刻み幅を半分にすると、十分小さい範囲では誤差がおおむね半分になります。

Pythonで確認する

def euler(f, t0, tf, y0, h):
    t = np.arange(t0, tf + h/2, h)
    y = np.empty((len(t), len(np.atleast_1d(y0))))
    y[0] = np.atleast_1d(y0)
    for n in range(len(t) - 1):
        y[n+1] = y[n] + h * f(t[n], y[n])
    return t, y

rows = []
plt.figure(figsize=(8, 4))
for h in [2.0, 1.0, 0.5]:
    t_e, y_e = euler(cooling_rhs, 0, 20, [T0], h)
    err = abs(y_e[-1, 0] - cooling_exact(20))
    rows.append({"step_min": h, "T_at_20_C": y_e[-1, 0], "abs_error_C": err})
    plt.plot(t_e, y_e[:, 0], marker="o", ms=3, label=f"Euler h={h}")
plt.plot(tt, cooling_exact(tt), "k--", label="Exact")
display(pd.DataFrame(rows))
plt.title("Euler method: effect of step size")
plt.xlabel("Time [min]"); plt.ylabel("Temperature [°C]")
plt.grid(True, alpha=0.3); plt.legend(); plt.tight_layout(); plt.show()

	step_min	T_at_20_C	abs_error_C
0	2.0	26.7870	2.4481
1	1.0	27.9283	1.3069
2	0.5	28.5644	0.6708

png

結果の読み取り

刻み幅を小さくすると20分後の誤差が減ります。粗い刻みでも曲線は滑らかに見えるため、見た目だけで採用してはいけません。たとえば「60°C以下で搬送可」という判定では、温度誤差を時間余裕へ換算し、安全側の許容誤差を先に決めます。

No.072：改良オイラー法

実務での意味

加熱・冷却速度が区間内で変わる場合、区間始点の傾きだけでは偏りが生じます。改良オイラー法（Heun法）は始点と予測終点の傾きを平均し、少ない追加計算で精度を上げます。

分析・モデル化の考え方

予測値 $\tilde y_{n+1}=y_n+h f(t_n,y_n)$ を作り、

$y_{n+1}=y_n+\frac{h}{2}\{f(t_n,y_n)+f(t_{n+1},\tilde y_{n+1})\}$

と補正します。大域誤差は $O(h^2)$ です。同じ刻み幅でオイラー法と比較し、手法変更の効果を定量化します。

Pythonで確認する

def heun(f, t0, tf, y0, h):
    t = np.arange(t0, tf + h/2, h)
    y = np.empty((len(t), len(np.atleast_1d(y0)))); y[0] = y0
    for n in range(len(t)-1):
        k1 = f(t[n], y[n])
        predictor = y[n] + h * k1
        k2 = f(t[n+1], predictor)
        y[n+1] = y[n] + h * (k1 + k2) / 2
    return t, y

h = 2.0
t_e, y_e = euler(cooling_rhs, 0, 20, [T0], h)
t_h, y_h = heun(cooling_rhs, 0, 20, [T0], h)
compare = pd.DataFrame({
    "method": ["Euler", "Improved Euler"],
    "T_at_20_C": [y_e[-1, 0], y_h[-1, 0]],
    "abs_error_C": [abs(y_e[-1, 0]-cooling_exact(20)), abs(y_h[-1, 0]-cooling_exact(20))]
})
display(compare)
plt.figure(figsize=(8, 4))
plt.plot(t_e, y_e[:, 0], "o-", label="Euler")
plt.plot(t_h, y_h[:, 0], "s-", label="Improved Euler")
plt.plot(tt, cooling_exact(tt), "k--", label="Exact")
plt.title("Euler and improved Euler at the same step size")
plt.xlabel("Time [min]"); plt.ylabel("Temperature [°C]")
plt.grid(True, alpha=0.3); plt.legend(); plt.tight_layout(); plt.show()

	method	T_at_20_C	abs_error_C
0	Euler	26.787	2.4481
1	Improved Euler	29.688	0.4529

png

結果の読み取り

同じ2分刻みでも改良オイラー法の誤差は大幅に小さくなります。ただし関数評価回数は1ステップ当たり2回です。組込み制御では、許容誤差、演算時間、制御周期の締切を並べて方式を決めます。

No.073：Runge–Kutta法

実務での意味

4次Runge–Kutta法（RK4）は、熱・反応・運動モデルで広く使われる固定刻み解法です。1ステップ内の複数地点で傾きを評価し、運転条件比較に必要な精度を得やすくします。

分析・モデル化の考え方

$k_1=f(t_n,y_n)$ 、 $k_2=f(t_n+h/2,y_n+hk_1/2)$ 、 $k_3=f(t_n+h/2,y_n+hk_2/2)$ 、 $k_4=f(t_n+h,y_n+hk_3)$ として、

$y_{n+1}=y_n+\frac{h}{6}(k_1+2k_2+2k_3+k_4)$

とします。大域誤差は $O(h^4)$ ですが、安定性が無条件に保証されるわけではありません。

Pythonで確認する

def rk4(f, t0, tf, y0, h):
    t = np.arange(t0, tf + h/2, h)
    y = np.empty((len(t), len(np.atleast_1d(y0)))); y[0] = y0
    for n in range(len(t)-1):
        k1 = f(t[n], y[n]); k2 = f(t[n]+h/2, y[n]+h*k1/2)
        k3 = f(t[n]+h/2, y[n]+h*k2/2); k4 = f(t[n]+h, y[n]+h*k3)
        y[n+1] = y[n] + h*(k1+2*k2+2*k3+k4)/6
    return t, y

method_rows = []
for name, solver, evaluations in [("Euler", euler, 1), ("Improved Euler", heun, 2), ("RK4", rk4, 4)]:
    t_m, y_m = solver(cooling_rhs, 0, 20, [T0], 2.0)
    method_rows.append({"method": name, "rhs_calls_per_step": evaluations,
                        "abs_error_at_20_C": abs(y_m[-1,0]-cooling_exact(20))})
display(pd.DataFrame(method_rows))

t_r, y_r = rk4(cooling_rhs, 0, 20, [T0], 2.0)
plt.figure(figsize=(8, 4))
plt.plot(t_r, y_r[:, 0], "o-", label="RK4 h=2")
plt.plot(tt, cooling_exact(tt), "k--", label="Exact")
plt.title("Fourth-order Runge–Kutta cooling trajectory")
plt.xlabel("Time [min]"); plt.ylabel("Temperature [°C]")
plt.grid(True, alpha=0.3); plt.legend(); plt.tight_layout(); plt.show()

	method	rhs_calls_per_step	abs_error_at_20_C
0	Euler	1	2.4481
1	Improved Euler	2	0.4529
2	RK4	4	0.0029

png

結果の読み取り

RK4は同じ刻み幅でさらに小さい誤差を示します。単純なモデルでは十分な選択肢ですが、入力が急変する時刻や安全判定の直前には出力点を細かくします。高次法でも、不連続な制御入力をまたぐ1ステップは分割する必要があります。

No.074：適応刻み幅法

実務での意味

炉の立上げ直後は温度変化が速く、定常付近では遅くなります。全区間を細かく刻む代わりに、誤差が大きい区間だけ刻み幅を縮めれば、精度と計算時間を両立できます。

分析・モデル化の考え方

埋込み型Runge–Kutta法は異なる次数の近似差から局所誤差を推定します。rtol は状態の大きさに比例する許容誤差、atol はゼロ近傍の絶対許容誤差です。許容尺度は温度計精度や判定マージンから決めます。

Pythonで確認する

adaptive_rows = []
solutions = {}
for rtol in [1e-3, 1e-6, 1e-9]:
    sol = solve_ivp(cooling_rhs, (0, 20), [T0], method="RK45", rtol=rtol, atol=rtol*0.1)
    solutions[rtol] = sol
    adaptive_rows.append({"rtol": rtol, "accepted_points": len(sol.t), "rhs_calls": sol.nfev,
                          "abs_error_at_20_C": abs(sol.y[0,-1]-cooling_exact(20))})
display(pd.DataFrame(adaptive_rows))
sol = solutions[1e-6]
steps = np.diff(sol.t)
plt.figure(figsize=(8, 4))
plt.step(sol.t[:-1], steps, where="post")
plt.title("Adaptive RK45 accepted step sizes")
plt.xlabel("Time [min]"); plt.ylabel("Accepted step size [min]")
plt.grid(True, alpha=0.3); plt.tight_layout(); plt.show()

	rtol	accepted_points	rhs_calls	abs_error_at_20_C
0	1.0000e-03	7	38	1.0039e-02
1	1.0000e-06	17	98	8.0528e-06
2	1.0000e-09	56	332	8.0192e-09

png

結果の読み取り

許容誤差を厳しくすると関数評価回数が増え、終点誤差が小さくなります。ソルバーが受理した内部時刻と、報告・制御に必要な出力時刻は別物です。品質判定時刻は t_eval やイベント検出で明示し、内部ステップ数だけをKPIにしません。

No.075：剛性方程式

実務での意味

洗浄反応や触媒反応では、瞬時に近い速い成分と、製造時間を支配する遅い成分が共存します。このような剛性系では、陽解法は精度ではなく安定性のために極小刻みを要求されます。

分析・モデル化の考え方

架空の連続反応 $A\to B\to C$ を

$A'=-k_1A,\quad B'=k_1A-k_2B,\quad C'=k_2B$

とします。 $k_1\gg k_2$ のとき時間尺度が大きく異なります。陽的RK法と、剛性向けの陰的Radau法を同じ許容誤差で比較します。

Pythonで確認する

k1, k2 = 1000.0, 0.12
def stiff_reaction(t, y):
    A, B, C = y
    return [-k1*A, k1*A-k2*B, k2*B]

stiff_rows = []
stiff_solutions = {}
for method in ["RK45", "Radau"]:
    sol_s = solve_ivp(stiff_reaction, (0, 20), [1, 0, 0], method=method,
                      rtol=1e-6, atol=1e-9, dense_output=True)
    stiff_solutions[method] = sol_s
    stiff_rows.append({"method": method, "rhs_calls": sol_s.nfev,
                       "accepted_points": len(sol_s.t), "mass_balance_error": abs(sol_s.y[:,-1].sum()-1)})
display(pd.DataFrame(stiff_rows))
t_plot = np.linspace(0, 20, 300)
ys = stiff_solutions["Radau"].sol(t_plot)
plt.figure(figsize=(8, 4))
for i, label in enumerate(["A", "B", "C"]): plt.plot(t_plot, ys[i], label=label)
plt.title("Stiff consecutive reaction solved by Radau")
plt.xlabel("Time [min]"); plt.ylabel("Normalized concentration")
plt.grid(True, alpha=0.3); plt.legend(); plt.tight_layout(); plt.show()

	method	rhs_calls	accepted_points	mass_balance_error
0	RK45	42542	6087	1.1102e-15
1	Radau	949	135	7.6688e-09

png

結果の読み取り

同じ精度条件でも陽的RK45は速い反応に拘束され、多数の関数評価を要します。Radauのような陰解法は1ステップが重い一方、剛性系では総計算量を抑えられます。濃度合計が1に保たれるか、負濃度が出ないかも受入試験に含めます。

No.076：Adams法

実務での意味

モデル評価が高価で、状態が滑らかに変化する連続運転では、過去の傾きを再利用する多段法が有効です。ここでは2段Adams–Bashforth法（AB2）で、毎ステップの新規評価を1回に抑えます。

分析・モデル化の考え方

$y_{n+1}=y_n+h\left(\frac32f_n-\frac12f_{n-1}\right)$

は陽的2段法です。最初の1点は別法で作る必要があり、ここではRK4を1ステップだけ使います。過去情報を使うため、刻み幅変更や不連続入力の扱いは1段法より複雑です。

Pythonで確認する

def adams_bashforth2(f, t0, tf, y0, h):
    t = np.arange(t0, tf+h/2, h)
    y = np.empty((len(t), len(np.atleast_1d(y0)))); y[0] = y0
    k1 = f(t[0], y[0]); k2 = f(t[0]+h/2, y[0]+h*k1/2)
    k3 = f(t[0]+h/2, y[0]+h*k2/2); k4 = f(t[0]+h, y[0]+h*k3)
    y[1] = y[0] + h*(k1+2*k2+2*k3+k4)/6
    f_prev, f_now = f(t[0], y[0]), f(t[1], y[1])
    for n in range(1, len(t)-1):
        y[n+1] = y[n] + h*(1.5*f_now-0.5*f_prev)
        f_prev, f_now = f_now, f(t[n+1], y[n+1])
    return t, y

ab_rows = []
for h in [1.0, 0.5, 0.25]:
    t_ab, y_ab = adams_bashforth2(cooling_rhs, 0, 20, [T0], h)
    ab_rows.append({"step_min": h, "steps": len(t_ab)-1,
                    "abs_error_at_20_C": abs(y_ab[-1,0]-cooling_exact(20))})
display(pd.DataFrame(ab_rows))
plt.figure(figsize=(8, 4))
plt.plot(t_ab, y_ab[:,0], "o-", ms=3, label="AB2 h=0.25")
plt.plot(tt, cooling_exact(tt), "k--", label="Exact")
plt.title("Two-step Adams–Bashforth method")
plt.xlabel("Time [min]"); plt.ylabel("Temperature [°C]")
plt.grid(True, alpha=0.3); plt.legend(); plt.tight_layout(); plt.show()

	step_min	steps	abs_error_at_20_C
0	1.00	20	0.2181
1	0.50	40	0.0530
2	0.25	80	0.0131

png

結果の読み取り

刻み幅を半分にすると、十分小さい領域では誤差がおおむね4分の1へ近づき、2次精度が確認できます。AB2は滑らかな非剛性問題向けです。設備ON/OFFや材料投入時刻では履歴を再初期化し、不連続点をまたいで過去傾きを外挿しません。

No.077：BDF法

実務での意味

BDF（後退差分公式）は、剛性のある反応・熱収支モデルで使われる陰的多段法です。長い定常運転の中に速い過渡が含まれる場合、安定性を確保しながら時間を進めます。

分析・モデル化の考え方

1次BDF（後退オイラー法）は

$y_{n+1}=y_n+h f(t_{n+1},y_{n+1})$

で、未知の $y_{n+1}$ が右辺にも現れます。各ステップで非線形方程式を解く必要があります。実務では自作せず、誤差制御、次数変更、ヤコビアン処理を備えた検証済み実装を使います。

Pythonで確認する

bdf_rows = []
bdf_solutions = {}
for method in ["Radau", "BDF"]:
    sol_b = solve_ivp(stiff_reaction, (0, 20), [1, 0, 0], method=method,
                      rtol=1e-6, atol=1e-9, dense_output=True)
    bdf_solutions[method] = sol_b
    bdf_rows.append({"method": method, "rhs_calls": sol_b.nfev, "jacobian_evals": sol_b.njev,
                     "lu_decompositions": sol_b.nlu, "final_C": sol_b.y[2,-1]})
display(pd.DataFrame(bdf_rows))
t_zoom = np.geomspace(1e-6, 1, 240)
plt.figure(figsize=(8, 4))
for method in ["Radau", "BDF"]:
    plt.semilogx(t_zoom, bdf_solutions[method].sol(t_zoom)[1], label=f"B ({method})")
plt.title("BDF and Radau during the fast transient")
plt.xlabel("Time [min, log scale]"); plt.ylabel("Intermediate B concentration")
plt.grid(True, which="both", alpha=0.3); plt.legend(); plt.tight_layout(); plt.show()

	method	rhs_calls	jacobian_evals	lu_decompositions	final_C
0	Radau	949	3	62	0.9093
1	BDF	385	1	38	0.9093

png

結果の読み取り

BDFとRadauがほぼ同じ濃度軌道を返すことは有力な相互確認です。ただし一致はモデルの正しさを証明しません。実測との比較に加え、反応速度を振った感度分析、物質収支、非負性を確認し、ヤコビアンを提供した場合の計算短縮も評価します。

No.078：振り子のシミュレーション

実務での意味

搬送アームや吊り荷の揺れは、停止精度、サイクルタイム、安全余裕へ影響します。振り子は2階の運動方程式を1階の連立ODEへ直す基本例であり、設備の残留振動評価へつながります。

分析・モデル化の考え方

角度 $\theta$ 、角速度 $\omega$ 、長さ $L$ 、減衰 $c$ に対し、

$\theta'=\omega,\qquad \omega'=-\frac{g}{L}\sin\theta-c\omega$

とします。小角近似 $\sin\theta\approx\theta$ は大振幅で誤差が増えるため、非線形モデルとの差を確認します。

Pythonで確認する

g, L, damping = 9.81, 1.2, 0.08
def pendulum(t, y): return [y[1], -(g/L)*np.sin(y[0])-damping*y[1]]
def pendulum_linear(t, y): return [y[1], -(g/L)*y[0]-damping*y[1]]

t_eval = np.linspace(0, 10, 501)
y0_p = [np.deg2rad(55), 0]
sol_nl = solve_ivp(pendulum, (0,10), y0_p, t_eval=t_eval, rtol=1e-9, atol=1e-11)
sol_li = solve_ivp(pendulum_linear, (0,10), y0_p, t_eval=t_eval, rtol=1e-9, atol=1e-11)
angle_diff = np.max(np.abs(np.rad2deg(sol_nl.y[0]-sol_li.y[0])))
print(f"Maximum angle difference from small-angle model: {angle_diff:.2f} deg")
plt.figure(figsize=(8, 4))
plt.plot(t_eval, np.rad2deg(sol_nl.y[0]), label="Nonlinear")
plt.plot(t_eval, np.rad2deg(sol_li.y[0]), "--", label="Small-angle approximation")
plt.title("Damped suspended-load swing")
plt.xlabel("Time [s]"); plt.ylabel("Angle [deg]")
plt.grid(True, alpha=0.3); plt.legend(); plt.tight_layout(); plt.show()

Maximum angle difference from small-angle model: 38.59 deg

png

結果の読み取り

初期角55度では小角近似との位相差が蓄積します。停止待ち時間や干渉判定には非線形モデルを使うべきです。実設備へ導入する際は吊り長さ、摩擦、走行加速度を同定し、荷重条件別に最大角度と整定時間を検証します。

No.079：SIRモデル

実務での意味

感染症流行による欠勤は、人員配置と納期リスクへ影響します。SIRモデルは診断ツールではなく、接触抑制や応援要員準備のシナリオ比較に使う簡略モデルです。

分析・モデル化の考え方

就業者を感受性あり $S$ 、感染中 $I$ 、回復 $R$ に分け、

$S'=-\beta SI/N,\quad I'=\beta SI/N-\gamma I,\quad R'=\gamma I$

とします。基本再生産数の単純指標は $R_0=\beta/\gamma$ です。ここでは接触率を25%下げるシナリオと比較します。

Pythonで確認する

N, gamma = 240, 1/7
def sir_rhs(beta):
    def rhs(t, y):
        S, I, R = y
        return [-beta*S*I/N, beta*S*I/N-gamma*I, gamma*I]
    return rhs

sir_rows = []
plt.figure(figsize=(8, 4))
days = np.linspace(0, 80, 321)
for label, beta in [("Baseline", 0.34), ("Contact rate -25%", 0.34*0.75)]:
    sol_sir = solve_ivp(sir_rhs(beta), (0,80), [238,2,0], t_eval=days, rtol=1e-8, atol=1e-10)
    peak_idx = np.argmax(sol_sir.y[1])
    sir_rows.append({"scenario": label, "R0": beta/gamma,
                     "peak_absent_people": sol_sir.y[1,peak_idx], "peak_day": days[peak_idx]})
    plt.plot(days, sol_sir.y[1], label=label)
display(pd.DataFrame(sir_rows))
plt.title("Illustrative SIR absence-risk scenarios")
plt.xlabel("Day"); plt.ylabel("People in infectious compartment")
plt.grid(True, alpha=0.3); plt.legend(); plt.tight_layout(); plt.show()

	scenario	R0	peak_absent_people	peak_day
0	Baseline	2.380	52.5647	25.5
1	Contact rate -25%	1.785	28.7663	37.5

png

結果の読み取り

接触率を下げる仮定ではピーク人数が減り、ピーク時期も変わります。これは施策効果の断定ではなく、応援要員や仕掛在庫の準備量を比較する材料です。実運用では個人情報を集計化し、公的ガイダンスを優先し、部署間接触、潜伏、検査、ワクチンなど省略した要因を明記します。

No.080：物理シミュレーション

実務での意味

炉内を一様温度とみなすと、表面センサーが目標へ達しても部品中心が未加熱という品質リスクがあります。二領域の集中定数モデルで、ヒーター出力候補ごとの昇温時間と温度差を比較します。

分析・モデル化の考え方

炉内空気温度 $T_a$ と部品中心温度 $T_c$ のエネルギー収支を

$C_aT_a'=P-h_a(T_a-T_{env})-h_c(T_a-T_c),\qquad C_cT_c'=h_c(T_a-T_c)$

とします。入力 $P$ を変更し、中心温度が品質下限へ達する時刻と最大温度差をKPIにします。

Pythonで確認する

C_air, C_core, h_loss, h_core, T_env = 18.0, 65.0, 0.45, 2.8, 25.0
def furnace_rhs(power):
    def rhs(t, y):
        T_air, T_core = y
        return [(power-h_loss*(T_air-T_env)-h_core*(T_air-T_core))/C_air,
                h_core*(T_air-T_core)/C_core]
    return rhs

minutes = np.linspace(0, 180, 721)
furnace_rows = []
plt.figure(figsize=(8, 4))
for power in [75.0, 90.0, 105.0]:
    sf = solve_ivp(furnace_rhs(power), (0,180), [25,25], t_eval=minutes, method="RK45",
                   rtol=1e-7, atol=1e-9)
    reached = np.flatnonzero(sf.y[1] >= 150)
    reach_time = minutes[reached[0]] if len(reached) else np.nan
    furnace_rows.append({"power_kJ_per_min": power, "core_150C_time_min": reach_time,
                         "max_air_core_gap_C": np.max(sf.y[0]-sf.y[1]), "core_at_180_C": sf.y[1,-1]})
    plt.plot(minutes, sf.y[1], label=f"Core, P={power:.0f}")
display(pd.DataFrame(furnace_rows))
plt.axhline(150, color="k", ls="--", label="Quality threshold")
plt.title("Two-zone furnace heating scenarios")
plt.xlabel("Time [min]"); plt.ylabel("Core temperature [°C]")
plt.grid(True, alpha=0.3); plt.legend(); plt.tight_layout(); plt.show()

	power_kJ_per_min	core_150C_time_min	max_air_core_gap_C	core_at_180_C
0	75.0	NaN	17.4573	121.3896
1	90.0	NaN	20.9487	140.6676
2	105.0	160.5	24.4402	159.9455

png

結果の読み取り

出力を上げると中心温度150°Cへの到達は早まりますが、空気と中心の最大温度差も確認すべきです。能力増強の判断は到達時間だけでなく、表面過熱、エネルギー原単位、温度均一性、設備上限を含む多目的評価にします。この二領域モデルは候補絞込みに使い、採用条件は実機試験で確認します。

対象ノックを通して見える実務上の示唆

解法は目的と時間尺度で選ぶ：概念確認にはオイラー法、滑らかな非剛性問題にはRK法、速い・遅い現象が混在する系にはRadauやBDFが候補です。
誤差許容値を業務単位で決める：rtol=1e-6 を慣例で選ぶのではなく、「温度±0.5°C」「到達時刻±1分」など判断の許容幅へ結び付けます。
数値誤差とモデル誤差を分ける：刻み幅収束で数値誤差を点検し、別ロット・別条件の実測でモデル誤差を評価します。
保存則と物理制約をテストする：物質収支、エネルギー収支、濃度非負、設備上下限は、グラフの見た目より強い検証材料です。
シミュレーションは比較の道具：一点予測を断定するのではなく、条件・パラメータを振り、意思決定が変わる境界を示します。

実務導入する場合に必要なこと

目的・KPIの合意：品質判定、サイクルタイム、エネルギー、安全余裕のどれを最適化するかを定義する
データ設計：時刻同期、センサー校正、入力・外乱・設備状態・ロット情報を保存する
パラメータ同定と検証分割：校正に使った運転データと、性能評価用データを分ける
不確実性・感度分析：熱容量、反応速度、初期状態を振り、結論が反転する条件を確認する
ソルバー受入試験：許容誤差、保存則、異常入力、イベント時刻、再現性を自動テストする
運用設計：モデル版、パラメータ、入力、実行時刻、結果、承認者を記録し、人が介入できる仕組みにする
段階導入：過去データ再現、シャドー運用、限定工程での試行を経て制御・計画へ接続する

まとめ

常微分方程式の数値解法は、時間変化を伴う製造現象を「条件を変えたらどうなるか」という比較可能な形へ変換します。オイラー法からRK法、適応刻み、多段法、剛性向け解法までの違いは、精度だけでなく安定性と計算コストにあります。実務価値を生むのは高度なソルバー単体ではなく、妥当な収支モデル、現場単位の許容誤差、実測検証、意思決定プロセスまでを一体で設計したときです。

法人向けのご相談

数理工房では、熱・反応・設備運動などのシミュレーション構築、実測データを用いたパラメータ同定、計算精度の検証、PoCから運用までの設計をご支援します。「既存のExcel計算を検証可能なモデルへ移したい」「実験回数を減らしながら条件探索したい」といった段階からご相談いただけます。

📩 お問い合わせ: surikobo.co.jp/contact
まずはお気軽にご相談ください。

製造業の常微分方程式入門｜加熱・反応・設備振動をPythonでシミュレーション

加熱・反応・搬送の時間変化を予測する：製造業の意思決定に効く常微分方程式10本ノック（No.071〜No.080）

はじめに：この記事で扱う製造業の実務課題

現場でよくある状況

なぜこの問題は判断が難しいのか

今回扱うノックの全体像

Python 環境の準備

架空データの作成

No.071：オイラー法

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.072：改良オイラー法

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.073：Runge–Kutta法

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.074：適応刻み幅法

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.075：剛性方程式

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.076：Adams法

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.077：BDF法

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.078：振り子のシミュレーション

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.079：SIRモデル

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

No.080：物理シミュレーション

実務での意味

分析・モデル化の考え方

Pythonで確認する

結果の読み取り

対象ノックを通して見える実務上の示唆

実務導入する場合に必要なこと

まとめ

法人向けのご相談

Related Posts

製造業の工程条件最適化をPythonで実践｜勾配降下法・L-BFGS・制約付き最適化

製造業の偏微分方程式入門｜熱伝導・炉内流れをPythonで可視化

製造業で学ぶ数値積分・微分｜連続炉の温度監視と電力量計算10本ノック

製造業の補間・近似入門｜疎な工程測定から品質曲線をPythonで推定

製造業の固有値問題入門｜設備振動・PCA・工程ネットワークをPythonで分析