非線形方程式の解法について
代表的な非線形方程式の数値計算解法についてまとめていきます。
・2分法
微分代数方程式 (Differential-algebraic system of equations)
逐次代入法
逐次代入法は、方程式 \(f(x)=0\) を \(x=g(x)\) に変形し、\(x_{k+1}=g(x_k)\) を繰り返して解を探す方法です。
収束の目安は、解の近傍で \(|g'(x)|<1\) が成り立つことです。変形の仕方によって収束・発散が大きく変わるため、\(g(x)\) の選び方が重要です。
配管の摩擦係数を求める際に使えます。
こちらに使い方をまとめています。
Steffensen法
Steffensen法は、逐次代入法を差分で加速した手法です。導関数を使わずに2次収束に近い挙動を示すため、Newton法の導関数計算を避けたい場合に有効です。
更新式の一例は \[ x_{k+1}=x_k-\frac{(g(x_k)-x_k)^2}{g(g(x_k))-2g(x_k)+x_k} \] です（分母が極小になる場合は不安定化に注意）。
BWRS式の解を得る際に使用してみました。
こちらに計算結果をまとめています。
2分法
2分法は、符号が異なる2点 \([a,b]\)（\(f(a)f(b)<0\)）を初期区間に取り、中点 \(m=(a+b)/2\) の符号で区間を半分ずつ狭める方法です。
収束は遅めですが、連続関数で初期区間が正しく取れていれば確実に収束します。1回ごとに区間幅が1/2になるため、必要反復回数の見積りが容易です。
はさみうち法
はさみうち法（Regula Falsi）は、区間 \([a,b]\) の端点を結ぶ直線と x 軸の交点を次点として使う方法です。2分法より速くなることが多く、初期区間の保証も維持できます。
一方で、関数形によっては片側端点が更新され続けず、収束が遅くなることがあります（改良法として Illinois 法などがあります）。
ニュートン法
ニュートン法は \[ x_{k+1}=x_k-\frac{f(x_k)}{f'(x_k)} \] で更新する方法です。初期値が適切で \(f'(x)\neq0\) なら非常に高速（2次収束）です。
ただし、初期値が悪いと発散・振動することがあり、導関数の評価が必要です。実務では反復上限・ステップ制限を併用して安定化します。
各数値計算法の特徴
２分法、はさみうち法、ニュートン法のそれぞれの特徴を状態方程式の３次式を求める例で確認してみます。
以下の条件で、Methaneの圧縮係数をPeng-Robinson(PR)で求めてみます。
温度：111.587[K]、圧力：101.325[kPa]。なお、臨界定数は以下の値を使用しています。
偏心因子：0.01155[-]、臨界温度：190.56[K]、臨界圧力：4599[kPa]。
詳細な計算は、いずれVBAのページにまとめていきますので、ここでは簡単に結果だけまとめます。
物性推算法のPRのページで紹介したように、zの３次方程式を求めると、今回の条件では以下の式となります。
$$f(z)=z^3-0.997z^2+0.0291z-9.41\times10^{-5}$$z対f(z)をグラフで表すと、以下のようになります。
各解法での収束状況をグラフとテーブルで確認してみます。
まずはテーブル
グラフは以下のようになります。
今回収束の条件として誤差10^-8とし、ほぼ同様の初期条件から計算を開始しておりましたが、2分法、はさみうち法は収束まで20回以上の計算が必要になっています。
２分法には更に注意点があります。
微分代数方程式
微分代数方程式（DAE）は、微分方程式と代数方程式が連成した系です。化学プロセスの動特性モデルでは、物質収支・エネルギー収支（微分）と平衡式・状態方程式（代数）が同時に現れます。
初期値は微分変数だけでなく代数拘束を満たす必要があり（consistent initial condition）、通常の ODE より設定が難しくなります。解法としては後退オイラーやBDF系の陰解法がよく使われます。
当サイトに不具合、ご意見等ございましたらCEsolutionにお知らせください。