Chkder Hybrd Hybrd1 Hybrd1_r Hybrd_r Hybrj Hybrj1 Hybrj1_r Hybrj_r Sos Sos_r

Sos_r() Sub Sos_r ( N As  Long, X() As  Double, RTolx As  Double, ATolx As  Double, Tolf As  Double, Info As  Long, K As  Long, XX() As  Double, YY As  Double, IRev As  Long, Optional Info2 As  Long, Optional Nprint As  Long = 0, Optional MaxIter As  Long = 0, Optional Iter As  Long, Optional RNorm As  Double  ) 連立非線形方程式 (ブラウンの方法) (リバースコミュニケーション版) 目的 本ルーチンはn本のn変数非線形方程式よりなる連立方程式 fi(x1, x2, ..., xn) = 0 (i = 1 〜 n) の解をブラウンの方法により求める. アルゴリズムはニュートン法を変形した反復法に基づいており, ガウス−ザイデル法の手順に似たガウス消去法を用いる. 収束はおおむね二次収束である. 本ルーチンの収束特性は方程式の順番に影響を受ける. 線形および弱い非線形の方程式を先に並べるとよい. IRevに従って関数値をユーザーが与える必要がある. アルゴリズムで必要なすべての偏微分係数は差分商で近似される. 引数 [in] N 未知数および方程式の数. (N > 0) [in,out] X() 配列 X(LX - 1) (LX >= N) [in] 初期近似解. [out] 求められた解ベクトル. [in] RTolx 収束判定に用いられる相対誤差許容値(RTolx >= 0). 次式により解の精度をチェックする.   Abs(X(I) - Xold(I)) <= RTolx*Abs(X(I)) + ATolx ここで, Xold(I)は前回の反復時の値を表す. [in] ATolx 収束判定に用いられる絶対誤差許容値(ATolx >= 0). 解ベクトルの要素のいくつかが0の可能性があるならば, 計算効率をよくするためにATolxを適切な値(残りの要素のスケールに依存)に設定しなければならない. [in] Tolf 残差の収束判定基準(Tolf >= 0). 全ての式の残差(関数値)がこれより小さくなれば収束と見なす. この収束判定は方程式のスケーリングに依存するため, Tolfに適切な値を与えるように細心の注意が必要である. 値が適切でない場合, 不十分な反復で終了してしまうことがある. [out] Info = 0: 正常終了. (サブコードをInfo2に返す) = -1: パラメータ N の誤り. (N < 1) = -2: パラメータ X() の誤り. = -3: パラメータ RTolx の誤り. (RTolx < 0) = -4: パラメータ ATolx の誤り. (ATolx < 0) = -5: パラメータ Tolf の誤り. (Tolf < 0) = 1: 反復回数の上限に達した. (収束が非常に遅い可能性がある. 誤差許容値が適切であるか確認せよ) = 2: 反復回数の上限に達した. (局所的な最小点にあるか計算精度不足の可能性がある) = 3: 反復が発散した. (残差ノルムが数回の連続した反復において増加した) = 4: ヤコビ行列に関連する行列が特異になった. [out] K IRev = 1: K番目の関数値を求めることを示す. IRev = 3: ここまでの反復回数. [out] XX() 配列 XX(LXX - 1) (LXX >= N) IRev = 1: 関数値を求めるべき点を返す. IRev = 3: 最新の近似解. [in,out] YY [in] IRev = 1: 再呼び出し時に関数値fi(XX()) (i = K)を与えること. [out] IRev = 3: 最新の近似解における残差ノルム. [in,out] IRev リバースコミュニケーションの制御変数. [in] 最初の呼び出し時に 0 に設定しておくこと. 2回目以降の呼び出し時には値を変更してはならない. [out] 0 以外の時には下記処理を行ってから再び本ルーチンを呼び出すこと.   = 0: 処理終了. 正常終了かどうかはInfoをチェックすること.   = 1: XX()におけるK番目の関数値を求め YYに設定すること. YY以外の変数を変更してはならない.   = 3: 途中結果(K: 反復回数, YY: XXにおける残差ノルムなど)を出力する(Nprint = 1の場合). どの変数も変更してはならない. [out] Info2 (省略可) Info = 0 のときのサブコード. = 1: 解ベクトルの各要素X(i)が誤差許容値判定条件を満たした: abs(X(i) - Xold(i)) <= RTolx*abs(X(i)) + ATolx. = 2: すべての残差がTolf以下になった: abs(F(X, i)) <= Tolf. = 3: 上の2つ共満たした. = 4: 収束と思われる(誤差許容値が小さすぎる, あるいは, 収束の遅い関数である). (残差ノルムが数回の連続した反復においてほぼ一定であった) [in] Nprint (省略可) 中間結果出力の制御. (省略時 = 0) = 0: 中間結果出力のために戻らない. = 1: 中間結果出力のために反復ごとにIRev = 50で戻る. (他の値であれば省略時の既定値とみなす) [in] MaxIter (省略可) 最大反復回数. (省略時 = 100) (MaxIter <= 0 であれば省略時の既定値とみなす) [out] Iter (省略可) 収束に要した反復回数. [out] RNorm (省略可) 求められた最終近似解における残差ノルム. 出典 SLATEC 使用例 次の非線形連立方程式を解く. x1^2 - x2 - 1 = 0 (x1 - 2)^2 + (x2 - 0.5)^2 - 1 = 0 初期値は, (x1, x2) = (0, 0) とする. Sub Ex_Sos_r() Const N = 2 Dim X(N - 1) As Double, RTolx As Double, ATolx As Double, Tolf As Double Dim Info As Long, K As Long, XX(N - 1) As Double, YY As Double, IRev As Long X(0) = 0: X(1) = 0 RTolx = 0.0000000001 '1.0e-10 ATolx = 0.0000000001 '1.0e-10 Tolf = 0 IRev = 0 Do Call Sos_r(N, X(), RTolx, ATolx, Tolf, Info, K, XX(), YY, IRev) If IRev = 1 Then If K = 1 Then YY = XX(0) ^ 2 - XX(1) - 1 ElseIf K = 2 Then YY = (XX(0) - 2) ^ 2 + (XX(1) - 0.5) ^ 2 - 1 End If End If Loop While IRev <> 0 Debug.Print "X1, X2 =", X(0), X(1) Debug.Print "Info =", Info End Sub 実行結果 X1, X2 = 1.06734608580667 0.13922766688682 Info = 0