フェーズ4数値解法 — 離散化、DAE の発生とその解き方
フェーズ2・3では solve_ivp を「ブラックボックス」として使い、SPM と SPMe を解いた。
本章ではその中身を開ける。偏微分方程式(PDE)を計算機で解ける形に落とす離散化、
時間を進める時間積分法とその安定性、そして DFN モデルで避けて通れない
微分代数方程式(DAE)の発生とその解き方を、テイラー展開だけを頼りに一歩ずつ組み立てる。
フェーズ5で DFN ソルバーを自作するための、数値解析の道具箱を完成させる章である。
4.1 離散化とは — テイラー展開から差分法へ
フェーズ1で導いた DFN の支配方程式は、固相濃度 $c_s(x,r,t)$ [mol/m³] のように 連続な場(空間の各点で値を持つ関数)についての偏微分方程式だった。 連続な関数は無限個の自由度を持つが、計算機が保持できる数値は有限個である。 そこで、解を求めたい領域に有限個の点(または小さな体積)を配置し、 その点での値だけを未知数として方程式を書き直す。 この操作を離散化(discretization)と呼ぶ。
$x_i = i\,\Delta x$:$i$ 番目の格子点の位置 [m]($i = 0, 1, 2, \dots$)。
$c_i$:格子点 $x_i$ における未知関数の値 $c(x_i)$ の近似 [mol/m³]。
$\Delta t$:時間刻み [s] — 時間方向の離散化の間隔(4.4 節で使う)。
4.1.1 テイラー展開から差分近似を作る
離散化の最も素朴な方法が差分法(finite difference method, FDM)である。 微分 $\partial c/\partial x$ を、格子点の値 $c_i$ の組み合わせで近似する。 その出発点は、既知としているテイラー展開だけである。 滑らかな関数 $c(x)$ を $x_i$ のまわりで展開する:
導出:前進・後退・中心差分と打ち切り誤差
$x_i$ の右隣の点 $x_{i+1} = x_i + \Delta x$ での値をテイラー展開で書く。 $c'$ は $\partial c/\partial x$、$c''$ は $\partial^2 c/\partial x^2$ の略記である:
$$ c_{i+1} = c_i + \Delta x\, c'_i + \frac{\Delta x^2}{2} c''_i + \frac{\Delta x^3}{6} c'''_i + \mathcal{O}(\Delta x^4) $$ここで $\mathcal{O}(\Delta x^4)$ は「$\Delta x \to 0$ のとき $\Delta x^4$ に比例して小さくなる項」 を表す記号(ランダウ記号)である。
(1) を $c'_i$ について解く。両辺から $c_i$ を引き、$\Delta x$ で割ると
$$ \frac{c_{i+1} - c_i}{\Delta x} = c'_i + \underbrace{\frac{\Delta x}{2} c''_i + \mathcal{O}(\Delta x^2)}_{\text{打ち切り誤差}} $$左辺を前進差分(forward difference)と呼ぶ。 厳密な微分 $c'_i$ との差(右辺第 2 項以降)を 打ち切り誤差(truncation error)と呼ぶ。 その主要項は $\Delta x$ の 1 乗に比例するので、前進差分は 1 次精度 $\mathcal{O}(\Delta x)$ である。
左隣の点 $x_{i-1} = x_i - \Delta x$ でも同様に展開する ($\Delta x \to -\Delta x$ と置き換えるだけ)。奇数次の項の符号が反転する:
$$ c_{i-1} = c_i - \Delta x\, c'_i + \frac{\Delta x^2}{2} c''_i - \frac{\Delta x^3}{6} c'''_i + \mathcal{O}(\Delta x^4) $$これを $c'_i$ について解けば後退差分(backward difference)が得られる。 これも 1 次精度である:
$$ \frac{c_i - c_{i-1}}{\Delta x} = c'_i - \frac{\Delta x}{2} c''_i + \mathcal{O}(\Delta x^2) $$(1) から (3) を引き算すると、偶数次の項($c_i$、$c''_i$、…)がすべて相殺する:
$$ c_{i+1} - c_{i-1} = 2\Delta x\, c'_i + \frac{\Delta x^3}{3} c'''_i + \mathcal{O}(\Delta x^5) $$$2\Delta x$ で割れば中心差分(central difference)となる。 打ち切り誤差の主要項が $\Delta x^2$ に比例するので、2 次精度である:
$$ \frac{c_{i+1} - c_{i-1}}{2\Delta x} = c'_i + \frac{\Delta x^2}{6} c'''_i + \mathcal{O}(\Delta x^4) $$今度は (1) と (3) を足し算すると、奇数次の項($c'_i$、$c'''_i$)が相殺する:
$$ c_{i+1} + c_{i-1} = 2 c_i + \Delta x^2 c''_i + \frac{\Delta x^4}{12} c''''_i + \mathcal{O}(\Delta x^6) $$整理すれば、拡散方程式に必要な2 階微分の中心差分が得られる。これも 2 次精度:
$$ \frac{c_{i+1} - 2c_i + c_{i-1}}{\Delta x^2} = c''_i + \frac{\Delta x^2}{12} c''''_i + \mathcal{O}(\Delta x^4) $$結果を一覧にまとめる:
| 名称 | 近似式 | 打ち切り誤差の主要項 | 精度 |
|---|---|---|---|
| 前進差分 | $(c_{i+1}-c_i)/\Delta x$ | $+\tfrac{\Delta x}{2}c''$ | $\mathcal{O}(\Delta x)$(1 次) |
| 後退差分 | $(c_i-c_{i-1})/\Delta x$ | $-\tfrac{\Delta x}{2}c''$ | $\mathcal{O}(\Delta x)$(1 次) |
| 中心差分(1 階) | $(c_{i+1}-c_{i-1})/(2\Delta x)$ | $+\tfrac{\Delta x^2}{6}c'''$ | $\mathcal{O}(\Delta x^2)$(2 次) |
| 中心差分(2 階) | $(c_{i+1}-2c_i+c_{i-1})/\Delta x^2$ | $+\tfrac{\Delta x^2}{12}c''''$ | $\mathcal{O}(\Delta x^2)$(2 次) |
物理的な意味:精度の次数は「グリッドを細かくしたときの見返り」
精度が $p$ 次($\mathcal{O}(\Delta x^p)$)であるとは、格子間隔を半分にすると誤差が $1/2^p$ になるということである。1 次精度なら誤差は半分に、2 次精度なら 4 分の 1 になる。 同じ誤差を達成するのに必要な格子点数が大きく違ってくるため、 精度の次数は計算コストに直結する。4.7 節では、実際に計算した解の誤差から この次数を「測定」して、実装が正しいことの証拠として使う。
4.2 有限体積法(FVM)と保存性 — 球座標の完全な導出
差分法は「微分を差分で置き換える」という分かりやすい方法だが、本教材の実装では 有限体積法(finite volume method, FVM)を採用する。 理由は一つ:DFN モデルの方程式はすべて保存則 (リチウムの物質保存、電荷の保存)であり、FVM は保存則を 離散化した後も厳密に保つように作られた方法だからである。
4.2.1 FVM の考え方:体積で積分してから離散化する
一般の保存則は「濃度 $c$ の時間変化 = フラックスの出入り + 生成」という形をしている。 フラックス(流束)$\boldsymbol{q}$ [mol/(m²·s)] を使うと
$$ \frac{\partial c}{\partial t} = -\nabla \cdot \boldsymbol{q} + s, $$ここで $s$ は単位体積あたりの生成率 [mol/(m³·s)] である。 FDM はこの式を「点」で近似したが、FVM は領域を小さな体積 コントロールボリューム(control volume, CV)に分割し、 各 CV について方程式を体積積分する。 $i$ 番目の CV(体積 $V_i$ [m³])で積分し、発散定理 $\int_{V} \nabla\cdot\boldsymbol{q}\, \mathrm{d}V = \oint_{\partial V} \boldsymbol{q}\cdot\boldsymbol{n}\, \mathrm{d}A$ ($\boldsymbol{n}$ は外向き単位法線)を使うと
$$ \frac{\mathrm{d}}{\mathrm{d}t}\int_{V_i} c \,\mathrm{d}V = -\oint_{\partial V_i} \boldsymbol{q}\cdot\boldsymbol{n}\,\mathrm{d}A + \int_{V_i} s \,\mathrm{d}V . $$ここまでは一切近似していない。この式は 「CV 内の総量の変化 = 境界面からの正味の流入 + 内部での生成」という、 保存則の積分形そのものである。未知数として CV 内の平均濃度
$$ \bar{c}_i \equiv \frac{1}{V_i}\int_{V_i} c\, \mathrm{d}V \quad [\text{mol/m}^3] $$を採用すると、近似が必要なのは境界面上のフラックスの値だけになる。 面のフラックスを隣り合うセル平均値の差分で近似する — これが FVM の全てである。
FDM と FVM の違い
- FDM:点の値 $c_i \approx c(x_i)$ が未知数。微分演算子を差分で置き換える。 離散化後の式を全格子点で足しても、総量の保存が厳密に成り立つ保証はない (特に係数が場所によって変わる場合や不等間隔グリッドで崩れやすい)。
- FVM:セル平均 $\bar{c}_i$ が未知数。保存則の積分形を離散化する。 ある面から出たフラックスは必ず隣のセルに入るため、 全セルを足すと内部の面のフラックスがすべて相殺し(望遠鏡和)、 総量の変化は領域の外側境界のフラックスだけで決まる。 つまり丸め誤差を除いて総量が厳密に保存する。
電池シミュレーションでは総リチウム量のわずかなドリフトが容量・電圧の系統誤差として蓄積するので、 この性質は実用上きわめて重要である。
4.2.2 球座標 FVM の導出 — 粒子内拡散を離散化する
DFN の固相拡散方程式(フェーズ1の式 1)を離散化しよう。 各位置 $x$ の活物質粒子(半径 $R_s$ [m])の内部で、 固相 Li 濃度 $c_s(r,t)$ [mol/m³] は球対称の拡散方程式に従う:
$$ \frac{\partial c_s}{\partial t} = \frac{1}{r^2}\frac{\partial}{\partial r}\!\left(D_s r^2 \frac{\partial c_s}{\partial r}\right), \qquad \left.\frac{\partial c_s}{\partial r}\right|_{r=0} = 0, \qquad \left.-D_s \frac{\partial c_s}{\partial r}\right|_{r=R_s} = \frac{j}{F}. $$粒子を玉ねぎのように $N_r$ 枚の球殻(シェル)に分割する。 シェルの境界(面)を等間隔 $\Delta r = R_s/N_r$ [m] に置き、 $i$ 番目のセル($i = 1, \dots, N_r$)の内側の面を $r_{i-1/2} = (i-1)\Delta r$、 外側の面を $r_{i+1/2} = i\,\Delta r$ とする。セルの代表点(中心)は $r_i = (i-\tfrac12)\Delta r$ である。
導出:球殻セルの収支式
セルの体積と面の面積。セル $i$ は 2 つの同心球面に挟まれた球殻なので、その体積は
$$ V_i = \frac{4\pi}{3}\left(r_{i+1/2}^3 - r_{i-1/2}^3\right) \quad [\text{m}^3], $$また面 $r_{i+1/2}$ の面積は球面の面積で
$$ A_{i+1/2} = 4\pi r_{i+1/2}^2 \quad [\text{m}^2]. $$保存則の積分形をセル $i$ に適用する。球対称なのでフラックスは半径方向成分 $q_r = -D_s\, \partial c_s/\partial r$ のみを持つ。面 $r_{i+1/2}$ を通って外向き ($r$ 増加方向)に流れるモル流量(単位時間あたりのモル数)を
$$ W_{i+1/2} \equiv A_{i+1/2}\, q_r(r_{i+1/2}) \quad [\text{mol/s}] $$と書くと、セル内総モル数 $V_i \bar{c}_i$ の収支は「内側の面からの流入 − 外側の面からの流出」:
$$ V_i \frac{\mathrm{d} \bar{c}_i}{\mathrm{d} t} = W_{i-1/2} - W_{i+1/2}. $$生成項はない(反応は表面のみで、境界条件として入る)。この式もまだ厳密である。
面フラックスを差分で近似する(ここが唯一の近似)。 面 $r_{i+1/2}$ はセル中心 $r_i$ と $r_{i+1}$ のちょうど中点なので、 4.1 節の中心差分(2 次精度)がそのまま使える:
$$ q_r(r_{i+1/2}) \approx -D_s \frac{c_{i+1} - c_i}{\Delta r} \;\;\Longrightarrow\;\; W_{i+1/2} \approx -D_s A_{i+1/2} \frac{c_{i+1} - c_i}{\Delta r}. $$以降、セル平均 $\bar{c}_i$ を単に $c_i$ と書く(セル中心の値との差は $\mathcal{O}(\Delta r^2)$ であり、スキームの精度を落とさない)。
中心セル($i=1$)の扱い。内側の面は $r_{1/2} = 0$、つまり点である。 面積が $A_{1/2} = 4\pi \cdot 0^2 = 0$ なので $W_{1/2} = 0$ が自動的に成り立つ:
$$ V_1 \frac{\mathrm{d} c_1}{\mathrm{d} t} = -W_{3/2}. $$対称条件 $\partial c_s/\partial r|_{r=0}=0$ を別途課す必要はない。 FDM で球座標を扱うと $r=0$ で $\tfrac{1}{r^2}\partial_r(r^2 \partial_r c)$ の $1/r^2$ が発散して特別な処理(ロピタルの定理など)が要るのと対照的に、 FVM では特異点が現れない。これも FVM を選ぶ実用上の理由の一つである。
表面セル($i=N_r$)の扱い。外側の面 $r_{N_r+1/2} = R_s$ では 境界条件がフラックスの値そのものを与えている: $q_r(R_s) = -D_s\,\partial c_s/\partial r|_{R_s} = j/F$。よって差分近似は不要で、
$$ W_{N_r+1/2} = A_{N_r+1/2}\, \frac{j}{F} = 4\pi R_s^2\, \frac{j}{F} $$を厳密に代入できる。$j > 0$(酸化、放電時の負極)なら $W_{N_r+1/2} > 0$、すなわち Li が粒子から流出し、表面セルの濃度が下がる。符号もこれで整合する。
以上を組み立てると、球拡散の FVM 半離散式(時間は連続のまま)が完成する:
展開すると、$\mathrm{d}c_i/\mathrm{d}t$ は $c_{i-1}, c_i, c_{i+1}$ の 1 次結合になる。 つまり右辺は三重対角行列 $\mathbf{A}$ を使って $\mathrm{d}\boldsymbol{c}/\mathrm{d}t = \mathbf{A}\boldsymbol{c} + \boldsymbol{b}$ と書ける($\boldsymbol{b}$ は表面フラックス由来の定数ベクトル)。 なお共通因子 $4\pi$ は分母・分子で約分されるので、実装では省いてよい (本章の JS・Python コードでも省いている)。
物理的な意味:離散的な保存性とファラデーの法則
半離散式を全セルについて $V_i$ 倍して足すと、内部の面のモル流量 $W_{i+1/2}$ は 「セル $i$ からの流出」と「セル $i+1$ への流入」として 1 回ずつ逆符号で現れ、 厳密に相殺する(望遠鏡和)。残るのは両端だけ:
$$ \frac{\mathrm{d}}{\mathrm{d}t}\sum_{i=1}^{N_r} V_i c_i = W_{1/2} - W_{N_r+1/2} = -\,4\pi R_s^2\, \frac{j}{F}. $$粒子内の総 Li 量の変化率が、表面を流れる電流 $4\pi R_s^2 j$ をファラデー定数で割ったものに 丸め誤差の範囲で厳密に一致する — 離散版のファラデーの法則である。 実際、本章の実装でこの収支を数値確認すると相対誤差は $10^{-15}$ 程度(倍精度の丸め誤差)に収まる。 4.7 節でこれを検証項目の一つとして使う。
補足:電解液(x 方向)の FVM と、係数が変わる場合の面フラックス
液相の物質保存(フェーズ1の式 2)も同じ手順で離散化できる。1 次元 $x$ 方向では CV は厚さ $\Delta x$ の薄い板で、$V_i = A_{\mathrm{cell}}\Delta x$、面の面積は一定 $A_{\mathrm{cell}}$ となる($A_{\mathrm{cell}}$:電極面積 [m²])。 ただし DFN では実効拡散係数 $D_e^{\mathrm{eff}} = D_e \varepsilon_e^{b}$ が 領域(負極・セパレータ・正極)ごとに不連続に変わるため、 面上の拡散係数には両側のセルの値の調和平均を使うのが定石である (直列合成抵抗と同じ考え方。算術平均だとフラックスの連続性が崩れる)。 この詳細はフェーズ5で DFN を組み立てるときに具体的に扱う。
4.3 Method of Lines — 空間だけ離散化して ODE 系にする
4.2 節で得た半離散式は、時間についてはまだ連続の微分方程式である。 このように空間だけを先に離散化し、時間方向は常微分方程式(ODE)として残す アプローチを Method of Lines(MOL、線の方法)と呼ぶ。 各格子点の値 $c_i(t)$ が時間軸に沿って伸びる 1 本の「線」となることが名前の由来である。
全セルの未知数を 1 本のベクトル $\boldsymbol{y}(t) = (c_1, c_2, \dots, c_{N_r})^{\mathsf T}$ に積むと、半離散式は
という連立 ODE の初期値問題になる。ここまで来れば、
あとは既存の ODE ソルバー(solve_ivp など)に渡すだけでよい。
フェーズ2の SPM は、まさにこれだった
フェーズ2では、負極と正極の代表粒子それぞれについて球拡散を離散化し、
状態ベクトル(負極 $N_r$ 個 + 正極 $N_r$ 個の殻濃度)を solve_ivp に渡した。
あれは「FVM で空間離散化 → MOL で ODE 系化 → 汎用ソルバーで時間積分」
という本章の枠組みをそのまま実行していたのである。
当時ブラックボックスだったのは最後の「時間積分」の部分 — それを次節で開ける。
用語を整理しておく。空間だけ離散化した $\mathrm{d}\boldsymbol{y}/\mathrm{d}t = \boldsymbol{f}(t,\boldsymbol{y})$ を 半離散系(semi-discrete system)、 さらに時間も離散化して $\boldsymbol{y}^0, \boldsymbol{y}^1, \dots$ と時刻ごとの値の列にしたものを全離散系(fully discrete system)と呼ぶ。 以降、上付き添字 $m$ は時刻 $t_m = m\,\Delta t$ の値を表す ($\boldsymbol{y}^m \approx \boldsymbol{y}(t_m)$)。
4.4 時間積分と安定性 — 陽解法・陰解法と「硬さ」
4.4.1 前進オイラー法と後退オイラー法
最も単純な時間積分法はテイラー展開の 1 次打ち切りである。 $\boldsymbol{y}(t_m + \Delta t) = \boldsymbol{y}(t_m) + \Delta t\, \boldsymbol{y}'(t_m) + \mathcal{O}(\Delta t^2)$ において $\boldsymbol{y}' = \boldsymbol{f}$ を代入し、$\mathcal{O}(\Delta t^2)$ を捨てると
$$ \boldsymbol{y}^{m+1} = \boldsymbol{y}^m + \Delta t\, \boldsymbol{f}(t_m, \boldsymbol{y}^m) \qquad \text{(前進オイラー法、forward Euler)} $$右辺は既知の $\boldsymbol{y}^m$ だけで計算できるので、代入一発で次の時刻に進める。 このような方法を陽解法(explicit method)と呼ぶ。 1 ステップあたりの誤差(局所誤差)は捨てた項の $\mathcal{O}(\Delta t^2)$、 終端時刻まで $T/\Delta t$ 回のステップで誤差が積み上がるため、 全体誤差(大域誤差)は $\mathcal{O}(\Delta t)$ の1 次精度である。
一方、右辺の $\boldsymbol{f}$ をまだ知らない新しい時刻で評価すると
$$ \boldsymbol{y}^{m+1} = \boldsymbol{y}^m + \Delta t\, \boldsymbol{f}(t_{m+1}, \boldsymbol{y}^{m+1}) \qquad \text{(後退オイラー法、backward Euler)} $$となる。未知数 $\boldsymbol{y}^{m+1}$ が両辺に現れるので、 毎ステップ(一般には非線形の)方程式を解かなければならない。 このような方法を陰解法(implicit method)と呼ぶ。 精度は同じく 1 次。「毎ステップ方程式を解く」余分なコストを払ってまで 陰解法を使う理由が、次に述べる安定性である。
4.4.2 テスト方程式と絶対安定領域
安定性を調べる標準の道具がテスト方程式(test equation)である:
$$ \frac{\mathrm{d}y}{\mathrm{d}t} = \lambda y, \qquad \lambda \in \mathbb{C}. $$なぜこの単純な式で十分なのか。線形の半離散系 $\mathrm{d}\boldsymbol{y}/\mathrm{d}t = \mathbf{A}\boldsymbol{y}$ は、 $\mathbf{A}$ の固有ベクトルを基底に取り直すと、固有値 $\lambda_k$ ごとに独立な スカラー方程式 $y_k' = \lambda_k y_k$ に分解できるからである。 つまり「$\mathbf{A}$ のすべての固有値についてテスト方程式が安定に解ければ、 系全体も安定」となる。
導出:前進/後退オイラーの絶対安定領域
前進オイラーをテスト方程式に適用する。$f(y) = \lambda y$ なので
$$ y^{m+1} = y^m + \Delta t\, \lambda y^m = (1 + \lambda \Delta t)\, y^m . $$毎ステップ、解は一定の係数 $G = 1 + \lambda\Delta t$ 倍される。 $G$ を増幅率(amplification factor)[-] と呼ぶ。 $m$ ステップ後には $y^m = G^m y^0$ である。
厳密解が減衰する($\mathrm{Re}\,\lambda \le 0$)のに数値解が成長したら破綻である。 数値解が有界であり続ける条件は
$$ |G| = |1 + \lambda \Delta t| \le 1 . $$複素平面上で $z = \lambda\Delta t$ とみると、これは点 $-1$ を中心とする半径 1 の円板である。 この領域を前進オイラーの絶対安定領域(region of absolute stability)と呼ぶ。
実負の固有値 $\lambda < 0$(拡散)の場合、条件 $|1+\lambda\Delta t| \le 1$ は $-2 \le \lambda\Delta t \le 0$、すなわち
$$ \Delta t \le \frac{2}{|\lambda|} . $$固有値の絶対値が大きい(速く減衰する)モードほど、小さい $\Delta t$ を強制する。 後で見るように、拡散問題ではこの制限が非常に厳しい。
後退オイラーでは $y^{m+1} = y^m + \Delta t\,\lambda y^{m+1}$ を $y^{m+1}$ について解いて
$$ y^{m+1} = \frac{1}{1 - \lambda\Delta t}\, y^m \qquad\Longrightarrow\qquad G = \frac{1}{1 - \lambda\Delta t} . $$$\mathrm{Re}\,\lambda \le 0$ なら $|1 - \lambda\Delta t| \ge 1$ となるので、 どんなに大きな $\Delta t$ でも $|G| \le 1$ が成り立つ。 左半平面全体($\mathrm{Re}\,z \le 0$)を安定領域に含む方法を A-安定(A-stable)と呼ぶ。後退オイラーは A-安定である。
まとめ:陽解法 vs 陰解法
- 前進オイラー(陽):1 ステップが安い(代入のみ)。ただし $\Delta t \le 2/|\lambda_{\max}|$ の安定限界があり、 これを 1% でも超えると解は指数的に発散する。
- 後退オイラー(陰):毎ステップ方程式を解く必要があるが A-安定。$\Delta t$ は精度の要求だけで選べる。
4.4.3 von Neumann 安定性解析 — FTCS の安定限界を導く
今度は具体的な全離散スキームの安定性を調べる。対象は 1 次元拡散方程式 $\partial c/\partial t = D\, \partial^2 c/\partial x^2$($D$:拡散係数 [m²/s])を 「時間: 前進差分、空間: 中心差分」で離散化した FTCS(Forward Time, Centered Space)スキームである:
$$ c_\ell^{m+1} = c_\ell^m + \mu \left( c_{\ell+1}^m - 2 c_\ell^m + c_{\ell-1}^m \right), \qquad \mu \equiv \frac{D\,\Delta t}{\Delta x^2} . $$$\mu$:拡散数(diffusion number)[-] — 時間刻みと格子幅の関係を表す無次元数。
$G$:増幅率 [-]、$\beta$:モードの位相角 [rad] — 隣の格子点に移るごとに波が進む位相。 波数 $\kappa_w$ [rad/m] の正弦波モードなら $\beta = \kappa_w \Delta x$ である。
導出:von Neumann 安定性解析(1 ステップずつ)
誤差も同じ漸化式に従う。数値解 $c_\ell^m$ と、丸め誤差などが乗った解 $\tilde{c}_\ell^m$ の差を $\varepsilon_\ell^m = \tilde{c}_\ell^m - c_\ell^m$ とおく。 FTCS は線形なので、2 つの解の差 $\varepsilon$ 自身も同じ漸化式を満たす:
$$ \varepsilon_\ell^{m+1} = \varepsilon_\ell^m + \mu\left(\varepsilon_{\ell+1}^m - 2\varepsilon_\ell^m + \varepsilon_{\ell-1}^m\right). $$「誤差が増幅されるか減衰するか」を調べればよい。
誤差をフーリエモードに分解する。周期境界のもとで、任意の誤差分布は 正弦波モードの重ね合わせで書ける。線形なのでモードごとに独立に追跡できる。 1 つのモードを
$$ \varepsilon_\ell^m = G^m e^{\mathrm{i} \beta \ell} $$とおく($\mathrm{i}$:虚数単位、$G^m$ は $G$ の $m$ 乗)。 $\beta \in [-\pi, \pi]$ が波の細かさを表し、$\beta = \pi$ が 格子で表現できる最も細かい波(1 点ごとに符号が反転するジグザグ)である。
漸化式に代入する。$\varepsilon_{\ell\pm1}^m = G^m e^{\mathrm{i}\beta \ell} e^{\pm\mathrm{i}\beta}$ なので
$$ G^{m+1} e^{\mathrm{i}\beta\ell} = G^m e^{\mathrm{i}\beta\ell}\left[ 1 + \mu\left( e^{\mathrm{i}\beta} - 2 + e^{-\mathrm{i}\beta} \right) \right]. $$両辺を $G^m e^{\mathrm{i}\beta\ell}$ で割ると、増幅率が求まる:
$$ G = 1 + \mu\left(e^{\mathrm{i}\beta} + e^{-\mathrm{i}\beta} - 2\right) = 1 + 2\mu\left(\cos\beta - 1\right). $$半角公式で整理する。$1 - \cos\beta = 2\sin^2(\beta/2)$ を使うと
$$ G = 1 - 4\mu \sin^2\!\frac{\beta}{2} . $$$G$ は実数で、$\mu > 0$ より常に $G \le 1$。問題は下限、 つまり $G < -1$ になって符号を交互に変えながら成長する(振動発散)可能性である。
すべてのモードで $|G| \le 1$ を要求する。$G \ge -1$ すなわち $4\mu\sin^2(\beta/2) \le 2$ が任意の $\beta$ で成り立つ必要がある。 左辺は $\beta = \pi$(最も細かいジグザグ波)で最大値 $4\mu$ をとるので、条件は $4\mu \le 2$:
$$ \mu \le \frac{1}{2} \qquad\Longleftrightarrow\qquad \boxed{\;\Delta t \le \frac{\Delta x^2}{2D}\;} $$この結果はテスト方程式の見方とも整合する。実は FTCS は 「MOL(空間中心差分)+ 前進オイラー」に他ならず、空間差分行列 $\mathbf{A}$ の固有値は $\lambda_\beta = -\frac{4D}{\Delta x^2}\sin^2(\beta/2)$、その最大絶対値は $|\lambda_{\max}| = 4D/\Delta x^2$(ジグザグモード)である。 前進オイラーの条件 $\Delta t \le 2/|\lambda_{\max}|$ に代入すると、 ちょうど $\Delta t \le \Delta x^2/(2D)$ が再現される。
4.4.4 拡散問題は「硬い」— stiffness
ここで固有値の分布に注目する。拡散の離散化行列の固有値は
- 最も遅いモード(粒子全体にわたるなだらかな濃度分布): $|\lambda| \sim D/R_s^2$ 程度 — 物理的に興味のある時間スケール(拡散時定数 $R_s^2/D_s$、 負極粒子なら $ (5.86\times10^{-6})^2 / 3.3\times10^{-14} \approx 1000$ s)
- 最も速いモード(格子幅のジグザグ):$|\lambda| \approx 4D/\Delta x^2 = 4N_r^2 D/R_s^2$ — 格子を細かくするほど2 乗で大きくなる
と、およそ $(2N_r/\pi)^2$ 倍($N_r = 40$ なら約 650 倍)も広がっている。 このように速い時定数と遅い時定数が混在し、安定性の制約(速いモード)が、 知りたい現象(遅いモード)に必要な精度よりはるかに小さい $\Delta t$ を強制する系を 硬い(stiff)系と呼ぶ。
数字で確かめよう。負極粒子($D_s = 3.3\times10^{-14}$ m²/s、$R_s = 5.86\times10^{-6}$ m)を $N_r = 40$ 分割すると $\Delta r = 1.47\times10^{-7}$ m、FTCS の安定限界は
$$ \Delta t \le \frac{\Delta r^2}{2 D_s} = \frac{(1.47\times10^{-7})^2}{2 \times 3.3\times10^{-14}} \approx 0.33 \text{ s}. $$1C 放電(3600 s)を解くには1 万ステップ以上が必要になる。 しかも濃度分布はゆっくりとしか変化しないので、精度のためなら $\Delta t = 10$ s でも十分なのに、である。 $N_r$ を 80 にすれば限界はさらに 1/4 の 0.08 s になる。 グリッドを細かくするほど陽解法はどんどん割に合わなくなる — これが 拡散型 PDE で陰解法が標準になっている理由である。
球拡散(負極粒子、1C 相当の一定流束で Li が流出)を前進オイラーで解く。 スライダーで $\Delta t$ を安定限界に対する比で調整し、▶ 再生で時間発展を見る。 比が 1.00 以下ならなだらかに解けるが、1.01 でもわずかに超えると ジグザグ振動が成長してやがて発散する(発散すると赤い注記が出る。軸は固定のまま)。 $N_r$ を増やすと安定限界の $\Delta t$ 自体が小さくなる(表示値に注目)ことも確認してほしい。
4.4.5 高次の方法:Crank–Nicolson と BDF
後退オイラーは A-安定だが 1 次精度しかない。実用では次の 2 系統がよく使われる。
Crank–Nicolson 法(台形則)
前進と後退の平均、すなわち右辺を区間の両端で平均する台形則である:
$$ \boldsymbol{y}^{m+1} = \boldsymbol{y}^m + \frac{\Delta t}{2}\left[ \boldsymbol{f}(t_m, \boldsymbol{y}^m) + \boldsymbol{f}(t_{m+1}, \boldsymbol{y}^{m+1}) \right]. $$テスト方程式での増幅率は $G = \dfrac{1 + \lambda\Delta t/2}{1 - \lambda\Delta t/2}$ で、 $\mathrm{Re}\,\lambda \le 0$ なら常に $|G| \le 1$(A-安定)、かつ 局所誤差 $\mathcal{O}(\Delta t^3)$・大域誤差 $\mathcal{O}(\Delta t^2)$ の 2 次精度である。線形問題なら毎ステップ三重対角の連立一次方程式を 解くだけでよく(Thomas 法、計算量 $\mathcal{O}(N)$)、本章の図4.2 や フェーズ2・3のデモの内部でもこの方法を使っている。
注意:Crank–Nicolson の「減衰しない」振動
CN の増幅率は $\lambda\Delta t \to -\infty$ の極限で $G \to -1$ に近づく。 発散はしないが、速いモードがほとんど減衰せずに符号だけ反転し続けるため、 初期条件が急峻(例:電流を階段状に印加した直後)だと解に人工的な振動が残ることがある。 $|\lambda\Delta t| \to \infty$ で $G \to 0$ となる性質(L-安定)を持つ 後退オイラーや後述の BDF は、この点で頑健である。
BDF(後退差分公式)
BDF(backward differentiation formula)は、過去の複数ステップの値で $\boldsymbol{y}'(t_{m+1})$ を後退差分近似し、それを $\boldsymbol{f}(\boldsymbol{y}^{m+1})$ と等置する多段の陰解法である。1 段の BDF1 は後退オイラーそのもの。2 段の BDF2 は
$$ \frac{3\boldsymbol{y}^{m+1} - 4\boldsymbol{y}^m + \boldsymbol{y}^{m-1}}{2\Delta t} = \boldsymbol{f}(t_{m+1}, \boldsymbol{y}^{m+1}) $$で 2 次精度かつ A-安定・L-安定。BDF3〜BDF6 はさらに高次だが、 安定領域が左半平面より少し欠ける(A($\alpha$)-安定)。 scipy や SUNDIALS の BDF ソルバーは、誤差推定に基づいて 次数(1〜5)と $\Delta t$ を自動調整しながら進む。
scipy の solve_ivp ではどれを選ぶか
| method | 種別 | 精度 | 向いている問題 |
|---|---|---|---|
RK45(既定) | 陽的 Runge–Kutta | 4(5)次 | 硬くない問題。1 ステップが安い。硬い問題に使うと安定限界に縛られ ステップ数が爆発する(コード1で実測する)。 |
BDF | 陰的多段(BDF1–5) | 可変次数 | 硬い問題の第一候補。ヤコビアン(または疎パターン)を渡すと速い。 |
Radau | 陰的 Runge–Kutta(Radau IIA) | 5 次 | 硬い問題で高精度が欲しいとき。L-安定で振動にも頑健。1 ステップは BDF より重い。 |
LSODA | 自動切替(Adams/BDF) | 可変次数 | 硬さが不明・途中で変わる問題。まず試す価値がある。 |
実際に「硬さ」を体感しよう。4.2 節の球拡散 FVM(これはフェーズ2の SPM の粒子と同じ問題)を
RK45 と BDF で解き、右辺の評価回数
nfev・ステップ数・所要時間を比べる。
# コード1: 球拡散 FVM を RK45 と BDF で解いて「硬さ」を体感する
import numpy as np, time
from scipy.integrate import solve_ivp
from scipy.sparse import diags
# --- パラメータ(LG M50 負極粒子、SPEC の値と同一)---
F = 96485.33212 # ファラデー定数 [C/mol]
Ds = 3.3e-14 # 固相拡散係数 [m^2/s]
Rs = 5.86e-6 # 粒子半径 [m]
cs_max = 33133.0 # 最大 Li 濃度 [mol/m^3]
I_1C = 48.69 # 1C 電流密度 [A/m^2]
a_n = 3*0.75/Rs # 比界面積 [1/m]
L_n = 85.2e-6 # 負極厚さ [m]
j = I_1C/(a_n*L_n) # 1C 相当の界面電流密度 [A/m^2]
# --- 球座標 FVM(4.2 節の式をそのままコードに)---
Nr = 40
r_f = np.linspace(0.0, Rs, Nr+1) # セル面の半径 r_{i±1/2}
V = (r_f[1:]**3 - r_f[:-1]**3)/3.0 # シェル体積(共通因子 4π は約分)
A = r_f**2 # 面の面積(4π は約分)
w_in = Ds*A[1:-1]/((Rs/Nr)*V[1:]) # 内側面からの結合係数
w_out = Ds*A[1:-1]/((Rs/Nr)*V[:-1]) # 外側面からの結合係数
diag = np.zeros(Nr); diag[:-1] -= w_out; diag[1:] -= w_in
Amat = diags([w_in, diag, w_out], [-1, 0, 1], format="csr")
b = np.zeros(Nr); b[-1] = -(j/F)*A[-1]/V[-1] # 表面 BC: -Ds dc/dr = j/F
def rhs(t, c):
return Amat @ c + b
c0 = np.full(Nr, 0.5*cs_max) # 初期: 一様(θ = 0.5)
t_end = 1800.0 # 30 分の放電
print(f"陽解法の安定限界の目安 Δt ≈ Δr²/(2Ds) = {(Rs/Nr)**2/(2*Ds):.3f} s")
print(f"必要ステップ数の目安 = {t_end/((Rs/Nr)**2/(2*Ds)):.0f}\n")
print(f"{'手法':6s} {'ステップ数':>8s} {'nfev':>8s} {'njev':>6s} {'時間 [s]':>9s}")
results = {}
for method in ["RK45", "BDF"]:
# BDF には定数ヤコビアン(= Amat そのもの)を渡す。行列を直接渡すと
# 評価が不要になるため njev は 0 と表示される。
kw = {"jac": Amat} if method == "BDF" else {}
t0 = time.perf_counter()
sol = solve_ivp(rhs, (0, t_end), c0, method=method,
rtol=1e-6, atol=1e-3, **kw)
el = time.perf_counter() - t0
results[method] = sol
print(f"{method:6s} {len(sol.t)-1:8d} {sol.nfev:8d} {sol.njev:6d} {el:9.3f}")
# 2 つの解が一致することも確認
d = np.max(np.abs(results["RK45"].y[:, -1] - results["BDF"].y[:, -1]))/cs_max
print(f"\n最終プロファイルの差 max|Δθ| = {d:.2e}(両者は同じ解に収束)")
r_c = (0.5*(r_f[:-1] + r_f[1:])/Rs).tolist()
plot_spec = {
"traces": [
{"x": r_c, "y": (results["RK45"].y[:, -1]/cs_max).tolist(),
"name": "RK45(陽的)", "mode": "lines"},
{"x": r_c, "y": (results["BDF"].y[:, -1]/cs_max).tolist(),
"name": "BDF(陰的)", "mode": "lines", "line": {"dash": "dash"}},
],
"layout": {"title": {"text": "t = 1800 s の濃度プロファイル(両手法は一致)"},
"xaxis": {"title": {"text": "r / R_s [-]"}},
"yaxis": {"title": {"text": "θ = c_s / c_s,max [-]"}}},
}
実行すると、RK45 は 3000 ステップ超・nfev 2 万回超を要するのに対し、
BDF は 100 ステップ弱で同じ解(最終プロファイルの差は $10^{-7}$ 以下)に到達する。
RK45 のステップ数は要求精度を緩めてもほとんど減らない — 精度ではなく
安定限界(RK45 の安定領域は実軸上で $|\lambda|\Delta t \lesssim 3.3$)が
刻み幅を決めているからで、実測値は $|\lambda_{\max}| T / 3.3 \approx 3400$ という見積もりとよく一致する。
これが「硬い問題に陽解法を使ってはいけない」の定量的な意味である。
4.5 DAE の発生 — 微分変数と代数変数
4.5.1 DFN の変数を数え上げる
ここまでの道具(FVM + MOL + 陰解法)で SPM は完全に解ける。 しかし DFN には、これだけでは扱えない構造上の特徴がある。 フェーズ1で導いた支配方程式を、時間微分を持つかどうかで仕分けしてみよう。
| 変数 | 支配方程式 | $\partial/\partial t$ は? | 分類 |
|---|---|---|---|
| $c_s(x,r,t)$ 固相濃度 | 固相拡散(式 1) | ある | 微分変数 |
| $c_e(x,t)$ 電解液濃度 | 液相物質保存(式 2) | ある($\varepsilon_e\,\partial c_e/\partial t$) | 微分変数 |
| $\phi_s(x,t)$ 固相電位 | 固相電荷保存(式 3) | ない | 代数変数 |
| $\phi_e(x,t)$ 液相電位 | 液相電荷保存(式 4) | ない | 代数変数 |
| $j(x,t)$ 界面電流密度 | Butler–Volmer(式 5) | ない | 代数変数 |
電位の式 3・4 は $x$ についての 2 階微分を含むだけの楕円型方程式、 Butler–Volmer(式 5)は微分すら含まない純粋な代数式である。 つまり $\phi_s, \phi_e, j$ には「時間発展の式」がなく、 各時刻において、その瞬間の $c_s, c_e$ と拘束条件を満たすように瞬時に決まる。
物理的な意味:なぜ電位には時間微分がないのか
電荷の緩和(電気二重層の充放電)の時定数はナノ〜マイクロ秒のオーダーで、 拡散(数百〜数千秒)より桁違いに速い。DFN では二重層容量を無視する近似 (フェーズ1の仮定)を置いたため、電位場は「無限に速い」変数となり、 時間微分項が消えて拘束条件に退化した。 これは 4.4 節の言葉でいえば「極端に硬い成分を極限まで硬くして、 代数方程式に置き換えた」ことに相当する。 硬さがモデルの近似(準静的電荷保存)として先取りされているのである。
空間離散化後の未知数の個数を具体的に数える。負極・セパレータ・正極の $x$ 方向ノード数を $N_n = 20,\ N_s = 10,\ N_p = 20$、各粒子の径方向シェル数を $N_r = 30$ とすると:
| 変数 | 存在する場所 | 個数 | 分類 |
|---|---|---|---|
| $c_s$ | 電極内の各 $x$ ノードに粒子 1 個 × $N_r$ シェル | $(20+20)\times 30 = 1200$ | 微分 |
| $c_e$ | 全領域の $x$ ノード | $20+10+20 = 50$ | 微分 |
| $\phi_s$ | 電極内の $x$ ノード | $20+20 = 40$ | 代数 |
| $\phi_e$ | 全領域の $x$ ノード | $50$ | 代数 |
| $j$ | 電極内の $x$ ノード | $40$ | 代数 |
| 合計 | $1380$ | 微分 1250 + 代数 130 |
4.5.2 質量行列が特異な系 = 微分代数方程式(DAE)
全 1380 個の未知数を 1 本のベクトル $\boldsymbol{y}$ に積んで方程式系を書くと、 $c_s, c_e$ の行には $\mathrm{d}/\mathrm{d}t$ が付き、$\phi_s, \phi_e, j$ の行には付かない。 これをまとめて書くには、時間微分の係数行列 $\mathbf{M}$(質量行列、mass matrix)を導入すればよい:
$\mathbf{M}$ は対角に 0 が並ぶ特異行列(逆行列が存在しない)なので、 両辺に $\mathbf{M}^{-1}$ を掛けて普通の ODE $\mathrm{d}\boldsymbol{y}/\mathrm{d}t = \mathbf{M}^{-1}\boldsymbol{f}$ に直す、という手が使えない。 このような系を微分代数方程式(differential-algebraic equations, DAE)と呼ぶ。 微分変数を $\boldsymbol{y}_d$($c_s, c_e$)、代数変数を $\boldsymbol{y}_a$($\phi_s, \phi_e, j$)と 分けて書くと、DFN は次の半陽的(semi-explicit)DAEである:
注意:solve_ivp は DAE を解けない
scipy の solve_ivp は $\mathrm{d}\boldsymbol{y}/\mathrm{d}t = \boldsymbol{f}$
の形しか受け付けず、特異な質量行列を渡すインターフェースがない。
フェーズ2・3で solve_ivp が使えたのは、SPM・SPMe では代数変数を
手で消去できた(電流一定なら $j$ が陽に決まり、電位は後処理で計算できた)からである。
DFN では $j(x)$ の分布そのものが $\phi_s, \phi_e$ と連立して初めて決まるため、
この消去ができない。だから DAE を DAE のまま解く方法(4.6 節)が必要になる。
4.5.3 指数(index)— この問題は index-1
DAE の「解きにくさ」を測る量が微分指数(differentiation index)である。 直感的には、代数拘束を時間微分していくと ODE 系に変換できるが、 そのために必要な微分の回数が指数である。
半陽的 DAE で拘束 $\boldsymbol{g}(\boldsymbol{y}_d, \boldsymbol{y}_a) = \boldsymbol{0}$ を 1 回時間微分してみる。連鎖律で
$$ \boldsymbol{0} = \frac{\mathrm{d}}{\mathrm{d}t}\boldsymbol{g} = \frac{\partial \boldsymbol{g}}{\partial \boldsymbol{y}_d} \underbrace{\frac{\mathrm{d}\boldsymbol{y}_d}{\mathrm{d}t}}_{=\,\boldsymbol{f}_d\ (\text{既知の式})} + \frac{\partial \boldsymbol{g}}{\partial \boldsymbol{y}_a} \frac{\mathrm{d}\boldsymbol{y}_a}{\mathrm{d}t}. $$もし行列 $\partial \boldsymbol{g}/\partial \boldsymbol{y}_a$ が正則(逆行列を持つ)なら、 これを $\mathrm{d}\boldsymbol{y}_a/\mathrm{d}t$ について解けて、系全体が ODE になる。 微分 1 回で済んだのでindex-1 である。 正則でなければさらに微分が必要で、index-2 以上の「高指数」DAE となり、 数値解法は格段に難しくなる(拘束のドリフト、初期化の複雑化など)。
ミニ例で確かめよう(このあとコード2で実際に解く例と同じもの)。 $y' = -(z + z^3)$、$0 = g(y,z) = z + z^3 - y$ という系では $\partial g/\partial z = 1 + 3z^2 > 0$ が常に成り立つので index-1 であり、 $g$ を 1 回微分すれば $z' = y'/(1+3z^2)$ という ODE が得られる。
DFN の場合、$\partial \boldsymbol{g}/\partial \boldsymbol{y}_a$ は 「$c_s, c_e$ を凍結したときの、$(\phi_s, \phi_e, j)$ についての回路方程式のヤコビアン」である。 Butler–Volmer は $\eta = \phi_s - \phi_e - U(\theta)$ について単調増加 ($\partial j/\partial \eta > 0$、指数関数の和なので)、電位方程式は 伝導率が正の拡散型演算子であり、電位の基準(接地)を 1 点固定すれば この代数系は一意に解ける = ヤコビアンは正則。 したがって DFN の空間離散化系は index-1 の DAE である。 これは幸運な性質で、index-1 なら標準的な BDF 型 DAE ソルバー(後述の IDA など)が そのまま適用でき、次節の 3 つの戦略もすべて機能する。
4.5.4 整合的な初期条件(consistent initialization)
ODE の初期値問題では $\boldsymbol{y}(0)$ を自由に与えられた。DAE ではそうはいかない。 初期値は代数拘束を満たさなければならない:
$$ \boldsymbol{g}\big(\boldsymbol{y}_d(0),\, \boldsymbol{y}_a(0)\big) = \boldsymbol{0}. $$つまり微分変数の初期値($c_s, c_e$ の初期分布)は自由に選べるが、 代数変数の初期値($\phi_s, \phi_e, j$)は拘束を解いて求めるしかない。 これを整合初期化(consistent initialization)と呼ぶ。 DFN の放電開始($t=0$ で電流 $I$ を階段状に印加)では:
- $c_s(x,r,0)$ は一様(100% SOC の $\theta_{n,100}, \theta_{p,100}$)、$c_e(x,0) = c_{e0}$ と与える。
- その $c_s, c_e$ を凍結したまま、式 3・4・5 の連立非線形方程式を ニュートン法で解いて $\phi_s(x,0), \phi_e(x,0), j(x,0)$ を求める。
印加電流を $0$ から $I$ に切り替えた瞬間、微分変数は連続のまま 代数変数($j$ の分布や端子電圧)だけが不連続にジャンプする。 「電流を流した瞬間に電圧が IR 降下と反応過電圧のぶんだけ即座に下がる」という フェーズ0以来の物理観察は、DAE の言葉では「代数変数のジャンプ」なのである。 整合初期化を怠る(たとえば $\phi, j$ を全部 0 で始める)と、 陰解法の最初のステップでニュートン法が収束しないか、非物理的な過渡が混入する。
4.6 DAE の解き方 — 3つの戦略とニュートン法
index-1 の半陽的 DAE を解く実用的な戦略は、大きく 3 つある。
4.6.1 戦略A:入れ子方式(nested approach)
時間積分は微分変数だけについて行い、 右辺 $\boldsymbol{f}_d(\boldsymbol{y}_d, \boldsymbol{y}_a)$ を評価するたびに、 その場で代数拘束 $\boldsymbol{g}(\boldsymbol{y}_d, \boldsymbol{y}_a) = \boldsymbol{0}$ を ニュートン法で解いて $\boldsymbol{y}_a$ を求める:
$$ \boldsymbol{y}_a = \boldsymbol{y}_a(\boldsymbol{y}_d) \quad\text{(拘束を解いて陰的に定まる関数とみなす)} \qquad\Longrightarrow\qquad \frac{\mathrm{d}\boldsymbol{y}_d}{\mathrm{d}t} = \boldsymbol{f}_d\big(\boldsymbol{y}_d,\, \boldsymbol{y}_a(\boldsymbol{y}_d)\big). $$
index-1 なら($\partial\boldsymbol{g}/\partial\boldsymbol{y}_a$ が正則なので)陰関数定理により
$\boldsymbol{y}_a(\boldsymbol{y}_d)$ は局所的に一意に決まり、この置き換えは正当である。
外側の時間積分には solve_ivp の BDF でも何でも使える。
実装が素直で理解しやすく、フェーズ5の教材コードの出発点もこの形にする。
欠点は、右辺評価 1 回ごとに内側のニュートン反復が走るためコストが高いこと、
内側の収束が甘いと外側の誤差制御が乱れることである。
4.6.2 戦略B:一体方式(monolithic approach)
時間積分そのものを陰解法(たとえば後退オイラー)にするなら、 どうせ毎ステップ非線形方程式を解く。それなら微分変数と代数変数を まとめて 1 つの非線形方程式にしてしまえばよい。 後退オイラーを $\mathbf{M}\,\mathrm{d}\boldsymbol{y}/\mathrm{d}t = \boldsymbol{f}$ に適用すると、 ステップ $k{+}1$ の未知数 $\boldsymbol{y}^{k+1}$ についての方程式
が得られる。微分変数の行は「後退オイラーの更新式」、代数変数の行は ($\mathbf{M}$ の該当行が 0 なので)「その時刻での拘束 $-\boldsymbol{g} = \boldsymbol{0}$」に 自動的になる。これを毎ステップ、ニュートン法で解く。 代数拘束は毎ステップ厳密に(ニュートン法の許容誤差まで)満たされ、 構造も見通しがよい。フェーズ5の完成版 DFN ソルバーはこの方式で実装する。
4.6.3 戦略C:専用 DAE ソルバー
戦略Bの「可変次数 BDF + 誤差制御 + ニュートン法」を高品質に実装した ライブラリを使う手もある。事実上の標準は SUNDIALS の IDA (Fortran 時代の DASSL の後継。$F(t, \boldsymbol{y}, \boldsymbol{y}') = 0$ 形式の index-1 DAE を可変次数 BDF で解く)。 PyBaMM は、モデルの式を CasADi(自動微分ライブラリ)で記号的に構築して 厳密なヤコビアンを生成し、IDAKLU(IDA + 疎行列 LU 分解ライブラリ KLU)で 解いている。フェーズ5で PyBaMM を参照解として使うとき、内部ではこれが動いている。
| A: 入れ子 | B: 一体 | C: 専用ソルバー | |
|---|---|---|---|
| 実装の手間 | 小(ODE ソルバー流用) | 中(Newton を自作) | 小(ただし導入・依存関係) |
| 計算効率 | 低い(右辺評価ごとに内側 Newton) | 良い(1 ステップ 1 Newton) | 最良(可変次数・最適化済み) |
| 制御・学習効果 | 中身が見える | 中身が全部見える | ブラックボックス |
| 本教材での役割 | フェーズ5の第一歩 | フェーズ5の完成形 | 検証相手(PyBaMM) |
4.6.4 ニュートン法の復習 — 導出・収束・失敗モード
どの戦略でも心臓部はニュートン法(Newton's method / Newton–Raphson 法)である。 非線形方程式 $\boldsymbol{F}(\boldsymbol{y}) = \boldsymbol{0}$($\boldsymbol{F}: \mathbb{R}^N \to \mathbb{R}^N$)を解きたい。
導出:テイラー展開の 1 次打ち切り
現在の近似解 $\boldsymbol{y}^{(k)}$ のまわりで、求める解 $\boldsymbol{y}^\ast$ ($\boldsymbol{F}(\boldsymbol{y}^\ast) = \boldsymbol{0}$)に向けてテイラー展開する。 修正量を $\boldsymbol{\delta} = \boldsymbol{y}^\ast - \boldsymbol{y}^{(k)}$ とおくと
$$ \boldsymbol{0} = \boldsymbol{F}(\boldsymbol{y}^\ast) = \boldsymbol{F}(\boldsymbol{y}^{(k)}) + \mathbf{J}(\boldsymbol{y}^{(k)})\,\boldsymbol{\delta} + \mathcal{O}(\|\boldsymbol{\delta}\|^2), $$ここで $\mathbf{J}$ はヤコビアン(Jacobian)行列、 成分は $J_{pq} = \partial F_p/\partial y_q$ である($N \times N$ 行列)。
2 次の項を捨てて $\boldsymbol{\delta}$ について解く。これがニュートン法の 1 反復:
$$ \mathbf{J}(\boldsymbol{y}^{(k)})\, \boldsymbol{\delta}^{(k)} = -\boldsymbol{F}(\boldsymbol{y}^{(k)}), \qquad \boldsymbol{y}^{(k+1)} = \boldsymbol{y}^{(k)} + \boldsymbol{\delta}^{(k)}. $$毎反復、連立一次方程式を 1 回解くのが主なコストである($N=1$ なら おなじみの $y^{(k+1)} = y^{(k)} - F/F'$、幾何学的には接線と $x$ 軸の交点)。
収束の速さ。捨てたのが 2 次の項なので、誤差 $\boldsymbol{e}^{(k)} = \boldsymbol{y}^{(k)} - \boldsymbol{y}^\ast$ は
$$ \|\boldsymbol{e}^{(k+1)}\| \le C\, \|\boldsymbol{e}^{(k)}\|^2 \qquad (C \text{ は } \mathbf{J} \text{ の変化率で決まる定数}) $$を満たす(解の近くで $\mathbf{J}$ が正則かつリプシッツ連続な場合)。 2 次収束(quadratic convergence)である。 誤差が $10^{-2} \to 10^{-4} \to 10^{-8} \to 10^{-16}$ と、 正しい桁数が反復ごとに倍増し、数回で倍精度の限界に達する。
失敗モードも知っておく必要がある:
- 初期値が遠い:テイラー展開の 1 次近似が成り立たない場所から始めると、 接線が見当違いの方向に飛ばし、振動・発散することがある。 対策は減衰(damping / line search: $\boldsymbol{y}^{(k+1)} = \boldsymbol{y}^{(k)} + \alpha\boldsymbol{\delta}$、 $0 < \alpha \le 1$ を残差が減るまで縮める)や、良い初期値の工夫。 時間発展計算では直前の時刻の解が極上の初期値になるため、 実際のステップ内 Newton は 2〜4 反復で収束するのが普通である。
- ヤコビアンが特異・悪条件:$\mathbf{J}$ が(ほぼ)特異だと $\boldsymbol{\delta}$ が求まらない、または巨大になる。 DFN では電位の基準を固定し忘れる($\phi$ 全体を定数シフトしても方程式が変わらない)と まさに特異になる。接地条件を 1 つ入れるのが正しい対処である。
- 非滑らかな $\boldsymbol{F}$:if 文による切替や $\sqrt{\cdot}$ の定義域外など、
微分が不連続な点をまたぐと 2 次収束が壊れる。$i_0$ の $\sqrt{c}$ には
下限クリップを入れる(params.js の
i0_funが実際にそうしている)。
図4.3 で、DFN の心臓部である Butler–Volmer 方程式を例にニュートン法を 1 反復ずつ動かしてみよう。 界面電流密度の目標値 $j_{\mathrm{target}}$ [A/m²] を与えて、それを実現する過電圧 $\eta$ [V] を求める:
$$ g(\eta) = j_{\mathrm{target}} - i_0\!\left[ \exp\!\left(\frac{F\eta}{2RT}\right) - \exp\!\left(-\frac{F\eta}{2RT}\right) \right] = j_{\mathrm{target}} - 2 i_0 \sinh\!\frac{F\eta}{2RT} = 0 $$($\alpha_a = \alpha_c = 0.5$ の Butler–Volmer。$i_0$:交換電流密度 [A/m²]、 $R$:気体定数 [J/(mol·K)]、$T$:温度 [K])。
「1 ステップ進める」を押すたびにニュートン法が 1 反復進む。曲線上の点から接線(橙)を引き、 $g=0$ との交点が次の反復点になる。表の誤差 $e_k$ の指数に注目 — 解に近づくと指数が反復ごとにほぼ倍増し(2 次収束)、 $e_k/e_{k-1}^2$ がほぼ一定になる。初期値を遠くに置くと、最初のうちは ほぼ一定幅(約 $2RT/F \approx 51$ mV)でしか進まず、2 次収束が「効き始める」までの 様子も観察できる。★ は厳密解 $\eta^\ast = \frac{2RT}{F}\,\mathrm{asinh}\,\frac{j_{\mathrm{target}}}{2 i_0}$。
4.6.5 疎なヤコビアン — 使わなければ $\mathcal{O}(N^3)$ で破綻する
ニュートン法の 1 反復は連立一次方程式 $\mathbf{J}\boldsymbol{\delta} = -\boldsymbol{F}$ を解くことだった。 $N \times N$ の密行列の LU 分解は $\mathcal{O}(N^3)$ の計算量である。 4.5 節の控えめなグリッドでも $N = 1380$、つまり 1 回の分解に約 $\tfrac{2}{3}N^3 \approx 1.7 \times 10^9$ 回の浮動小数点演算。 1 回なら耐えられるが、時間ステップ数 × ステップあたり Newton 反復数(数千回)を掛けると $10^{12}$ 回超 — 数分〜数十分待たされる計算になり、 グリッドを倍にすれば($N^3$ なので)8 倍に膨らむ。破綻である。
救いは、$\mathbf{J}$ の成分のほとんどがゼロという事実にある。 $J_{pq} \ne 0$ となるのは「方程式 $p$ が変数 $q$ に直接依存する」ときだけで、 離散化された方程式はどれも局所的(自分と隣のセル、同じ場所の他の変数しか参照しない)だからである。 DFN のヤコビアンの構造をブロックごとに整理すると:
- $c_s$ ブロック:各粒子の拡散は $r$ 方向に三重対角。 粒子同士は $x$ 方向に直接結合しない(隣の粒子の濃度は式に現れない)ので、 「小さな三重対角が粒子の数だけ並んだブロック対角」になる。 表面シェルの行だけが、境界条件を通じて同じ場所の $j$ の列に非ゼロを持つ。
- $c_e$, $\phi_e$, $\phi_s$ ブロック:$x$ 方向の 2 階差分なので三重対角の帯。 $\phi_e$ の式は濃度過電位項($\partial \ln c_e/\partial x$)を通じて $c_e$ の 3 点にも依存する。
- Butler–Volmer の行:その場所の $j, \phi_s, \phi_e, c_{s,\mathrm{surf}}, c_e$ の 5 変数のみに非ゼロ。
非ゼロ要素の割合(密度)は数%以下で、$N$ を増やすほど下がる。 この疎性(sparsity)を活かす解法 — 帯行列用の LU($\mathcal{O}(N w^2)$、 $w$:帯幅 = 対角から最も遠い非ゼロまでの距離)や汎用の疎 LU — を使えば、 計算量はほぼ $N$ に比例するまで落ちる。 ただし帯幅 $w$ は変数の並べ方で激変する。図4.4 で確かめよう。
全 96 変数の小さな DFN のヤコビアン非ゼロパターン(色は方程式の由来)。 「フィールド順」は $[c_s \dots, c_e \dots, \phi_s \dots, \phi_e \dots, j \dots]$ と 変数の種類ごとに並べた順序で、ブロック構造が読み取りやすい。 ボタンで「ノード順」($x$ ノードごとに全変数をまとめる順序)に並べ替えると、 同じ行列なのに非ゼロが対角付近に集まり、帯幅が 82 → 11 に縮む。 帯 LU の計算量は帯幅の 2 乗に比例するので、並べ替えだけで約 50 倍速くなる勘定である。
scipy でも solve_ivp(..., method="BDF", jac_sparsity=S) と
疎パターン $S$ を渡すだけで、ヤコビアンの有限差分近似(色付き差分)と疎 LU が有効になり、
大規模問題が桁違いに速くなる。フェーズ5の自作ソルバーでは、
ノード順に並べた帯行列に対して scipy.linalg.solve_banded を使う。
4.6.6 コードで確かめる:同じ DAE を戦略Aと戦略Bで解く
仕上げに、4.5.3 節のミニ DAE($y' = -(z+z^3)$、$0 = z + z^3 - y$、厳密解 $y = e^{-t}$)を 戦略A(入れ子)と戦略B(一体)の両方で解いて、 同じ解に到達することを確認する。整合初期化($z(0)$ を Newton で求める)も入っている。 DAE ソルバーは魔法ではない — この 60 行が本質のすべてである。
# コード2: 半陽的 index-1 DAE を 2 通りの戦略で解いて一致を確認する
# 微分方程式: dy/dt = -(z + z^3) … 微分変数 y
# 代数拘束 : 0 = z + z^3 - y … 代数変数 z
# 拘束上では z + z^3 = y なので dy/dt = -y、つまり厳密解は y = e^{-t}。
import numpy as np
def g(y, z): return z + z**3 - y # 代数拘束 g(y,z) = 0
def dgdz(z): return 1.0 + 3.0*z**2 # ∂g/∂z(index-1: 常に > 0)
def solve_z(y, z_guess, tol=1e-12):
"""ニュートン法で g(y,z)=0 を z について解く(戦略 A の内側ループ)"""
z = z_guess
for _ in range(20):
dz = -g(y, z)/dgdz(z)
z += dz
if abs(dz) < tol:
return z
raise RuntimeError("Newton が収束しない")
h, t_end = 0.01, 3.0
N = round(t_end/h)
t = np.linspace(0.0, t_end, N+1)
# --- 整合的な初期値: y(0)=1 を与えたら z(0) は勝手に選べない ---
y0 = 1.0
z0 = solve_z(y0, 0.5) # z0^3 + z0 = 1 → z0 ≈ 0.6823
print(f"整合初期値: z(0) = {z0:.6f}(g(y0,z0) = {g(y0,z0):.1e})\n")
# --- 戦略 A: 入れ子方式(前進オイラー + 各ステップで代数式を Newton)---
yA = np.empty(N+1); yA[0] = y0
z = z0
for k in range(N):
z = solve_z(yA[k], z) # まず代数変数を解く
yA[k+1] = yA[k] - h*(z + z**3) # それを使って ODE を 1 歩進める
# --- 戦略 B: 一体方式(後退オイラー、[y,z] をまとめて Newton)---
# F(u) = [ y_{k+1} - y_k + h (z+z^3) ; z + z^3 - y_{k+1} ] = 0
yB = np.empty(N+1); yB[0] = y0
u = np.array([y0, z0])
for k in range(N):
u_new = u.copy() # 前ステップの値を初期推定に
for _ in range(20):
y_, z_ = u_new
Fvec = np.array([y_ - yB[k] + h*(z_ + z_**3), g(y_, z_)])
J = np.array([[1.0, h*dgdz(z_)],
[-1.0, dgdz(z_)]])
du = np.linalg.solve(J, -Fvec)
u_new += du
if np.max(np.abs(du)) < 1e-12:
break
u = u_new
yB[k+1] = u[0]
# --- 比較 ---
y_exact = np.exp(-t)
print(f"max|y_A - y_B| = {np.max(np.abs(yA - yB)):.3e}(差は O(h)。h→0 で一致)")
print(f"max|y_A - e^(-t)| = {np.max(np.abs(yA - y_exact)):.3e}(前進オイラー: 1 次精度)")
print(f"max|y_B - e^(-t)| = {np.max(np.abs(yB - y_exact)):.3e}(後退オイラー: 1 次精度)")
plot_spec = {
"traces": [
{"x": t.tolist(), "y": yA.tolist(), "name": "A: 入れ子(前進オイラー+Newton)", "mode": "lines"},
{"x": t.tolist(), "y": yB.tolist(), "name": "B: 一体(後退オイラー+Newton)",
"mode": "lines", "line": {"dash": "dash"}},
{"x": t.tolist(), "y": y_exact.tolist(), "name": "厳密解 e^(-t)",
"mode": "lines", "line": {"dash": "dot"}},
],
"layout": {"title": {"text": "index-1 DAE: 2 つの戦略と厳密解"},
"xaxis": {"title": {"text": "t [-]"}},
"yaxis": {"title": {"text": "y [-]"}}},
}2 つの戦略の差は $\mathcal{O}(h)$($h$:時間刻み)で、$h \to 0$ で同じ厳密解 $e^{-t}$ に収束する (前進・後退オイラーの誤差は符号が逆なので、差は各誤差の和程度になる)。 $h$ を半分にして再実行し、誤差が半分になること(1 次精度)も確かめてみてほしい。
4.7 検証の作法 — 収束次数と保存量
数値解は必ず「もっともらしい嘘」をつく。バグがあっても、それらしい曲線は出てくるからである。 だから数値計算には検証(verification)の作法がある。 本教材では次の 3 本柱を使う(フェーズ5で DFN ソルバーを組んだら、この 3 つを必ず実行する)。
4.7.1 解析解・参照解との比較
解析解が知られている特殊ケースに問題を退化させて比べるのが第一歩である。 球拡散なら、一定流束での長時間漸近解(放物線分布 + 線形降下)や、 フェーズ2で使った表面濃度の級数解が使える。 解析解がない一般の場合は、同じスキームの十分細かいグリッドの解を 「参照解(reference solution)」として使う。参照解は比較対象より 1 桁以上細かくし、時間刻みなど他の誤差要因は共通に揃えるのが作法である (時間誤差を測りたいときは逆に空間を固定して $\Delta t$ を振る)。
4.7.2 グリッド収束 — 次数を「測定」する
4.1〜4.2 節の理論によれば、本章の FVM は空間 2 次精度のはずである。 ならば、グリッド幅 $\Delta r$ を系統的に変えて誤差 $e(\Delta r)$ を測り、 両対数プロットの傾きから観測次数(observed order)
$$ p = \frac{\ln\big(e(\Delta r_1)/e(\Delta r_2)\big)}{\ln\big(\Delta r_1/\Delta r_2\big)} $$を計算して、理論値 2 と一致するか確かめられる。 次数まで一致して初めて「実装が正しい」と言える — 係数を 1 箇所書き間違えただけでも、収束はするが次数が 1 に落ちる、 といった形でバグが露見するからである。これは数値計算で最も強力なバグ検出器の一つである。
一定流束の球拡散を $N_r = 5$〜$80$ の 8 種類のグリッドで解き(時間積分は Crank–Nicolson、$\Delta t = 1$ s で共通)、参照解($N_r = 400$、同じ CN 法)との 相対 L2 誤差を両対数プロットした。破線は傾き 2($e \propto \Delta r^2$)の参照線。 データ点がこの傾きに乗っていること、最小二乗フィットの観測次数が $p \approx 2.0$ になることを確認してほしい。スライダーでどの評価時刻でも 次数が保たれることを見られる(全ケースは初期化時に計算済み)。
4.7.3 保存量のチェック
4.2 節で示したとおり、FVM では総 Li 量の収支が厳密に成り立つ:
$$ \frac{\mathrm{d}}{\mathrm{d}t}\left[\sum_i V_i c_i\right]_{\text{粒子}} = -\,4\pi R_s^2\,\frac{j}{F}, \qquad \frac{\mathrm{d}}{\mathrm{d}t}\big[\text{セル全体の総 Li 量}\big] = 0 . $$シミュレーション中にこの量を監視し、初期値との相対差が丸め誤差レベル ($10^{-12}$ 以下が目安。本章の実装では $10^{-15}$ 程度)に留まることを確認する。 保存量のドリフトは、境界条件の符号ミス・面フラックスの取り違え・ ソース項の付け忘れを高感度に検出する。 逆に、時間積分の許容誤差を緩めすぎたときに最初に現れる異常もこれである。
検証プロトコル(フェーズ5で実施する)
- 退化テスト:DFN のパラメータを退化させて SPM/SPMe(フェーズ2・3)を再現する。
- グリッド収束:$x$・$r$ 方向それぞれで観測次数 ≈ 2 を確認(図4.2 の手順)。
- 保存量:総 Li 量と電荷収支($i_s + i_e = I$)の恒等性を監視する。
- 外部参照:PyBaMM(Chen2020 パラメータ、戦略Cの IDAKLU)による 放電曲線・内部分布と比較する。プロ品質のソルバーとの突き合わせが最終試験である。
4.8 章末 — 理解度チェックと次章への橋渡し
理解度チェック
-
中心差分 $(c_{i+1} - c_{i-1})/(2\Delta x)$ が前進差分より 1 次だけ高精度になる理由を、
テイラー展開の言葉で説明せよ。
解答
$c_{i\pm1}$ のテイラー展開は $\Delta x$ の奇数次項の符号だけが異なる。 引き算すると偶数次項($c_i$ と $\tfrac{\Delta x^2}{2}c''$)が厳密に相殺し、 誤差の主要項が前進差分の $\tfrac{\Delta x}{2}c''$($1$ 次)から $\tfrac{\Delta x^2}{6}c'''$($2$ 次)に繰り上がるため。 対称性によるキャンセルが「タダで」1 次ぶんの精度を稼いでいる。 -
FVM の半離散式を全セルで足し合わせると総量の変化が境界フラックスだけで書ける理由と、
これが FDM では一般に保証されない理由を述べよ。
解答
FVM では内部の面 $r_{i+1/2}$ のモル流量 $W_{i+1/2}$ が、セル $i$ の式に $-W_{i+1/2}$、 セル $i+1$ の式に $+W_{i+1/2}$ として同じ値で現れる。 $V_i$ 倍して足すと望遠鏡和になり内部項がすべて相殺、 残るのは両端の境界フラックスのみ($\mathrm{d}/\mathrm{d}t \sum V_i c_i = W_{\text{内端}} - W_{\text{外端}}$)。 FDM は点での微分近似を各点で独立に行うため、隣接する式の打ち切り誤差が 相殺する保証がなく、足し合わせに $\mathcal{O}(\Delta x)$ 程度の残差が残り得る (係数が空間変化する場合や不等間隔グリッドで顕著)。 -
$D_s = 4 \times 10^{-15}$ m²/s(正極粒子)、$\Delta r = 1 \times 10^{-7}$ m のとき、
FTCS の安定限界 $\Delta t$ を求めよ。また、グリッドを 2 倍細かくすると限界は何倍になるか。
解答
$\Delta t \le \dfrac{\Delta r^2}{2 D_s} = \dfrac{(10^{-7})^2}{2 \times 4\times10^{-15}} = \dfrac{10^{-14}}{8\times10^{-15}} = 1.25$ s。 $\Delta r \to \Delta r/2$ とすると限界は $\Delta r^2$ に比例するので $1/4$(約 0.31 s)になる。1 時間の放電に必要な最少ステップ数は 2880 → 11520 と 4 倍に増える。 -
DFN の未知数 $c_s, c_e, \phi_s, \phi_e, j$ を微分変数と代数変数に分類し、
この DAE が index-1 である理由を一言で述べよ。
解答
微分変数:$c_s, c_e$(物質保存に $\partial/\partial t$ がある)。 代数変数:$\phi_s, \phi_e, j$(電荷保存は楕円型、Butler–Volmer は代数式)。 index-1 である理由:$c_s, c_e$ を固定したときの代数系 $\boldsymbol{g}(\boldsymbol{y}_d, \boldsymbol{y}_a) = \boldsymbol{0}$ のヤコビアン $\partial\boldsymbol{g}/\partial\boldsymbol{y}_a$ が(電位の基準を 1 点固定すれば)正則 — つまり拘束を 1 回時間微分するだけで $\mathrm{d}\boldsymbol{y}_a/\mathrm{d}t$ について 解けて ODE 系に直せるから。 -
ニュートン法が 2 次収束しているとき、ある反復で残差(誤差)が $10^{-2}$、
次の反復で $10^{-4}$ だった。さらに 2 反復後の誤差はおよそいくらか。
また、実際の計算でその値に到達できない理由は何か。
解答
2 次収束では $e_{k+1} \approx C e_k^2$($C \approx 1$ とみなす)なので $10^{-4} \to 10^{-8} \to 10^{-16}$。およそ $10^{-16}$。 ただし倍精度浮動小数点の相対精度(マシンイプシロン $\approx 2.2\times10^{-16}$)と 丸め誤差のため、実際には $10^{-15}$ 前後で頭打ちになる (図4.3 の表でも最後の反復で誤差が「ほぼ 0」に張り付くのが見える)。
次章への橋渡し:道具は揃った
本章で手に入れた道具を並べてみる:
- FVM(4.2)— $c_s$(球)、$c_e, \phi_s, \phi_e$($x$ 方向)をすべて保存形で離散化する。
- MOL + 陰解法(4.3–4.4)— 硬い系でも大きな $\Delta t$ で安定に進める。
- DAE の構造(4.5)— $\mathbf{M}\,\mathrm{d}\boldsymbol{y}/\mathrm{d}t = \boldsymbol{f}$、 index-1、整合初期化。
- 後退オイラー + ニュートン法 + 疎ヤコビアン(4.6)— 戦略Bの実装部品一式。
- 検証の作法(4.7)— 収束次数・保存量・参照解比較。
フェーズ5では、これらを全部組み立てて DFN を解く。 変数ベクトルの設計(図4.4 のノード順)、残差関数 $\boldsymbol{F}(\boldsymbol{y}^{k+1})$ の実装、 整合初期化、時間ループ、そして PyBaMM との突き合わせ検証 — 本章の各節が そのままフェーズ5の実装手順書になっていることに気づくはずである。