◆ Chkder()

Sub Chkder	(	M As	Long,
		N As	Long,
		X() As	Double,
		Fvec() As	Double,
		Fjac() As	Double,
		Xp() As	Double,
		Fvecp() As	Double,
		Mode As	Long,
		Err() As	Double,
		Info As	Long
	)

微分係数計算値のチェック

目的: 本ルーチンはXにおいて計算されたN変数M非線形連立方程式の微分係数と関数値の整合性をチェックする. Hybrj, Hybrj1などのユーザーサブルーチンをチェックするのに使用することができる.

ユーザーはChkderを, 最初は Mode = 1, 次に Mode = 2 として2回呼び出さなければならない.
Mode = 1:
  入力: Xにチェックを行う点を入れる.
  出力: Xpに近傍の点が設定される.
Mode = 2:
  入力: FvecにXにおける関数値, FjacにXにおける各関数の微分係数値, また, FvecpにXpにおける関数値を入れる.
  出力: Errに各微分係数値の正しさの尺度が入る.

関数値を計算するときに桁落ちあるいは丸め誤差により重大な情報落ちが生じた場合, 本ルーチンは正しく動作しない. 従って, Xの要素のどれも著しく小さい値(特に0)か情報落ちが生じる値であってはならない.

引数

[in]	M	方程式の数. (M > 0)
[in]	N	未知数の数. (N > 0)
[in]	X()	配列 X(LX - 1) (LX >= N) チェックする点の座標.
[in]	Fvec()	配列 Fvec(LFvec - 1) (LFvec >= M) X における関数値 (Mode = 2 の場合).
[in]	Fjac()	配列 Fjac(LFjac1 - 1, LFjac2 - 1) (LFjac1 >= M, LFjac2 >= N) X におけるヤコビ行列(∂fi/∂xj) (i = 0 〜 M-1, j = 0 〜 N-1) の値 (Mode = 2 の場合).
[out]	Xp()	配列 Xp(LXp - 1) (LXp >= N) チェックする点 X の近傍の点の座標 (Mode = 1 の場合).
[in]	Fvecp()	配列 Fvecp(LFvecp - 1) (LFvecp >= M) Xp における関数値 (Mode = 2 の場合).
[in]	Mode	制御変数. (Mode = 1 または 2) 本ルーチンは次のように2回呼び出さなければならない. 1回目: M, N, X()を設定して Mode = 1 として呼び出すとXp()が返る. 2回目: Fvec(), Fjac(), Fvecp()を設定して Mode = 2 として呼び出すと, Err()にそれぞれの式の微分係数の正しさを表す数値が返る.
[out]	Err()	配列 Err(LErr - 1) (LErr >= M) 微分係数のチェック結果 (Mode = 2 の場合). (0 <= Err(i) <= 1) 1 であれば正しく, 0 であれば正しくないことを表す. その中間の値の場合, 確実な判定ができないが, 0.5以上であれば概ね正しいことを表す.
[out]	Info	= 0: 正常終了. = -1: パラメータ M の誤り. (M <= 0) = -2: パラメータ N の誤り. (N <= 0) = -3: パラメータ X() の誤り. = -4: パラメータ Fvec() の誤り. = -5: パラメータ Fjac() の誤り. = -6: パラメータ Xp() の誤り. = -7: パラメータ Fvecp() の誤り. = -8: パラメータ Mode の誤り. (Mode <> 1 かつ Mode <> 2) = -9: パラメータ Err() の誤り.

出典: netlib/minpack

使用例: 次の非線形連立方程式を解くために使ったサブルーチンによる微分係数計算値をチェックする.
x1^2 - x2 - 1 = 0

(x1 - 2)^2 + (x2 - 0.5)^2 - 1 = 0

チェックする点は, (x1, x2) = (1, 1) とする.
Sub FHybrj(N As Long, X() As Double, Fvec() As Double, Fjac() As Double, IFlag As Long)

If IFlag = 1 Then

Fvec(0) = X(0) ^ 2 - X(1) - 1

Fvec(1) = (X(0) - 2) ^ 2 + (X(1) - 0.5) ^ 2 - 1

ElseIf IFlag = 2 Then

Fjac(0, 0) = 2 * X(0)

Fjac(1, 0) = 2 * X(0) - 4

Fjac(0, 1) = -1

Fjac(1, 1) = 2 * X(1) - 1

End If

End Sub

Sub Ex_Chkder()

Const M = 2, N = 2

Dim X(N - 1) As Double, Fvec(N - 1) As Double, Fjac(N - 1, N - 1) As Double

Dim Xp(N - 1) As Double, Fvecp(N - 1) As Double, Err(N - 1) As Double

Dim Mode As Long, Info As Long

X(0) = 1: X(1) = 1

'-- Check 1

Mode = 1

Call Chkder(M, N, X(), Fvec(), Fjac(), Xp(), Fvecp(), Mode, Err(), Info)

'-- Check 2

Mode = 2

Call FHybrj(N, X(), Fvec(), Fjac(), 1)

Call FHybrj(N, X(), Fvec(), Fjac(), 2)

Call FHybrj(N, Xp(), Fvecp(), Fjac(), 1)

Call Chkder(M, N, X(), Fvec(), Fjac(), Xp(), Fvecp(), Mode, Err(), Info)

Debug.Print "Err =", Err(0), Err(1)

Debug.Print "Info =", Info

End Sub

実行結果: Err = 1 0.923072860548928

Info = 0