数値連続法と分岐

FEMモデルの離散化から生じる代数的問題を,以下の形式で書くことができます.

\[F(U) = 0.\]

以下の説明では,モデルが追加のスカラーパラメータ \(\lambda\) に依存していると仮定します.そのため \(F(U) = F(U, \lambda)\) です.

数値連続法

数値連続法は,システムの解を追跡するために役立ちます

\[F(U, \lambda) = 0, \quad F\colon \mathbb{R}^{N} \times \mathbb{R} \to \mathbb{R}^{N}.\]

GetFEM には,区分的 \(C^{1}\) (\(PC^{1}\)) の解曲線の連続手法が実装されています(詳細は [Li-Re2014] を参照).状態変数 \(U\) とパラメータ \(\lambda\) との間に陽な違いがないので,簡潔にするため, \(Y := (U, \lambda)\) と表記します.併せて,例えば,接線を計算する際のスケーリングが悪くなるのを避けるために,以下の加重スカラ積とノルムを使用します.

\[\langle Y, \tilde{Y} \rangle_{w} := \kappa \langle U, \tilde{U} \rangle + \lambda \tilde{\lambda},\quad \lVert Y \rVert_{w} := \sqrt{\kappa \lVert U \rVert^{2} + \lambda^{2}},\qquad Y = (U, \lambda),\, \tilde{Y} = (\tilde{U}, \tilde{\lambda}).\]

ここで, \(\kappa\)\(\kappa \langle U, \tilde{U} \rangle\) が対応する空間変数のスカラー積に比例するように選ぶ必要があり,通常 \(L^{2}\) です.たとえば, \(\kappa = h^{d}\) 次のようにすることができます,ここで, \(h\) はメッシュサイズで, \(d\) は対称の問題の次元を表します.あるいは, \(\kappa\) は,簡単にするために \(1/N\) として設定することもできます.

連続法の考え方は,古典的な予測子修正子法によって滑らかな解曲線を追跡し,滑らかな部分を連続的に結合することです.

使用される特定の予測子修正子法は,MATCONT [Dh-Go-Ku2003] に実装されている 不正確なMoore-Penrose 連続をわずかに変更して使用します.これは,解曲線上におおよそ対応する単位接線ベクトル \(T_{j}\) のシーケンスである連続する点 \(Y_{j}\) の連続を計算します.

\[\lVert F(Y_{j}) \rVert \leq \varepsilon,\quad F'(Y_{j}; T_{j}) = 0,\quad \lVert T_{j} \rVert_{w} = 1,\quad j = 0, 1,\dotsc.\]

それを説明するために,上記の処理の関係を満たすカップル \((Y_{j}, T_{j})\) を考えます. 予測 では,次のような \((Y_{j+1}, T_{j+1})\) の初期近似が成り立ちます.

\[Y_{j+1}^{0} := Y_{j} + h_{j} T_{j},\quad T_{j+1}^{0} := T_{j},\]

ここで, \(h_{j}\) はステップサイズです.その選択については後で説明します.

補正 では,次のようなシーケンスが計算されます \(\{(Y_{j+1}^{l}, T_{j+1}^{l})\}\), ここで \(T_{j+1}^{l} := \tilde{T}_{j+1}^{l} / \lVert \tilde{T}_{j+1}^{l} \rVert_{w}\) とカップル \((Y_{j+1}^{l}, \tilde{T}_{j+1}^{l})\) は方程式 \(F^{l}(Y, T) = 0\) に適用されるNewton法のイテレーションにより与えられます.

\[\begin{split}F^{l}(Y, T) := \begin{pmatrix}F(Y)\\ (T_{j+1}^{l-1})^{\top}(Y - Y_{j+1}^{l-1})\\ \nabla F(Y_{j+1}^{l-1})T\\ \langle T_{j+1}^{l-1}, T \rangle_{w} - \langle T_{j+1}^{l-1}, T_{j+1}^{l-1} \rangle_{w}\end{pmatrix}\end{split}\]

初期近似 \((Y_{j+1}^{l-1}, T_{j+1}^{l-1})\) は次のようになります.\(F\) のポテンシャルの非微分可能性により,Newton法の微分的に滑らかな変形が使用されます( [Fa-Pa2003] のアルゴリズム7.2.14).

../_images/getfemusercorrection.png

補正.

カップル \((Y_{j+1}^{l}, T_{j+1}^{l})\)\((Y_{j+1}, T_{j+1})\) を満たす,もし \(\lVert F(Y_{j+1}^{l})\rVert \leq \varepsilon\) \(\lVert Y_{j+1}^{l} - Y_{j+1}^{l-1}\rVert_{w} \leq \varepsilon'\), と \(T_{j+1}^{l}\)\(T_{j}\) の間の角度の余弦は \(c_{\mathrm{min}}\) のものより大きいか等しい(または,非微分可能性の場合にはその選択関数の1つ)の部分勾配 \(F\) が,解析的に構築られているのに対して, \(\ lambda\) は1e-8に等しい増分で前方有限差分法によって評価されます.

次の予測におけるステップサイズ \(h_{j+1}\) は,Newton補正がどのように成功したかに依存します.次の式で必要な反復回数 \(l_{\mathrm{it}}\) を示します.

\[\begin{split}h_{j+1} := \begin{cases}\max\{h_{\mathrm{dec}} h_{j}, h_{\mathrm{min}}\}& \text{if no new couple has been accepted},\\ \min\{h_{\mathrm{inc}} h_{j}, h_{\mathrm{max}}\}& \text{if a new couple has been accepted and } l_{\mathrm{it}} < l_{\mathrm{thr}},\\ h_{j}& \text{otherwise},\end{cases}\end{split}\]

ここで, \(0 < h_{\mathrm{dec}} < 1 < h_{\mathrm{inc}}\), \(0 < l_{\mathrm{thr}}\)\(0 < h_{\mathrm{min}} < h_{\mathrm{max}}\) は定数となります.最初は,いくつかの \(h_{\mathrm{min}} \leq h_{\mathrm{init}} \leq h_{\mathrm{max}}\) に対して \(h_{1} := h_{\mathrm{init}}\) のように設定します.

ここで,滑らかな振る舞いの1つのサブ領域に対応する解曲線の一部を近似したとします. \(F\) の滑らかな振る舞いの別のサブ領域に対応する部分を回復したいと考えます.\((Y_{j},T_{j})\) は最後に計算されたカップルとします.

../_images/getfemusertransition.png

解曲線の滑らかな区分間の遷移.

他の滑らかな区分の接線を近似するには,まず,次のような点 \(Y_{j} + h T_{j}\) をとります.ここで, \(h\)\(h_{\mathrm{min}}\) よりも少し大きいです,そのためこの点が円滑な動作の他のサブ領域の内部に属するようにします.それから, \(\tilde{T}\) は以下となります

\[\nabla F(Y_{j} + h T_{j}) \tilde{T} = 0,\quad \lVert \tilde{T} \rVert_{w} = 1,\]

そしてこれはこのベクトルの適切な方向を決定するために残っています.これは,以下の観察に基づいて行うことができます.第1に, \(Y_{j} - r \tilde{h} \tilde{T}\)\(Y_{j}\) の任意の正の \(\tilde{h}\) と同じサブ領域内にあるような \(r \in \{\pm 1\}\) が存在します.これは, \(\nabla F(Y_{j} - r \tilde{h} \tilde{T}) T_{-} = 0\)\(\frac{\lvert T_{-}^{\top} \tilde{T}\rvert}{\lVert T_{-} \rVert \lVert \tilde{T} \rVert}\) が1よりもかなり小さいという事実によって特徴づけられます.第2に, \(Y_{j} + r \tilde{h} \tilde{T}\) は,ある正の閾値よりも大きい \(\tilde{h}\) のために,他のサブ領域に現れ,このような値の場合, \(\frac{\lvert T_{+}^{\top} \tilde{T}\rvert}{\lVert T_{+} \rVert \lVert \tilde{T} \rVert}\)\(T_{+}\)\(\nabla F(Y_{j} + r \tilde{h} \tilde{T}) T_{+} = 0\) の場合に1に近いです.

これは, \(\tilde{T}\): の所望の方向を選択するための以下の手順を提案します. \(\tilde{h}\) の値を \(h_{\mathrm{min}}\) から連続的に増加させ, \(\tilde{h}\)\(r \in \{\pm 1\}\) に到着したときに

\[\frac{\lvert T^{\top} \tilde{T}\rvert}{\lVert T \rVert \lVert \tilde{T} \rVert} \approx 1\quad \text{if}\ \nabla F(Y_{j} + r \tilde{h} \tilde{T}) T = 0,\]

他の滑らかな区分の接線の近似値として \(r \tilde{T}\) を取ります.

この近似を使って,予測子修正子を \((Y_{j}, r \tilde{T})\) のように再始動します.

GetFEM では,モデルパラメータ化の2つの連続の方法が実装されています.

  1. パラメータ \(\lambda\) は,モデルが依存するスカラーデータです.

  2. このモデルはモデルが依存するベクトルデータ \(P\) によって スカラーパラメータ \(\lambda\) によりパラメータ化されます.この場合,線形経路を使用します

    \[\lambda \mapsto P(\lambda) := (1 - \lambda)P^{0} + \lambda P^{1},\]

    ここで, \(P^{0}\)\(P^{1}\) には, \(P\) の値が与えられ,問題の解の集合をトレースします.

    \[F(U, P(\lambda)) = 0.\]

限界点の検出

システム \(F(U,\lambda) = 0\) の解を追跡する場合, 限界点 (折り返し点または折り返し点とも呼ばれます)に注目します,これにより同じ値に対する解の数 \(\lambda\) は変わります.これらの点は,試験関数 \(\tau_{\mathrm{LP}}\) の符号変化によって検出することができます.

\[\tau_{\mathrm{LP}}(T_{j}) \tau_{\mathrm{LP}}(T_{j+1}) < 0,\]

ここで, \(\tau_{\mathrm{LP}}\) は次のように定義されます.

\[\tau_{\mathrm{LP}}(T) := T_{\lambda},\quad T = (T_{U},T_{\lambda}) \in \mathbb{R}^{N} \times \mathbb{R}.\]
../_images/getfemuserlimitpoint.png

限界点.

数値分岐

\(\bar{Y}\) はシステム \(F(Y) = 0\)分岐点 と呼ばれます \(F(\bar{Y}) = 0\) ならば, 2つ以上の異なる解曲線がそれを通過します.次の結果は, 滑らかな 分岐点のテストを与えます(例えば, [Georg2001] を参照).

Let \(s \mapsto Y(s)\) be a parameterization of a solution curve and \(\bar{Y} := Y(\bar{s})\) be a bifurcation point. Moreover, let \(T^{\top} \dot{Y}(\bar{s}) > 0\), \(B \notin \mathrm{Im}(J(\bar{Y}))\), \(C \notin \mathrm{Im}(J(\bar{Y})^{\top})\), \(d \in \mathbb{R}\) and

\[\begin{split}J(Y) := \begin{pmatrix}\nabla F(Y)\\ T^{\top}\end{pmatrix}.\end{split}\]

以下で \(\tau_{\mathrm{BP}}(Y)\) を定義します.

\[\begin{split}\begin{pmatrix}J(Y)& B\\ C^{\top}& d\end{pmatrix} \begin{pmatrix}V(Y)\\ \tau_{\mathrm{BP}}(Y)\end{pmatrix} = \begin{pmatrix}0\\ 1\end{pmatrix}.\end{split}\]

次に \(\tau_{\mathrm{BP}}(Y(s))\) は, \(s = \bar{s}\) で記号を変換します.

自明ですが, \(B\)\(C\)\(d\) をランダムに取る場合,それらが上記の要件を満たす可能性が高くなります.したがって,数値連続法は,各連続ステップでの補正によって提供されたベクトル \(Y\)\(T\) を取り,次の \(\tau_{\mathrm{BP}}\) の符号を監視することによって分岐点を検出することができます.

分岐点 \(\bar{Y}\)\(\tau_{\mathrm{BP}}(Y_{j}) \tau_{\mathrm{BP}}(Y_{j+1}) < 0\) の符号の変化によって検出されると,これは,特別なステップ長適応([Al-Ge1997] の8.1節を参照)を用いて上記の予測子修正子のステップによってより正確に近似することができます.すなわち,次のステップ長を

\[h_{j+1} := -\frac{\tau_{\mathrm{BP}}(Y_{j+1})}{\tau_{\mathrm{BP}}(Y_{j+1}) - \tau_{\mathrm{BP}}(Y_{j})}h_{j}\]

\(\lvert h_{j+1} \rvert < h_{\mathrm{min}}\) のようになるまで,関数 \(s \mapsto \tau_{\mathrm{BP}}(Y(s))\) のゼロを見つけるための割線のメソッドに対応させます.

最後に,解の分岐を切り替えたいとしましょう.この目的のために,我々は,2つの異なる解曲線が交差する,いわゆる 単純な分岐点 の場合を考えます.

\(\tilde{Y}\) は,与えられた \(\bar{Y}\) の近似であり, \(V(\tilde{Y})\) は,試験関数 \(\tau_{\mathrm{BP}}(\tilde{Y})\) を計算するための拡張システムの解の最初の部分です. [Georg2001] で提案されているように, \(V(\tilde{Y})\) を予測変数として使用し, \((\tilde{Y}, V(\tilde{Y}))\) から始まる1つの連続ステップを実行して,新しい分岐上の点を取得できます.この連続ステップが正常に実行され,新しい分岐上の点がリカバリされた後,通常の予測子修正子ステップを続行して,この分岐をトレースすることができます.

最近では, \(PC^{1}\) -分岐のための数値計算ツールが GetFEM で開発されました.\(J\) は,次の式で定義される実数パラメータの行列関数となります.

\[\begin{split}J(\alpha) := (1-\alpha)\begin{pmatrix}\nabla F(Y_{j})\\ T_{j}^{\top}\end{pmatrix} + \alpha\begin{pmatrix}\nabla F(Y_{j+1})\\ T_{j+1}^{\top}\end{pmatrix}.\end{split}\]

[Li-Re2014hal] で提案されているように,以下のテストは, \(Y_{j}\)\(Y_{j+1}\) の間の分岐点 \(PC^{1}\) の検出に使用できます.

\[\det J(0) \det J(1) < 0.\]

このテストを数値的に実行するために,以下を導入します.

\[\begin{split}M(\alpha) := \begin{pmatrix}J(\alpha)& B\\ C^{\top}& d\end{pmatrix}\end{split}\]

\(\tau_{\mathrm{BP}}(\alpha)\) で上記と同様に

\[\begin{split}M(\alpha) \begin{pmatrix}V(\alpha)\\ \tau_{\mathrm{BP}}(\alpha)\end{pmatrix} = \begin{pmatrix}0\\ 1\end{pmatrix}.\end{split}\]

これは,Cramer法から,次の通りとなります.

\[\tau_{\mathrm{BP}}(\alpha) = \frac{\det J(\alpha)}{\det M(\alpha)}\]

ただし, \(\det M(\alpha)\) はゼロではありません.したがって, \(\det J(\alpha)\) がゼロのときに \(\det M(\alpha)\) が0でないように \(B\)\(C\) および \(d\) が選択された場合, \(\det J(\alpha)\) の符号変更は \(\tau_{\mathrm{BP}}(\alpha)\) から0の通過によって特徴付けられ,\(\det M(\alpha)\) の符号は \(\tau_{\mathrm{BP}}(\alpha)\) の符号変化によります. \(\det M(\alpha)\) の符号変化が特異点に起因して変化することによります.結論として, \(\det J(0)\det J(1)\) の符号は \(\alpha\)\([0,1]\) を通過するときの \(\tau_{\mathrm{BP}}(\alpha)\) の挙動を追跡し, \(\det J(\alpha)\) の符号変化を監視することによって決定されます.

[Li-Re2014hal] で示されているように, \(B\)\(C\)\(d\) はランダムに選択することができます.現在の値 \(\alpha\) のインクリメント \(\delta\) は適宜変更され,\(\tau_{\mathrm{BP}}\) の特異点が有効に扱われます. \(\tau_{\mathrm{BP}}(\alpha)\) のそれぞれの計算の後, \(\delta\) は次のように設定されます.

\[\begin{split}\delta := \begin{cases}\min\{2\delta, \delta_{\mathrm{max}}\}& \text{if $\lvert \tau_{\mathrm{BP}}(\alpha) - \tau_{\mathrm{BP}}(\alpha - \delta) \rvert < 0.5 \tau_{\mathrm{fac}} \tau_{\mathrm{ref}}$,}\\ \max\{0.1\delta, \delta_{\mathrm{min}}\}& \text{if $\lvert \tau_{\mathrm{BP}}(\alpha) - \tau_{\mathrm{BP}}(\alpha - \delta) \rvert > \tau_{\mathrm{fac}} \tau_{\mathrm{ref}}$,}\\ \delta& \text{otherwise},\end{cases}\end{split}\]

ここで, \(\delta_{\mathrm{max}} > \delta_{\mathrm{min}} > 0\) で, \(\tau_{\mathrm{fac}} > 0\) は既定の定数で, \(\tau_{\mathrm{ref}} := \max\{\lvert \tau_{\mathrm{BP}}(1) - \tau_{\mathrm{BP}}(0) \rvert, 10^{-8}\}\) です.

\(Y_{j}\)\(Y_{j+1}\) の間で分岐点 \(PC^{1}\) が検出されると,それは二分法のような手順でより正確に近似されます.得られた近似は, \(Y_{j},\) と同じ滑らかな枝上にあり,滑らかさの対応する領域から対応する単位接線も計算されます.

連続な場合とは対照的に, \(PC^{1}\) 分岐点からどれくらいの分岐があるのか,またそれらが求められる方向は明確ではありません.このため,新しい分岐上の点を見つけるために,予測方向のシーケンス全体の連続ステップが試されます.

分岐点の近似値を表す \(\tilde{Y}\), \(\tilde{T}\) と対応する接線はそれぞれ次のようになります.いくつかの参照ベクトル \(\tilde{V}_{1}\)\(\tilde{V}_{2}\)\(\pm V\) となり \(V\) は次を満たします.

\[\nabla F(\tilde{Y}+h_{\mathrm{min}}\tilde{V}) V = 0, \quad \lVert V \rVert_{w} = 1,\]

\(\tilde {V}\) は,\(\tilde{V}_{1}\)\(\tilde{V}_{2}\) のような線形結合で表されます.線形結合の総数は, \(n_{\mathrm{dir}},\) で与えられ,参照ベクトルは次の方針に従って連続して選択されます.

  1. 次のように計算します \(\tilde{V}_{1} := -\tilde{T}\)\(\tilde{V}_{2}\)

    \[\nabla F(\tilde{Y}+h_{\mathrm{min}}\tilde{T}) \tilde{V}_{2} = 0, \quad \lVert \tilde{V}_{2} \rVert_{w} = 1.\]
  2. 次の点に対応する単位接線の集合を表します. \(\{\tilde{T}_{1},\dotsc\tilde{T}_{n_{\mathrm{br}}}\}\) 分岐点から枝分かれする方向に向いている.それから \(\tilde{V}_{2}\) は, \(\{\tilde{T}_{1},\dotsc\tilde{T}_{n_{\mathrm{br}}}\}\) から異なる組み合わせとして逐次取られます.

  3. これまで利用可能だったすべての組み合わせがすでに使用されている場合, \(\tilde{V}_{1}\) は変更せずに \(\tilde{V}_{2} := \tilde{V}_{2}^{+}\)\(\tilde{V}_{2}^{+}\) が満たすようにしてください.

    \[\nabla F\Bigl(\tilde{Y}+h_{\mathrm{min}}\Bigl(\tilde{V}_{2}^{-} + 0.1\frac{\tilde{V}_{3}}{\lVert \tilde{V}_{3} \rVert_{w}}\Bigr)\Bigr) \tilde{V}_{2}^{+} = 0, \quad \lVert \tilde{V}_{2}^{+} \rVert_{w} = 1.\]

    ここで, \(\tilde{V}_{2}^{-}\) は以前に採用されたベクトル \(\tilde{V}_{2}\) と, \(\tilde{V}_{3}\) がランダムに選択されます.

\(\tilde{V}_{1}\)\(\tilde{V}_{2}\) の合計数は,式 \(n_{\mathrm{span}}\) で与えられます.

\(PC^1\) 数値分岐のより詳細な情報は, [Li-Re2015hal] にあります.

モデルの解曲線の近似

数値連続法は getfem/getfem_continuation.h で定義されています.それを使用するには,最初に対応するオブジェクトを介してそれを設定する必要があります

getfem::cont_struct_getfem_model S(model, parameter_name, sfac, ls, h_init, h_max, h_min, h_inc, h_dec,
                                   maxit, thrit, maxres, maxdiff, mincos, maxres_solve, noisy, singularities,
                                   non-smooth, delta_max, delta_min, thrvar, ndir, nspan);

ここで parameter_name\(\lambda\) を表すモデルデータの名前で, sfac はスケールファクタ \(\kappa\) を表し, ls はプロセスに組み込まれた線形システムに使用されるソルバの名前です(例えば,getfem::default_linear_solver<getfem::model_real_sparse_matrix, getfem::model_real_plain_vector>(model)).実数 h_inith_maxh_minh_inch_dec\(h_{\mathrm{init}}\)\(h_{\mathrm{max}}\)\(h_{\mathrm{min}}\)\(h_{\mathrm{inc}}\) ,および \(h_{\mathrm{dec}}\) を表し,整数 maxitthritmaxresmaxdiffmincosmaxres_solve\(l_{\mathrm{thr}}\)\(\varepsilon\)\(\varepsilon'\)\(c_{\mathrm{min}}\) であり,線形システムの目標残差値をそれぞれ求めるために設定されます.非負の整数 noisy は,連続プロセスの過程でどの程度詳細な情報を表示するかを決定します(詳細なほど値が大きくなります).整数 singularities は,特異点を指定する必要があります(完全に無視する場合は0,限界点を検出する場合は1,分岐点を検出して処理する場合は2).ブール値 nosmooth は,滑らかな連続法か分岐かを指定する必要があります.また,滑らかでないツールでも手法を指定する必要があります.実数 delta_maxdelta_min ,および thrvar は,\(\delta_{\mathrm{max}}\), \(\delta_{\mathrm{min}}\)\(\tau_{\mathrm{fac}}\), を表します. ndirnspan\(n_{\mathrm{dir}}\)\(n_{\mathrm{span}}\) をそれぞれ示します.

必要に応じて,ベクトルデータによるパラメータ化は,次のように宣言されます.

S.set_parametrised_data_names(initdata_name, finaldata_name, currentdata_name);

ここで,データ名 initdata_namefinaldata_name は,それぞれ, \(P^{0}\)\(P^{1}\) を表す必要があります. currentdata_name の下には, \(P(\lambda)\) の値が格納されていなければなりません.すなわち,モデルが依存するデータの実際の値です.

次に,連続法は次の式で初期化されます

S.init_Moore_Penrose_continuation(U, lambda, T_U, T_lambda, h);

ここで, U\(Y_{0}=\) (U,lambda) におけるパラメータの値の解でなければなりません.この初期化の間,初期値は T_lambda の符号に従って計算され, \(Y_{0}\) に対応する初期単位接線 \(T_{0}\) が返され, T_UT_lambda が返されます.さらに, h は初期ステップサイズ h_init に設定されます.

その後,連続の際の1ステップは次のように呼び出されます

S.Moore_Penrose_continuation(U, lambda, T_U, T_lambda, h, h0);

各呼び出しの後,解曲線上の新しい点および対応する接線が,変数 Ulambda および T_UT_lambda で返されます.次の予測のステップサイズは h で返されます.現在のステップのサイズはオプションの引数 h0 に返されます.選択された 特異点 の値によれば,各連続ステップの終わりに,限界点および分岐点の試験関数が評価されます.さらに,滑らかな分岐点が検出されれば,数値分岐の手順が実行され,分岐点の近似および両方の分岐曲線に対する接線が連続オブジェクト S に保存されます.そこから,次のいずれかのカーブを連続的にトレースするため初期化できるように,S のメンバ関数で簡単に回復することができます.

平滑問題を使用するテストプログラムの例全体は tests/test_continuation.cc, interface/tests/matlab/demo_continuation.minterface/src/scilab/demos/demo_continuation.sce にあります, 一方 interface/src/scilab/demos/demo_continuation_vee.sceinterface/src/scilab/demos/demo_continuation_block.sce は平滑でない例になります.