I would greatly appreciate some help in figuring out what I'm doing wrong in trying to get mnewt to run. I'm trying to use it to solve 3 sets of simultaneous nonlinear equations known as the Fleischman equations, which are used to impute specified distributional characteristics on randomly generated data. The 3 equations are as follows:

b2 + 6bd + 2c2 + 15d2 - 1 = 0
2c(b2 + 24bd + 105d2 + 2) - skewness = 0
24[bd + c2(1 + b2 + 28bd) + d2(12 + 48bd + 141c2 + 225d2)] - kurtosis = 0

where skewness and kurtosis are specified by the user. In the code below, I'm using skewness = -1 and kurtosis = 2.5.

The problem is that after no more than 3 iterations the fjac() values get huge and the ludcmp routine crashes with an overflow error in the a() matrix, which holds the latest fjac() upon input. Instead of converging on anything reasonable, the intermediate root estimates just rapidly balloon to astronomical size.

Below is my code in VB6 (which is easily translatable to Fortran or C). I do not supply the code for ludcmp or lubksb since they are unchanged from that given in NR.

Private Sub Main()

Dim d As Integer, i As Integer, j As Integer, k As Integer, np As Integer
Dim f() As Double, Ntrial As Integer, tolf As Double, tolx As Double, errf As Double, _
errx As Double
Dim Abterm As Boolean

tolf = 0.00000001
tolx = 0.00000001

n = 3

np = 3 'this is the number of functions

ReDim f(n), fjac(np, np), fvec(np) As Double, p(np) As Double, x(n) As Double, _
indx(np) As Integer

x(1) = 1
x(2) = 0
x(3) = 0

fvec(1) = X1 ^ 2 + 6 * X1 * x3 + 2 * X2 ^ 2 + 15 * x3 ^ 2 - 1
fvec(2) = 2 * X2 * (X1 ^ 2 + 24 * X1 * x3 + 105 * x3 ^ 2 + 2) + 1
fvec(3) = 24 * (X1 * x3 + X2 ^ 2 * (1 + X1 ^ 2 + 28 * X1 * x3) + x3 ^ 2 * (12 + 48 * X1 * x3 + 141 * X2 ^ 2 + 225 * x3 ^ 2)) - 2.5

Ntrial = 40

Abterm = False

'mnewt(ntrial, x, n, tolx, tolf) routine follows:
'Given an initial guess x for a root in n dimensions, take Ntrial Newton-Raphson steps to improve
'the root. Stop if the root converges in either summed absolute variable increments tolx or
'summed absolute function values tolf.

For k = 1 To Ntrial

Call fdjac(n, x, f, np, fvec, fjac) 'User subroutine supplies function values at x in
'fvec and returns Jacobian matrix in fjac and
'updated f() with x() values from previous iteration.
errf = 0#
For i = 1 To n 'Check function convergence.
errf = errf + Abs(f(i))
Next
If errf <= tolf Then Exit For
For i = 1 To n
p(i) = -f(i)
Next

Call ludcmp(fjac, n, np, indx, d, Abterm) 'Solve linear equations using LU decomposition.
If Abterm = True Then Exit For

Call lubksb(fjac, n, np, indx, p)

errx = 0#
For i = 1 To n
errx = errx + Abs(p(i))
x(i) = x(i) + p(i)
Next
If errx <= tolx Then Exit For

Next

MsgBox "b = " & x(1) & vbCrLf & "c = " & x(2) & vbCrLf & "d = " & x(3) & vbCrLf

End Sub

Private Sub fdjac(n, x, f, np, fvec, fjac) 'I slightly modified the arguments

'Computes forward-difference approximation to Jacobian. On input,
'x(1:n) is the 'point at which the Jacobian is to be evaluated,
'fvec(1:n) is the vector of function values at the point, and
'np is the physical dimension of the Jacobian array fjac '(1:n,1:n) which is output.

'n = number of functions
'x(n) = the next values of each of the n functions

'subroutine funcv(n,x,fvec) is a fixed-name, user-supplied routine that returns the
'vector of functions at x.

'Parameters: NMAX is the maximum value of n; EPS is the approximate square root of the
'machine precision.

Dim h As Double, temp As Double, EPS As Double
ReDim temparray(n)
EPS = 0.0001

'Compute fvec() with current x() values; Note that the x() values are updated prior to each call to fdjac
Call Funcv(x, n, f)
For j = 1 To n
fvec(j) = f(j)
Next

For j = 1 To n
temp = x(j)
h = EPS * Abs(temp)
If (h <= EPS) Then h = EPS
x(j) = temp + h 'Trick to reduce finite precision error.
h = x(j) - temp

Call Funcv(x, n, f) 'computes f with forward differenced x values

x(j) = temp
For j = 1 To n
For i = 1 To n 'Forward difference formula.
fjac(i, j) = (f(i) - fvec(i)) / h
Next
Next

End Sub

Private Sub Funcv(x, n, f)

f(1) = x(1) ^ 2 + 6 * x(1) * x(3) + 2 * x(2) ^ 2 + 15 * x(3) ^ 2 - 1
f(2) = 2 * x(2) * (x(1) ^ 2 + 24 * x(1) * x(3) + 105 * x(3) ^ 2 + 2) + 1
f(3) = (24 * (x(1) * x(3) + x(2) ^ 2 * (1 + x(1) ^ 2 + 28 * x(1) * x(3)) + x(3) ^ 2 * (12 + 48 * x(1) * x(3) + 141 * x(2) ^ 2 + 225 * x(3) ^ 2)) - 2.5)

End Sub

I found one error: the initial setting of the fvec() values at the beginning of the main program did not use array syntax for the X values. I changed this, but it made no difference: same progressive increase in error until overflow is reached in ludcmp.