Claude Codeでsolve_ivpの微分方程式を解析解で検算する

solve_ivpは、エラーなく終わっても答えが正しいとは限りません。許容誤差の既定値はrtol=1e-3、atol=1e-6で、結果は設定した精度の範囲でしか合っていないからです。Claude Codeに微分方程式のコードを書かせるなら、解析解が分かる系で検算するテストを同時に置くのが近道です。

この記事では、減衰、振動子、硬い系の3つを解析解と突き合わせます。許容誤差とメソッドを変えると誤差がどう動くか、scipy 1.18.1 / numpy 2.5.4で走らせた数値で並べます。

解析解が分かる系をテストオラクルにする

数値解の正否は、正解を知っている系で確かめます。この正解役がテストオラクルです。選ぶ条件は、手で解けて、本番のコードと同じ関数を通せることです。

硬い系の解析解は、微分して元の式に戻せば確かめられます。Claudeが出した解析解をそのまま信じず、sympy.dsolveや手計算で自分でも一度確認します。オラクル自体が誤っていれば、検算は意味を失います。

指数減衰で許容誤差の効き方を見る

最小の系から始めます。区間は0〜5、評価点はt_evalで51点に固定しました。

import numpy as np
from scipy.integrate import solve_ivp
 
def decay(t, y):
    return -2.0 * y
 
t = np.linspace(0, 5, 51)
exact = np.exp(-2.0 * t)
 
sol = solve_ivp(decay, (0, 5), [1.0], t_eval=t)  # 既定の許容誤差
print(np.abs(sol.y[0] - exact).max(), sol.nfev)

既定のままだと、最大の絶対誤差は約4.1e-4、右辺の評価回数(nfev)は74でした。rtol=1e-8, atol=1e-10にすると誤差は約7.4e-10、評価回数は512です。精度を上げた分だけ、右辺の評価は約7倍に増えました。

もう一つ見落としやすい点があります。t=5での真値は約4.5e-5と小さく、既定のatol=1e-6が相対許容のrtol * abs(y)に近づきます。実際、この点の相対誤差は約0.71%でした。公式の説明にある「atolがrtol * abs(y)より大きいと、正しい桁数は保証されない」の例です。

atolは、解が取りうる最小の大きさより十分小さく置きます。成分ごとにスケールが違う連立系では、atolに配列を渡して成分ごとに決められます。

単振動で長時間の誤差を見る

減衰は時間とともに誤差が縮みますが、振動は誤差が溜まります。y'' = -yを連立にして、0〜20πの10周期を400点で比べました。

def osc(t, y):
    return [y[1], -y[0]]
 
T = 20 * np.pi
t = np.linspace(0, T, 400)
sol = solve_ivp(osc, (0, T), [1, 0], method="DOP853",
                t_eval=t, rtol=1e-9, atol=1e-12)
print(np.abs(sol.y[0] - np.cos(t)).max())

既定の設定では、振幅1の波に対して誤差が約1.8%出ています。見た目の波形は滑らかで、グラフだけでは気づきません。許容誤差を絞ると、RK45は評価回数が約19倍に増えましたが、DOP853は約5倍で済みました。公式の説明は、高精度が必要なときのRunge-KuttaにはDOP853を勧めています。この結果はその記述と整合します。

次数の低いRK23も同じ条件で試しました。既定の許容誤差では最大誤差が約5.1e-2、評価回数は974でRK45の約2.8倍の誤差です。rtol=1e-9, atol=1e-12に絞ると、誤差は約5.2e-8まで下がりますが、評価回数は82010に膨らみます。RK45の7904回と比べて10倍を超えるので、精度を詰める検算には向きません。

解析解を持たない系では、保存量が代わりの目安になります。振動子ではy² + y'² が常に1なので、np.abs(sol.y[0]**2 + sol.y[1]**2 - 1).max()で確かめられます。既定の設定のRK45では、この値が約1.3e-3でした。保存量のずれが閾値内かを見るテストなら、解析解がない本番の系にも流用できます。

硬い系でRK45とRadauを比べる

硬い系は、速く消える成分と遅い成分が混在します。-1000yの項が、すぐに消える成分を作ります。

def stiff(t, y):
    return -1000.0 * y + 3000.0 - 2000.0 * np.exp(-t)
 
def exact(t):
    return 3 - 0.998 * np.exp(-1000 * t) - 2.002 * np.exp(-t)
 
t = np.linspace(0, 10, 101)
for m in ["RK45", "Radau", "BDF"]:
    s = solve_ivp(stiff, (0, 10), [0.0], method=m, t_eval=t,
                  rtol=1e-6, atol=1e-9)
    print(m, np.abs(s.y[0] - exact(t)).max(), s.nfev, s.success)

RK45もsuccessはTrueで、答えは合っています。違うのは右辺の評価回数で、RK45はRadauの約43倍でした。誤差だけを見ればRK45も合格ですが、評価回数が桁違いに増えるのが硬い系の兆候です。公式の説明も、RK45が異常に多くの反復をするなら硬い問題を疑い、RadauかBDFを使うよう書いています。

Radauの誤差がRK45より一桁大きい点には注意が要ります。手法が違えば、同じrtolでも実際の誤差は同じになりません。許容誤差は目標値であって、保証ではありません。

ヤコビアンを渡すと評価回数が減る、LSODAは自動で切り替える

Radau、BDF、LSODAは、右辺をyで微分した行列(ヤコビアン)を使います。jacを渡さないと有限差分で近似されます。ここではjac=lambda t, y: [[-1000.0]]を渡して、同じ条件で比べました。

誤差はほとんど変わらず、nfevだけが減りました。BDFで約16%、LSODAで約1割減り、Radauでは3%ほどでした。1変数の系なので差は小さめですが、成分が多いほど有限差分の評価が増えるため、差は開きます。ヤコビアンは「渡したほうがよい」とされており、検算でも誤差が変わらないことを確かめたうえで渡せます。

LSODAは、系の硬さを自動で判定して手法を切り替えます。この系ではBDFと同程度の誤差で、nfevは最少でした。RK45とRadauのどちらにするか決めきれないときの候補です。ただし、同じ振動子(0〜20π、rtol=1e-9)に使うとnfevは1853で、DOP853の2246とほぼ同じ規模でした。硬くない系で特別に得をするわけではありません。

解析解の任意の点と、到達時刻まで検算する

t_evalで決めた点以外でも比べたいときは、dense_output=Trueを使います。sol.sol(t)が、任意の時刻での値を返します。指数減衰を既定の許容誤差で解き、0〜5を1001点で解析解と比べた最大誤差は約4.3e-4でした。ステップ数は13点だけです。許容誤差をrtol=1e-8, atol=1e-10に絞ると、誤差は約2.4e-9まで下がりました。補間の精度は、ソルバーの許容誤差にそのまま従います。

もう一つ、eventsは「値がある閾値を切る時刻」の検算に使えます。減衰が0.01になる時刻は、解析解ならln(100)/2 ≒ 2.302585です。

ev = lambda t, y: y[0] - 0.01
ev.terminal = True      # 到達したら積分を止める
ev.direction = -1       # 下向きに横切るときだけ検出する
sol = solve_ivp(decay, (0, 5), [1.0], events=ev)
print(sol.t_events[0][0], np.log(100) / 2)

既定の許容誤差では、検出された時刻は約2.30391で、解析解との差は約1.3e-3でした。rtol=1e-8, atol=1e-10なら差は約5.5e-9です。値の誤差が4e-4程度でも、時刻に直すとその3倍ほどのずれになります。「いつ超えるか」を出す計算では、値の検算より厳しい許容誤差が要ると分かります。directionを指定しないと、上向きの横切りも拾われます。振動する系では、同じ閾値を何度も横切るため、検出した時刻が検算したい回かどうかも確かめます。terminalを付けなければ、t_eventsに全部の時刻が並びます。

pytestに落として回帰テストにする

検算を一度の実験で終わらせず、テストに残します。パラメータ化で3メソッドを回し、successの確認も入れます。

import numpy as np
import pytest
from scipy.integrate import solve_ivp
 
@pytest.mark.parametrize("method", ["RK45", "Radau", "BDF"])
def test_stiff_matches_closed_form(method):
    f = lambda t, y: -1000.0 * y + 3000.0 - 2000.0 * np.exp(-t)
    exact = lambda t: 3 - 0.998 * np.exp(-1000 * t) - 2.002 * np.exp(-t)
    t = np.linspace(0, 10, 101)
    sol = solve_ivp(f, (0, 10), [0.0], method=method, t_eval=t,
                    rtol=1e-6, atol=1e-9)
    assert sol.success
    np.testing.assert_allclose(sol.y[0], exact(t), atol=1e-4)

減衰と振動子のテストを加えた3関数で、scipy 1.18.1では5件が通りました。閾値は実測誤差の数倍に置きます。厳しすぎると環境差で落ち、緩すぎるとバグを見逃します。assert_allcloseのrtolとatolは、solve_ivpの許容誤差とは別物なので、混同しないでください。

parametrizeの組み方は、pytestのparametrizeとfixtureの記事に詳しくあります。

Claude Codeに任せる範囲と指示の置き方

公式のワークフロー集には、テスト追加では「検証したい挙動を具体的に伝える」と書かれています。数値計算では、これを次のような依頼に変えられます。

dy/dt = -2y を solve_ivp で解く関数を書いて。
解析解 exp(-2t) との最大絶対誤差を計算するテストも pytest で追加。
実行して、誤差と nfev を表で報告して。

実行結果はClaudeに読ませますが、閾値の根拠は人が決めます。CLAUDE.mdに規約として置くなら、例えば次のような形です。

## 数値計算のルール
- solve_ivp には rtol と atol を明示する。既定値のまま使わない
- 解析解が分かる系を必ずテストに入れる(減衰・単振動など)
- 結果の success を assert する。method を変えたら誤差と nfev を残す

Claudeは既存のテストの書き方に合わせてコードを書く、とワークフロー集にあります。最初の1本を人が書いて型を示しておくと、以降のテストがその型に揃います。

Claudeが書いた検算テストで見るべき3か所

テストが緑でも、検算として働いているとは限りません。Claudeの出力で、次の3か所を確かめます。

1. 解析解の係数と符号: dy/dt = -2yの解をexp(-t)と書くと、最大誤差は約0.25で、閾値が1e-2でも落ちます。怖いのは小さなずれです。exp(-2.01t)と書いた場合の誤差は約1.8e-3で、閾値を1e-2と緩く置いたテストは通ってしまいます。実測した数値誤差(約4.1e-4)の数倍に閾値を置けば、このずれは検出できます
2. 比較する点のそろえ方: t_evalを渡さずに解くと、この系では解が13点で返ります。51点の解析解とは形が合わず、そのまま引き算するとエラーになります。無理に補間して合わせると、原因の切り分けが難しくなります。t_evalかdense_outputのどちらかに統一させます
3. assert_allcloseの許容値: rtol=1e-3だけを指定してatolを省くと、atolは0です。減衰の末尾のように真値が小さい点では、相対誤差が大きいため、既定の許容誤差で解いた結果はこのテストを通りません。solve_ivpのatolを理由にassert_allcloseのatolを緩める前に、どの点で落ちているかを見ます

依頼文にも「実行して誤差とnfevを報告して」と入れておくと、閾値を人が決める材料が揃います。

検算でつまずきやすい点

• t_evalを省いた比較: t_evalがないと、解はソルバーが選んだ点でしか得られません。解析解と比べる点が合わず、補間を挟んで誤差の原因が混ざります。比較点が決まっているならt_evalを渡します
• successを見ない: ステップ失敗でも例外は出ません。statusが-1なら途中で止まっています。assert sol.successを最初に置きます
• 連立系の取り違え: 2階の式は連立の1階に直します。sol.yは行が成分、列が時刻です。sol.y[0]が位置、sol.y[1]が速度だと確かめます
• 硬い系にRK45を使い続ける: 正しい答えは出ますが、評価回数が膨らみます。nfevも一緒に記録しておくと気づけます

同じ「真値を先に置く」手法を、道路線形の計算に使った例がクロソイドの検算の記事にあります。ピーク検出の面積を検算する流れはクロマトグラムの記事が扱っています。

まとめ

solve_ivpの検算で効くのは、既定の許容誤差を信じないことと、評価回数を一緒に残すことの2点です。誤差が小さくてもnfevが桁違いなら、メソッドを替える合図になります。解析解のある3系をpytestに置いておけば、scipyの更新やコードの書き換えでも数値の変化に気づけます。