Option Explicit Function MyFx(ByVal sfX As String, ByVal X As Double) As Double sfX = Replace(sfX, "$X", "(" & CStr(X) & ")") MyFx = Evaluate(sfX) End Function Function MySlope(ByVal sfX As String, ByVal X As Double, ByVal Fx As Double) As Double Dim h As Double h = 0.001 * (1 + Abs(X)) MySlope = (MyFx(sfX, X + h) - Fx) / h End Function Function MySecDeriv(ByVal sfX As String, ByVal X As Double) As Double Dim h As Double, Fp As Double, Fx As Double, Fm As Double h = 0.001 * (1 + Abs(X)) Fx = MyFx(sfX, X) Fp = MyFx(sfX, X + h) Fm = MyFx(sfX, X - h) MySecDeriv = (Fp - 2 * Fx + Fm) / h / h End Function Function NewtonRoot(ByVal sfX As String, ByVal X As Double, ByVal Toler As Double) As Double Dim h As Double, Fx As Double, Diff As Double Do h = 0.001 * (1 + Abs(X)) Fx = MyFx(sfX, X) Diff = h * Fx / (MyFx(sfX, X + h) - Fx) X = X - Diff Loop Until Abs(Diff) < Toler NewtonRoot = X End Function Function NewtonRootEx(ByVal sfX As String, ByVal A As Double, ByVal B As Double, ByVal Toler As Double) As Double Dim X As Double, h As Double, Fx As Double, Diff As Double X = (A + B) / 2 Do h = 0.001 * (1 + Abs(X)) Fx = MyFx(sfX, X) Diff = h * Fx / (MyFx(sfX, X + h) - Fx) X = X - Diff If X < A Or X > B Then X = (A + B) / 2 If MyFx(sfX, A) * MyFx(sfX, X) > 0 Then A = X Else B = X End If Diff = B - A End If Loop Until Abs(Diff) <= Toler NewtonRootEx = X End Function Function NewtonMinMax(ByVal sfX As String, ByVal X As Double, ByVal Toler As Double) As Double Dim h As Double, Fx As Double, Diff As Double Dim Fp As Double, Fm As Double, Drv1 As Double, Drv2 As Double Do h = 0.001 * (1 + Abs(X)) Fx = MyFx(sfX, X) Fp = MyFx(sfX, X + h) Fm = MyFx(sfX, X - h) Drv1 = (Fp - Fm) / 2 / h Drv2 = (Fp - 2 * Fx + Fm) / h / h Diff = Drv1 / Drv2 X = X - Diff Loop Until Abs(Diff) < Toler NewtonMinMax = X End Function Function NewtonMinMaxEx(ByVal sfX As String, ByVal A As Double, ByVal B As Double, _ ByVal Toler As Double) As Double Dim h As Double, Fx As Double, Fa As Double, Diff As Double Dim Fp As Double, Fm As Double, Drv1 As Double, Drv2 As Double X = (A + B) / 2 Do h = 0.001 * (1 + Abs(X)) Fx = MyFx(sfX, X) Fp = MyFx(sfX, X + h) Fm = MyFx(sfX, X - h) Drv1 = (Fp - Fm) / 2 / h Drv2 = (Fp - 2 * Fx + Fm) / h / h Diff = Drv1 / Drv2 X = X - Diff If X < A Or X > B Then X = (A + B) / 2 Fx = MyFx(sfX, X) Fa = MyFx(sfX, A) If MySlope(sfX, A, Fa) * MySlope(sfX, X, Fx) > 0 Then A = X Else B = X End If Diff = B - A End If Loop Until Abs(Diff) < Toler NewtonMinMaxEx = X End Function Function QuartileRoot(ByVal sfX As String, ByVal A As Double, ByVal B As Double, _ ByVal Toler As Double) As Double Const ALPHA = 0.275 Dim Fa As Double, Fb As Double, Fx As Double, X As Double Fa = MyFx(sfX, A) Fb = MyFx(sfX, B) Do If Abs(Fa) < Abs(Fb) Then X = A + ALPHA * (B - A) Else X = A + (1 - ALPHA) * (B - A) End If Fx = MyFx(sfX, X) If Fx * Fa > 0 Then A = X Fa = Fx Else B = X Fb = Fx End If Loop Until Abs(B - A) <= Toler QuartileRoot = (A + B) / 2 End Function Sub ScanRoot1() ' Scan a range and uses Newton's method to find a root when the sign of the ' funtion changes. Const MAX_ERROR = 5 Dim sfX As String Dim R As Integer, NumSteps As Long, ErrCount As Integer Dim A As Double, B As Double, Stp As Double, Toler As Double, FxToler As Double Dim Xa As Double, Fa As Double, SFa As Integer, Da As Double, SDa As Integer Dim Xb As Double, Fb As Double, SFb As Integer, Db As Double, SDb As Integer Dim Fx As Double, h As Double, Diff As Double, X As Double, Drv2 As Double Dim bMoveOneExtraStep As Boolean sfX = UCase([A2].Value) A = [A4].Value B = [A6].Value Stp = [A8].Value Toler = [A10].Value FxToler = [A12].Value Range("B2:D1000").Value = "" Xa = A Fa = MyFx(sfX, Xa) Da = MySlope(sfX, Xa, Fa) SFa = Sgn(Fa) SDa = Sgn(Da) R = 2 NumSteps = 0 ErrCount = 0 On Error GoTo ErrHandler Do NumSteps = NumSteps + 1 bMoveOneExtraStep = False Xb = A + NumSteps * Stp Fb = MyFx(sfX, Xb) Db = MySlope(sfX, Xb, Fb) SFb = Sgn(Fb) SDb = Sgn(Db) ' landed on a root?? If Abs(Fb) = 0 Then Cells(R, 2) = Xb Cells(R, 3) = 0 Drv2 = MySecDeriv(sfX, Xb) If SDa * SDb < 0 Then If Abs(Drv2) > FxToler Then If Drv2 > 0 Then Cells(R, 4) = "Root & Minimum" Else Cells(R, 4) = "Root & Maximum" End If Else Cells(R, 4) = "Root & Saddle point" End If End If R = R + 1 bMoveOneExtraStep = True ' Found a range that contains a root? ElseIf SFb * SFa < 0 Then X = NewtonRoot(sfX, (Xa + Xb) / 2, Toler) Cells(R, 2) = X Cells(R, 3) = MyFx(sfX, X) R = R + 1 bMoveOneExtraStep = True ' located a range that has a minimum/maximum/root? ElseIf SFa * SFb > 0 And SDa * SDb < 0 Then X = NewtonMinMax(sfX, (Xa + Xb) / 2, Toler) Drv2 = MySecDeriv(sfX, X) ' found a root? Cells(R, 2) = X Cells(R, 3) = MyFx(sfX, X) bMoveOneExtraStep = True If Abs(Fx) < FxToler Then Cells(R, 4) = "Root & " Else Cells(R, 4) = "" End If If Drv2 > 0 Then Cells(R, 4) = Cells(R, 4) & "Minimum" Else Cells(R, 4) = Cells(R, 4) & "Maximum" End If R = R + 1 End If ResumeHere: If bMoveOneExtraStep Then NumSteps = NumSteps + 1 Xa = A + NumSteps * Stp Fa = MyFx(sfX, Xa) Da = MySlope(sfX, Xa, Fa) SFa = Sgn(Fa) SDa = Sgn(Da) Else Xa = Xb Fa = Fb Da = Db SFa = SFb SDa = SDb End If Loop Until Xa >= B ExitProc: Exit Sub ErrHandler: ErrCount = ErrCount + 1 If ErrCount > MAX_ERROR Then MsgBox "Reached maximum limit or runtime errors!" Resume ExitProc Else bMoveOneExtraStep = True Resume ResumeHere End If End Sub