Methods that use the Derivative
The derivative of a function gives more information which can be used to solve equations more efficiently. Two methods that use the derivative are Newton's method and Halley's method.
Convergence of Newton's method is quadratic. This means that once the algorithm is close to the solution, the number of correct digits roughly doubles with each iteration.
Halley's method is a generalization of Newton's method that uses the second derivative of the function. It has cubic convergence, which means that the number of correct digits roughly triples with each iteration. This can make it significantly faster than Newton's method. However, this advantage is often offset by the increased complexity of the calculations.
The Newton Raphson method
Newton devised a method to solve equations using values of the function and its first derivative. It is implemented by the NewtonRaphsonSolver class.
The NewtonRaphsonSolver class has two constructors. The first takes no arguments. The function to solve and its derivative must be set using the TargetFunction and DerivativeOfTargetFunction properties, which are both Func<T, TResult> delegates. The initial guess must be set using the InitialGuess property.
The second constructor takes these three properties as its three parameters.
// Use the parameterless constructor:
NewtonRaphsonSolver solver = new NewtonRaphsonSolver();
// Then set the target function and its derivative:
solver.TargetFunction = Math.Sin;
solver.DerivativeOfTargetFunction = Math.Cos;
// and the initial guess:
solver.InitialGuess = 4;
// Or specify all three directly in the constructor:
NewtonRaphsonSolver solver2 = new NewtonRaphsonSolver(Math.Sin, Math.Cos, 4);The iteration can be controlled as outlined before by setting the AbsoluteTolerance, RelativeTolerance, ConvergenceCriterion, MaxIterations properties. The Solve method calculates the zero of the target function:
Console.WriteLine("Newton-Raphson Solver: sin(x) = 0");
Console.WriteLine(" Initial guess: 4");
double result = solver.Solve();
Console.WriteLine(" Result: {0}", solver.Status);
Console.WriteLine(" Solution: {0}", solver.Result);
Console.WriteLine(" Estimated error: {0}", solver.EstimatedError);
Console.WriteLine(" # iterations: {0}", solver.IterationsNeeded);Halley's method
Edmond Halley, famous for the comet that bears his name, devised a method to solve equations using values of the function and its first and second derivatives. It is implemented by the HalleySolver class.
The HalleySolver class has several constructors. The first takes no arguments. The function to solve and its derivatives must be set using the TargetFunction, DerivativeOfTargetFunction and SecondDerivativeOfTargetFunction properties, which are all Func<T, TResult> delegates.
Alternatively, the TargetFunctionWithDerivatives property can be set to a function that returns the value of the target function and its first and second derivatives as a tuple. This is often preferred because it tends to avoid duplicate calculations.
For example, the Bessel functions satisfy a second order differential equation. This means that the second derivative can be expressed as a relatively simple expression involving the function value and the first derivative.
// Use the parameterless constructor:
HalleySolver solver = new HalleySolver();
// Then set the target function and its derivatives:
solver.TargetFunctionWithDerivatives = x => {
var J0 = Special.BesselJ0(x);
var J1 = Special.BesselJ1(x);
return (J1, J0 - J1 / x, ((2 - x * x) * J1 - x * J0) / (x * x));
};The initial guess must be set using the InitialGuess property.
The iteration can be controlled as outlined before by setting the AbsoluteTolerance, RelativeTolerance, ConvergenceCriterion, MaxIterations properties. The Solve method calculates the zero of the target function:
Console.WriteLine("Halley Solver: J1(x) = 0");
Console.WriteLine(" Initial guess: 4");
double result = solver.Solve();
Console.WriteLine(" Result: {0}", solver.Status);
Console.WriteLine(" Solution: {0}", solver.Result);
Console.WriteLine(" Estimated error: {0}", solver.EstimatedError);
Console.WriteLine(" # iterations: {0}", solver.IterationsNeeded);When you don't have the derivative...
You can still use these classes if you don't have the derivative of the target function. In this case, use the static GetDifferentiator(Func<Double, Double>, Int32) method of the FunctionMath class to create a Func<T, TResult> delegate that represents the numerical derivative of the target function.
It should be noted, however, that this technique requires significantly more function evaluations. If it is possible to find an interval that contains the root, the preferred method is to use one of the root bracketing solvers from the previous section.
A better alternative, if the function is simple enough to express in a single expression, is to use automatic differentiation to have the derivative calculated for you. To do this, you pass a lambda expression that evaluates the function to the GetDerivative(Expression<Func<Double, Double>>) method, which also creates a Func<T, TResult> delegate.
Both techniques are illustrated below:
solver.TargetFunction = Special.BesselJ0;
var derivative1 = FunctionMath.GetNumericalDifferentiator(Special.BesselJ0);
var derivative2 = SymbolicMath.GetDerivative(x => Special.BesselJ0(x));
solver.DerivativeOfTargetFunction = derivative2;
solver.InitialGuess = 5;
Console.WriteLine("Zero of Bessel function near x=5:");
result = solver.Solve();
Console.WriteLine(" Solution: {0}", solver.Result);
Console.WriteLine(" # iterations: {0}", solver.IterationsNeeded);Bounds on the solution
In some cases, a straightforward application of Newton's algorithm can break down. For example, when an approximation ends up in a region where the function is not defined. In such cases, it can be beneficial to define bounds that the algorithm should stay within.
You can do this by setting the LowerBound and UpperBound properties. You can set only one of these properties or both.
The example below finds tries to solve the equation \( \sqrt(1-x)=1 \), starting with \( x = 1 \). Without an upper bound, the algorithm fails with a status of BadFunction. When the upper bound is set to 10, the algorithm finds the solution quickly:
solver.TargetFunction = Special.BesselJ0;
var derivative1 = FunctionMath.GetNumericalDifferentiator(Special.BesselJ0);
var derivative2 = SymbolicMath.GetDerivative(x => Special.BesselJ0(x));
solver.DerivativeOfTargetFunction = derivative2;
solver.InitialGuess = 5;
Console.WriteLine("Zero of Bessel function near x=5:");
result = solver.Solve();
Console.WriteLine(" Solution: {0}", solver.Result);
Console.WriteLine(" # iterations: {0}", solver.IterationsNeeded);Reference
Newton, I. Methodus fluxionum et serierum infinitarum. 1664-1671.
Halley, Edmond Methodus nova accurata & facilis inveniendi radices æqationum quarumcumque generaliter, sine praviæ reductione. 1694.