Nonlinear Solver Lab

W4法の使い方

  • このマニュアルでは簡単に連立非線形方程式ソルバー、 W4法の説明をします。 さらに詳しい情報が必要であれば、 参考文献 をご覧ください。

W4 利用ガイド

  1. 解きたい連立非線形方程式
  2. ニュートン・ラフソン法
  3. W4法
  4. W4の使い方

1. 解きたい連立非線形方程式

一般的に連立非線形方程式は、$\textbf{F}(\textbf{x})=\textbf{0}$のように表せる。
例えば、方程式の次元が３の場合、 $\textbf{F}$ は以下の３つの成分がある。

$f_{1}(x_1,x_2,x_3)=0$,
$f_{2}(x_1,x_2,x_3)=0$,
$f_{3}(x_1,x_2,x_3)=0$.

2. ニュートン・ラフソン法

$\textbf{F}(\textbf{x})=\textbf{0}$ という連立非線形方程式を解くためのニュートン・ラフソン法では、まず解 $\textbf{x}^{*}=\textbf{x}+\Delta\textbf{x}$の周りで方程式を線形化する。
$\textbf{F}(\textbf{x}^{*}) \sim \textbf{F}(\textbf{x}) +J\Delta\textbf{x}$
というテーラー展開において、 $J = \frac{\partial \textbf{F}}{\partial\textbf{x}}$ はヤコビアン行列である。

解の定義により、$\textbf{F}(\textbf{x}^{*})=\textbf{0}$であり、 次ステップの数値解は以下で定義できる。
$\textbf{x}^{n+1} = \textbf{x}^{n} -J^{-1}\textbf{F}(\textbf{x}^n)$.

ニュートン・ラフソン法は、連続極限で以下の微分方程式と見做すことができる。
$\dot{\textbf{x}} = -J^{-1}\textbf{F}(\textbf{x})$

3. W4法

W4法は２階微分項をニュートン・ラフソン法に付け加えたものとして以下のように定義できる。
$\ddot{\textbf{x}} +2\dot{\textbf{x}} -\dot{X}X^{-1}\dot{\textbf{x}} +XY\textbf{F}=\textbf{0}$.

自由度を一つ付け加え、 $X$ と $Y$という行列を適切に選ぶことで、より安全に連立非線形方程式を解けるようにしたものがW4法である。 通常は行列として$XY=J^{-1}$のように選ぶことで局所収束性を担保できるが、$X$の選び方には任意性があり、それによってW4の種類がいくつか存在する。

２階微分方程式は２つの１階微分方程式に分けられるが、具体的なW4法の方程式は以下で与えられる。
$\textbf{x}^{n+1} = \textbf{x}^{n} +\Delta\tau X\textbf{p}$
$\textbf{p}^{n+1} = (1-2\Delta\tau)\textbf{p}^{n} -\Delta\tau Y\textbf{F}$.

4. W4法の使い方

W4法の呼び出し関数には、どの種類のW4でも、８つの引数がある。 (v,F,Fa,J,ini,dt,itermax,errmax).

  • v :: 変数 (n次元 numpy array)
  • F :: n個の非線形方程式
  • Fa :: 非線形方程式の大きさ (収束判定用)
  • J :: Fから定義されるヤコビアン行列
  • ini:: vの初期条件
  • dt :: 時間ステップサイズ
  • itermax:: 最大イテレーション数
  • errmax :: 収束判定のためのエラー最大値

(1) 連立非線形方程式を定義する。

まず 方程式の次元を決める。dim
Sympy 形式で非線形方程式を定義する。
変数は x[0] から x[dim-1] で定義されている。

f = [];
f.append( 10*(x[1]-x[0]**2) );
f.append( 1 - x[0] );

(2) Numpy変数をSympy変数から定義する。

fから自動的に以下の関数でNumpy変数に変換できる。

vn, Fn, Fan, Jacn = calc_source(x,f);

(3) 変数 v の初期条件を決める。

Numpy形式で初期条件を決める。

vini = np.array([]);
vini = np.append(vini, 1.2);
vini = np.append(vini, 1.0);

(4) W4関数を呼び出す。

W4関数を呼び出して初期条件viniからイテレーションを始める。
もし、errやiterが errmaxやitermax 以下であれば問題が解けている。
最終的な答えは vans に格納されている。

vans, iter, err = w4ul(vn,Fn,Fan,Jacn,vini,dtau,itermax,errmax);
In [ ]: