Solving Systems of Non-Linear Equations
The EquationSystemSolver class is the abstract base class for all classes that solve systems of nonlinear equations. Each algorithm is implemented by a different class, derived from EquationSystemSolver.
EquationSystemSolver inherits from the ManagedIterativeAlgorithm class. All properties and methods of this class are also available.
Available Algorithms
Newton's Method
Newton's method is a generalization of the Newton-Raphson method for solving equations in one variable to the multi-dimensional case. The Jacobian of the system takes the place of the derivative.
Newton's method is implemented by the NewtonRaphsonSystemSolver class
Powell's Dogleg Method
Powell's dogleg method, also called Powell's hybrid method, attempts to minimize the sum of the squares of the function values. It does this using a combination of Newton's method and the steepest descent method. This is a so-called trust region method. This means that every step moves the current point to within a finite region. This makes the method more stable than Newton's method.
On the other hand, the fact that the method is, in essence, a specialized minimizer means that the algorithm can get stuck in a local minimum that does not correspond to a solution of the system of equations.
The DoglegSystemSolver class implements the dog leg algorithm. This class uses derivatives when available, and uses an optimized variant using numerical derivatives when they are not. It has one additional property, TrustRegionRadius, which allows you to set the initial trust region radius. The default is 1.
Defining the system of equations
The equations that make up the system of equations can be defined in two ways. Viewed as a set of equations, the system is a collection of multivariate functions that are each set equal to a real number (usually 0). The functions must be supplied as an array of Func<T, TResult> delegates that map a Vector<T> to a real number. These can be supplied to the constructor, or set later using the SetTargetFunctions method. The example below defines the system of equations that is the basis for the Rosenbrock function:
1 - x = 0
10(y - x2) = 0
Func<Vector<double>, double>[] rosenbrock =
{
x => 1 - x[0],
x => 10 * (x[1] - x[0] * x[0])
};A system of equations can also be looked at as a single vector equation, where a multivariate function returns a vector that is set equal to a vector (usually all 0's). This function is supplied as one Func<T1, T2, TResult> delegate. The first argument is a Vector<T> that contains the x-values. The second argument is a vector that is to hold the result. The delegate can be passed to the constructor. It can also be set and retrieved through the TargetFunction property. This is the default representation of the system of equations. The example below illustrates how to define the same system of equations using a single delegate that evaluates all equations simultaneously.
Vector<double> Rosenbrock(Vector<double> x, Vector<double> f) {
if (f == null)
f = Vector.Create<double>(2);
f[0] = 1 - x[0];
f[1] = 10 * (x[1] - x[0] * x[0]);
return f;
}The RightHandSide property is a Vector<T> that contains the right-hand sides of the equations. If no value is specified, all right-hand sides are assumed to be zero.
Some algorithms require the Jacobian (matrix of partial derivatives) of the system. Each row of the Jacobian corresponds to the gradient of the corresponding target function. There are two ways to supply this function. The first is through the JacobianFunction property. This property is a Func<T1, T2, TResult>. This is a method that takes a Vector<T> argument and returns the matrix result in a second Matrix<T>. The return value is also a reference to this matrix.
Matrix<double> RosenbrockJacobian(Vector<double> x, Matrix<double> J)
{
if (J == null)
J = Matrix.Create<double>(2, 2);
J[0, 0] = -1;
J[0, 1] = 0;
J[1, 0] = -20 * x[0];
J[1, 1] = 10;
return J;
}Alternatively, the individual gradients can be specified using the SetGradientFunctions property. This method takes an array of either Func<T, TResult>. or Func<T1, T2, TResult> delegates. If the delegate has a second argument, it represents the vector that is to hold the result.
Func<Vector<double>, Vector<double>>[] rosenbrockGradients =
{
x => Vector.Create(-1.0, 0.0),
x => Vector.Create(-20.0 * x[0], 10)
};The following code puts this all together. It creates a DogLegSystemSolver and defines the system of equations using the functions defined in earlier examples:
DoglegSystemSolver dogleg = new DoglegSystemSolver();
dogleg.InitialGuess = Vector.Create(0.5, 0.5);
dogleg.TargetFunction = Rosenbrock;
dogleg.SetGradientFunctions(rosenbrockGradients);If the functions are of a simple enough form, then automatic differentiation can be used. This completely eliminates the need to define derivatives, gradients, or Jacobians by hand. The code is also very concise:
dogleg.SetSymbolicTargetFunctions(
x => 1 - x[0],
x => 10 * (x[1] - x[0] * x[0]));Termination Criteria
Algorithms for solving systems of nonlinear equations typically have two criteria that may each signal that a solution has been found.
The ValueTest property returns a VectorConvergenceTest<T> object that represents the convergence test based on the function values. The test is successful when some norm of the vector of function values is less than the tolerance. By default, the function value with the largest absolute value is used. The test returns Divergent if one of the function values is infinite, and BadFunction when the value is NaN.
The SolutionTest property returns a VectorConvergenceTest<T> object that represents the convergence test based on the estimated solution. The test is successful if the change in the approximate solution is less than the tolerance. The test returns a value of Divergent if one of the elements of the solution is infinite, and BadFunction when one of the values is NaN.
The VectorConvergenceTest<T> has many options to specify the details of the convergence test. By default, the error is calculated as the maximum of the errors in each of the elements of the change of the solution from the previous iteration compared to the corresponding elements of the solution.
The following example specifies that the algorithm is to end when the sum of the absolute values of the function values is less than 0.00001. The solution test is not used:
dogleg.ValueTest.Tolerance = 0.00001;
dogleg.ValueTest.Norm = VectorConvergenceNorm.SumOfAbsoluteValues;
dogleg.SolutionTest.Enabled = false;Calculating the solution
You must specify an initial guess for the solution by setting the InitialGuess property to a suitable Vector<T> value.
The Solve method does the actual work of solving the system of equations. When called without parameters, it returns the approximation to the zero of the target function. When called with a single Vector<T> argument, it attempts to find a point where the target functions equal the corresponding components of the vector.
The Solve method always returns the best estimate for the solution of the system. Successive calls to the Result property will also return this value, until the next call to Solve. An optional parameter lets you specify the right-hand side of the equation TargetFunction(x) = rightHandSide that is to be solved. You can also set this value through the RightHandSide property.
If the ThrowExceptionOnFailure property is set to true, an exception is thrown if the algorithm has failed to converge to a solution within the desired accuracy. If false, the Solve method returns the best approximation to the zero, regardless of whether it is within the requested tolerance.
Interpreting the results
The Status property indicates how the algorithm terminated. Its possible values and their meaning are listed below.
Value | Description |
|---|---|
NoResult | The algorithm has not been executed. |
Converged | The algorithm ended normally. The desired accuracy has been achieved. |
IterationLimitExceeded | The number of iterations needed to achieve the desired accuracy is greater than MaxIterations. |
EvaluationLimitExceeded | The number of function evaluations needed to achieve the desired accuracy is greater than MaxEvaluations. |
RoundOffError | Round-off prevented the algorithm from achieving the desired accuracy. |
BadFunction | Bad behavior of one or more target functions prevented the algorithm from achieving the desired accuracy. |
Divergent | The algorithm converged to a point that is not a solution of the system. |
ConvergedToFalseSolution | The algorithm converged but the solution that was obtained is not a solution to the system of equations. |
The normal status for successful completion is Converged. Of particular interest in this instance is the value ConvergedToFalseSolution. This status is returned when the algorithm converges, but the solution that was found is not a solution of the system of equations. This can happen when an algorithm that solves the system by minimizing the sum of the squares of the function values converges to a local minimum that does not coincide with the global minimum of the sum of squares.
Several properties provide more information about the algorithm. The IterationsNeeded property returns the number of iterations the algorithm needed to complete. The EvaluationsNeeded property returns the number of times the target functions were evaluated. The convergence tests, ValueTest and SolutionTest, each have an Error property that gives an estimate for the error in the solution.
In the next example, we print a summary of the results of the algorithm:
Console.WriteLine(" Initial guess: {0:F2}", dogleg.InitialGuess);
Console.WriteLine(" Status: {0}", dogleg.Status);
Console.WriteLine(" Solution: {0}", dogleg.Result);
Console.WriteLine(" Estimated error: {0}", dogleg.ValueTest.Error);
Console.WriteLine(" # iterations: {0}", dogleg.IterationsNeeded);
Console.WriteLine(" # evaluations: {0}", dogleg.EvaluationsNeeded);References
J.J. Moré, M.Y. Cosnard, Numerical Solution of Nonlinear Equations", ACM Transactions on Mathematical Software , Vol 5, No 1, 1979.
M. J. D. Powell, A Hybrid Method for Nonlinear Equations. Numerical Methods for Nonlinear Algebraic Equations, P. Rabinowitz, editor. Gordon and Breach, 1970.