Chkder Hybrd Hybrd1 Hybrd1_r Hybrd_r Hybrj Hybrj1 Hybrj1_r Hybrj_r Sos Sos_r

Sos() Sub Sos ( F As  LongPtr, N As  Long, X() As  Double, RTolx As  Double, ATolx As  Double, Tolf As  Double, Info 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) の解をブラウンの方法により求める. アルゴリズムはニュートン法を変形した反復法に基づいており, ガウス−ザイデル法の手順に似たガウス消去法を用いる. 収束はおおむね二次収束である. 本ルーチンの収束特性は方程式の順番に影響を受ける. 線形および弱い非線形の方程式を先に並べるとよい. 関数値を計算するユーザールーチンが必要である. アルゴリズムで必要なすべての偏微分係数は差分商で近似される. 引数 [in] F 関数fi(x)の値を求めるユーザーサブルーチンで, 次のように定義すること. Sub F(N As Long, X() As Double, K As Long, Fval As Double) K = 1〜N の場合: K番目の関数値Fk(X)を求め, Fvalに返すこと. K <= 0 の場合: Nprint = 1 の場合に, 中間結果出力のために反復ごとに呼び出される. -K回目の反復における近似解 X(), 残差ノルム Fvalが渡される. End Sub X()のインデックスは0〜N-1である. 関数値を計算するとき, Kの値は1〜Nであることに注意せよ. Fval以外の変数を変更してはならない. [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に返す) = -2: パラメータ N の誤り. (N < 1) = -3: パラメータ X() の誤り. = -4: パラメータ RTolx の誤り. (RTolx < 0) = -5: パラメータ ATolx の誤り. (ATolx < 0) = -6: パラメータ Tolf の誤り. (Tolf < 0) = 1: 反復回数の上限に達した. (収束が非常に遅い可能性がある. 差許容値が適切であるか確認せよ) = 2: 反復回数の上限に達した. (局所的な最小点にあるか計算精度不足の可能性がある) = 3: 反復が発散した. (残差ノルムが数回の連続した反復において増加した) = 4: ヤコビ行列に関連する行列が特異になった. [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: 中間結果出力のためのFの呼び出し(K <= 0)を行わない. = 1: 中間結果出力のために反復ごとにFを(K <= 0として)呼び出す. (他の値であれば省略時の既定値とみなす) [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 FSos(N As Long, X() As Double, K As Long, Fval As Double) If K = 1 Then Fval = X(0) ^ 2 - X(1) - 1 ElseIf K = 2 Then Fval = (X(0) - 2) ^ 2 + (X(1) - 0.5) ^ 2 - 1 End If End Sub Sub Ex_Sos() Const N = 2 Dim X(N - 1) As Double, RTolx As Double, ATolx As Double, Tolf As Double Dim Info As Long X(0) = 0: X(1) = 0 RTolx = 0.0000000001 '1.0e-10 ATolx = 0.0000000001 '1.0e-10 Tolf = 0 Call Sos(AddressOf FSos, N, X(), RTolx, ATolx, Tolf, Info) Debug.Print "X1, X2 =", X(0), X(1) Debug.Print "Info =", Info End Sub 実行結果 X1, X2 = 1.06734608580667 0.13922766688682 Info = 0