◆ Hybrj1_r()

Sub Hybrj1_r	(	N As	Long,
		X() As	Double,
		Fvec() As	Double,
		XTol As	Double,
		Info As	Long,
		XX() As	Double,
		YY() As	Double,
		YYpd() As	Double,
		IRev As	Long
	)

Solution of a system of nonlinear equations by Powell hybrid method (simple driver) (reverse communication version)

Purpose: This routine finds a zero of a system of n nonlinear functions in n variables
fi(x1, x2, ..., xn) = 0 (i = 1 to n)

by a modification of the Powell hybrid method.

The user must provide the calculated function values and the Jacobian according to IRev.

Hybrj1_r is the simple driver for standard use, which calls the Hybrj_r with default parameters.

Parameters

[in]	N	Number of functions and variables. (N > 0)
[in,out]	X()	Array X(LX - 1) (LX >= N) [in] An initial estimate of the solution vector. [out] IRev = 0: The obtained solution vector. IRev = 2: The abscissa where the Jacobian shoule be evaluated.
[out]	Fvec()	Array Fvec(LFvec - 1) (LFvec >= N) The function values evaluated at the solution vector X().
[in]	XTol	Target relative tolerance. Termination occurs when the relative error between two consecutive iterations is at most XTol. (XTol >= 0)
[out]	Info	= 0: Successful exit. (Relative error between two consecutive iterates is at most XTol) = -1: The argument N had an illegal value. (N <= 0) = -2: The argument X() is invalid. = -3: The argument Fvec() is invalid. = -4: The argument XTol had an illegal value. (XTol < 0) = -6: The argument XX() is invalid. = -7: The argument YY() is invalid. = -8: The argument YYpd() is invalid. = 1: Number of calls to F with IFlag = 1 has reached the limit (100*(N + 1). = 2: XTol is too small. No further improvement in the approximate solution X is possible. = 3: Iteration is not making good progress, as measured by the improvement from the last five Jacobian evaluations. = 4: Iteration is not making good progress, as measured by the improvement from the last ten iterations.
[out]	XX()	Array XX(LXX - 1) (LXX >= N) When returned with IRev = 1, XX() contains the abscissa where the function value should be evaluated and given in the next call.
[in]	YY()	Array YY(LYY - 1) (LYY >= N) When returned with IRev=1, the function value fi(XX()) (i = 1 to N) should be given in YY() in the next call.
[in,out]	YYpd()	Array YYpd(LYYpd1 - 1, LYYpd2 - 1) (LYYpd1 >= N, LYYpd2 >= N) [in] IRev = 2: Jacobian(dfi/dxj) at X() should be given in YYpd() in the next call. [out] IRev = 0: The orthogonal matrix Q produced by the QR factorization of the final approximate Jacobian.
[in,out]	IRev	Control variable for reverse communication. [in] Before first call, IRev should be initialized to zero. On succeeding calls, IRev should not be altered. [out] If IRev is not zero, complete the following process and call this routine again. = 0: Computation finished. See return code in Info. = 1: User should set the function values at XX() in YY(). Do not alter any variables other than YY(). = 2: User should set the Jacobian at X() in YYpd(). Do not alter any variables other than YYpd().

Reference: netlib/minpack

Example Program: Solve the following system of nonlinear equations.
x1^2 - x2 - 1 = 0

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

The initial approximation (x1, x2) = (0, 0) is used.
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_Hybrj1_r()

Const N = 2

Dim X(N - 1) As Double, Fvec(N - 1) As Double, XTol As Double, Info As Long

Dim XX(N - 1) As Double, YY(N - 1) As Double, YYpd(N - 1, N - 1) As Double

Dim IRev As Long

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

XTol = 0.0000000001 '1.0e-10

IRev = 0

Do

Call Hybrj1_r(N, X(), Fvec(), XTol, Info, XX(), YY(), YYpd(), IRev)

If IRev = 1 Or IRev = 2 Then

Call FHybrj(N, XX(), YY(), YYpd(), IRev)

End If

Loop While IRev <> 0

Debug.Print "X1, X2 =", X(0), X(1)

Debug.Print "Info =", Info

End Sub

Example Results: X1, X2 = 1.06734608580669 0.139227666886861

Info = 0