Python : 非線形（連立）方程式の解法

pythonで数値的に非線形連立方程式の解法のメモです。

問題

以下の非線形連立方程式を例に解き方を整理します。

$$
\left\{ \ \begin{aligned}
& 2x^2 - y^2 = 0 \\
& y = 0.5 \left(  {\rm sin}(x) + {\rm cos}(y) \right) 
\end{aligned} \right.
$$

以下の関数にして、関数のゼロ（root）を探すことで解を得えます。

$$
\left\{ \ \begin{aligned}
& f_0 \left(x,y \right) = 2x^2 - y^2  \\
& f_1 \left(x, y \right) = y - 0.5 \left(  {\rm sin} (x) + {\rm cos} (y) \right) 
\end{aligned} \right.
$$

解法

今回は、Scipy の optimize パッケージにあるrootで解いてみます。以下は、マニュアルです。（マニュアルの中に例題もあります。）

scipy.optimize.root — SciPy v1.10.0 Manual

methodは、hybr としておきます。
MINPACK に実装されているものを利用しているようです。
MINPACKの詳細は以下にのリファレンスがありました。
More, Jorge J., Burton S. Garbow, and Kenneth E. Hillstrom. 1980. User Guide for MINPACK-1.
検索すると、PDFがでてきます。ここまでは読んでいません。

 $${ x }$$ -> x[0],  $${ y }$$ ->x[1]です。$${ f_0(x,y) }$$ -> f[0], $${ f_1(x,y) }$$ -> f[1] と置いています。

import numpy as np
import scipy as sp
def res(x):
    f = np.zeros(2, dtype=np.float64)
    f[0] = 2*x[0]**2 - x[1]**2
    f[1] = x[1] - 0.5 * (np.sin(x[0]) + np.cos(x[1]))
    return f
x0 = np.array([0.0, 0.0])
ans = sp.optimize.root(res, x0, method='hybr')
print(ans)
print('--------')

以下が、出力結果

    fjac: array([[-0.95931118,  0.28235096],
       [-0.28235096, -0.95931118]])
     fun: array([-7.21644966e-15, -1.11022302e-16])
 message: 'The solution converged.'
    nfev: 18
     qtf: array([-1.41162713e-11, -4.36299709e-12])
       r: array([-1.6134865 ,  1.02076612, -1.052239  ])
  status: 1
 success: True
       x: array([0.43781417, 0.61916273])

一番最後の行が答えです。この値の時の関数の値は、
fun: array([-7.21644966e-15, -1.11022302e-16])　です。ほぼゼロで、収束しているようです。

例題：状態方程式を解く

連立ではないですが、化学工学で出てくる非線形方程式を例にといてみます。
Soave-Redlich-Kwong（SRK）equation of state を解いてある状態のモル体積を求めてみます。

SRK EOSは、

$$
P = \frac{RT}{\hat{V} - b} - \frac{\alpha a}{ \hat{V} (\hat{V} + b ) }   
$$

$$
a = 0.42747 \frac{(RT_c)^2}{P_c} \\
b = 0.08664 \frac{R T_c}{P_c} \\
m = 0.48505 + 1.55171 \omega - 0.1561 \omega^2 \\
\alpha = \left[ 1 + m \left( 1 - \sqrt{T / T_c} \right) \right]^2 
$$

と定義されています。
$${ P }$$ : 圧力[atm]、$${T}$$ : 絶対温度[K]、$${\hat{V}}$$ : モル体積[L/mol]、$${T_c}$$ : 臨界圧力[K]、$${P_c}$$ : 臨界圧力[atm]、$${\omega}$$ : 偏心因子[-]、$${R}$$ : 気体定数[L atm/(mol-K)]です。

例として、プロパン（臨界圧力 41.93437[atm]、臨界温度 369.82[K]、偏心因子 0.152[-]）の350K, 15 atm におけるモル体積を求めてみます。状態方程式を移行して、以下の関数がゼロになる$${ \hat{V} }$$を探します。

$$
f(\hat{V}) = P - \frac{RT}{\hat{V} - b} + \frac{\alpha a}{ \hat{V} (\hat{V} + b ) }
$$

import numpy as np
import scipy as sp
def f(V):
    # Properties of Propane
    T = 350 # K
    P = 15 # atm
    Tc = 369.82 # K
    Pc = 41.93437 # atm
    w = 0.152 # Acentric factor for CO
    R = 0.08206 # L atm / (mol K)
    a = 0.042747*(R*Tc)**2 / (Pc)
    b = 0.08664*(R*Tc/Pc)
    m = 0.48508 + 1.55171*w - 0.1561*w**2
    alpha = (1 + m*(1 - np.sqrt(T/Tc)))**2
    term1 = R*T/(V-b)
    term2 = alpha*a/(V*(V+b))
    return P - term1 + term2

V = 2.0 # L/mol, initial guess
V = sp.optimize.root(f, V) #, maxiter=100, f_tol=1e-6)
print(V)

結果は、
    fjac: array([[-1.]])
     fun: array([5.55111512e-17])
 message: 'The solution converged.'
    nfev: 7
     qtf: array([6.03536665e-12])
       r: array([-7.84436269])
  status: 1
 success: True
       x: array([1.94609414])

答えは、x: array([1.94609414])　です。fの値もほぼゼロなので問題なさそうです。
1.946 L/mol と求まりました。

まとめ

Scipy.optimize.rootを使って非線形方程式を解いてみました。
root関数の中身までは終えていません。

所感

MINPACKは、1980年にまとめられたプロジェクトのようです。研究者たちの蓄積のおかげで簡単に方程式が扱えるのが、有難です。