Sub main() '----------------------------------------------------------------------- ' Demonstrate how to use Newton's routine to solve a set of nonlinear equations. ' Use MNEWT of "Numerical Recipes (1st Edition)" to solve f(x)=0 ' Instructor: Nam Sun Wang '----------------------------------------------------------------------- ' Declare subroutine/function N = 3 Dim xvect(3) ' Read/set the initial guess ------------------------------------------- ' (Since the logarithmic function is involved, 0 is not a good guess.) xvect(1) = Cells(8, 2) xvect(2) = Cells(9, 2) xvect(3) = Cells(10, 2) ' Call a routine from Numerical Recipes to solve f(x)=0 ---------------- ' Usage: CALL MNWET(ntrial ,x, n, tolx, tolf) ' where ' ntrial ... number of tries ' x ... initial guess ' n ... dimension ' tolx ... error in the horizontal (x) direction ' tolf ... error in the vertical (f) direction ' MNEWT calls USRFUN(x(), dfdx(), f()), which the user must provide. Call MNEWT(100, xvect(), N, 0.0001, 0.0001) ' Output results ------------------------------------------------------- X = xvect(1) y = xvect(2) z = xvect(3) Cells(8, 2) = X Cells(9, 2) = y Cells(10, 2) = z ' Double Check ---------------------------------------------------- ' Cells(4, 2) = Sin(X) + y ^ 2 + Log(z) - 7 'log in VBA is the natural log" ' Cells(5, 2) = 3 * X + 2 ^ y - z ^ 3 + 1 ' Cells(6, 2) = X + y + z - 5 End Sub '----------------------------------------------------------------------- Sub usrfun(xvect(), dfdx(), f()) '----------------------------------------------------------------------- ' Note: notation for Version 1.0 has function f and its Jacobian dfdx inverted ' and less number of arguments. ' Usage: Version 1 ' sub usrfun(xvect(), dfdx(), beta()) where dfdx=Jacobian, beta=-f ' Usage: Version 2 ' sub usrfun(xvect(), n, np, f(), dfdx()) ' Variables: ' xvect... vector of independent variables (n variables) ' n ... dimension (unnecessary variable -- I think the authors messed up.) ' np ... physical dimension (needed for the multi-dimensional array dfdx) ' f ... multi-dimensional function (n functions) ' dfdx ... Jacobian matrix ... df/dx (See Eqn 9.6.4 of Numerical Recipes.) ' A set of nonlinear algebraic equations: ' f1(x,y,z) = 0 = sin(x) + y*y + ln(z) - 7. ' f2(x,y,z) = 0 = 3*x + 2^y -z*z*z + 1. ' f3(x,y,z) = 0 = x + y + z - 5. ' The name of the subroutine must be USRFUN, since the name is not provided ' when calling MNEWT. '----------------------------------------------------------------------- ' Note that in the 1st Edition of Numerical Recipes, MNEWT fixes the first ' dimension of alpha to be 15! ' dim xvect(3), dfdx(15,1), f(3) 'No need to dimension because they are arguments ' Use the original variables to reduce confusion ----------------------- X = xvect(1) y = xvect(2) z = xvect(3) ' Specify the function: f(x)=0 ----------------------------------------- f(1) = Sin(X) + y ^ 2 + Log(z) - 7 f(2) = 3 * X + 2 ^ y - z ^ 3 + 1 f(3) = X + y + z - 5 ' Note: Numerical Recipes' Edition 1 is NEGATIVE of the function. f(1) = -f(1) f(2) = -f(2) f(3) = -f(3) ' Define derivatives: dfdx --------------------------------------------- dfdx(1, 1) = Cos(X) dfdx(1, 2) = 2 * y dfdx(1, 3) = 1 / z dfdx(2, 1) = 3 dfdx(2, 2) = Log(2) * 2 ^ y dfdx(2, 3) = -3 * z ^ 2 dfdx(3, 1) = 1 dfdx(3, 2) = 1 dfdx(3, 3) = 1 End Sub '----------------------------------------------------------------------- ' MNEWT of "Numerical Recipes" (Edition 1) ' Need LUDCMP, LUBKSB, and USRFUN (which the user provides) '----------------------------------------------------------------------- Sub MNEWT(NTRIAL, X(), N, TOLX, TOLF) ReDim ALPHA(N, N), BETA(N), INDX(N) For K = 1 To NTRIAL Call usrfun(X(), ALPHA(), BETA()) ERRF = 0! For I = 1 To N ERRF = ERRF + Abs(BETA(I)) Next I If ERRF <= TOLF Then Exit For Call LUDCMP(ALPHA(), N, N, INDX(), D) Call LUBKSB(ALPHA(), N, N, INDX(), BETA()) ERRX = 0! For I = 1 To N ERRX = ERRX + Abs(BETA(I)) X(I) = X(I) + BETA(I) Next I If ERRX <= TOLX Then Exit For Next K Erase INDX, BETA, ALPHA End Sub '----------------------------------------------------------------------- ' LUDCMP is called by MNEWT of "Numerical Recipes" (Edition 1) '----------------------------------------------------------------------- Sub LUDCMP(A(), N, NP, INDX(), D) TINY = 1E-20 ReDim VV(N) D = 1! For I = 1 To N AAMAX = 0! For J = 1 To N If Abs(A(I, J)) > AAMAX Then AAMAX = Abs(A(I, J)) Next J ' If AAMAX = 0! Then Print "Singular matrix.": Exit Sub If AAMAX = 0! Then Call MsgBox("Singular matrix."): Exit Sub 'No printing from Excel VV(I) = 1! / AAMAX Next I For J = 1 To N For I = 1 To J - 1 Sum = A(I, J) For K = 1 To I - 1 Sum = Sum - A(I, K) * A(K, J) Next K A(I, J) = Sum Next I AAMAX = 0! For I = J To N Sum = A(I, J) For K = 1 To J - 1 Sum = Sum - A(I, K) * A(K, J) Next K A(I, J) = Sum DUM = VV(I) * Abs(Sum) If DUM >= AAMAX Then IMAX = I AAMAX = DUM End If Next I If J <> IMAX Then For K = 1 To N DUM = A(IMAX, K) A(IMAX, K) = A(J, K) A(J, K) = DUM Next K D = -D VV(IMAX) = VV(J) End If INDX(J) = IMAX If A(J, J) = 0! Then A(J, J) = TINY If J <> N Then DUM = 1! / A(J, J) For I = J + 1 To N A(I, J) = A(I, J) * DUM Next I End If Next J Erase VV End Sub '----------------------------------------------------------------------- ' LUBKSB is called by MNEWT of "Numerical Recipes" (Edition 1) '----------------------------------------------------------------------- Sub LUBKSB(A(), N, NP, INDX(), B()) II = 0 For I = 1 To N LL = INDX(I) Sum = B(LL) B(LL) = B(I) If II <> 0 Then For J = II To I - 1 Sum = Sum - A(I, J) * B(J) Next J ElseIf Sum <> 0! Then II = I End If B(I) = Sum Next I For I = N To 1 Step -1 Sum = B(I) For J = I + 1 To N Sum = Sum - A(I, J) * B(J) Next J B(I) = Sum / A(I, I) Next I End Sub