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.
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, 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 FuncT, TResult delegates
that map a VectorT 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])
};
Dim rosenbrock() As Func(Of Vector(Of Double), Double) =
{
Function(x) 1 - x(0),
Function(x) 10 * (x(1) - x(0) * x(0))
}
No code example is currently available or this language may not be supported.
let rosenbrock : Func<Vector<float>, float>[] =
[|
Func<_,_>(fun x -> 1.0 - x.[0]);
Func<_,_>(fun x -> 10.0 * (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 FuncT1, T2, TResult delegate.
The first argument is a VectorT
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;
}
Shared Function Rosenbrock(x As Vector(Of Double), f As Vector(Of Double)) As Vector(Of Double)
If (f Is Nothing) Then f = Vector.Create(Of Double)(2)
f(0) = 1 - x(0)
f(1) = 10 * (x(1) - x(0) * x(0))
Return f
End Function
No code example is currently available or this language may not be supported.
let Rosenbrock (x : Vector<double>) (f : Vector<double>) =
f.[0] <- 1.0 - x.[0]
f.[1] <- 10.0 * (x.[1] - x.[0] * x.[0])
f
The RightHandSide
property is a VectorT 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 FuncT1, T2, TResult. This is a method that
takes a VectorT argument and returns the matrix result in
a second MatrixT. 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;
}
Function RosenbrockJacobian(x As Vector(Of Double), J As Matrix(Of Double)) As Matrix(Of Double)
If (J Is Nothing) Then J = Matrix.Create(Of Double)(2, 2)
J(0, 0) = -1
J(0, 1) = 0
J(1, 0) = -20 * x(0)
J(1, 1) = 10
Return J
End Function
No code example is currently available or this language may not be supported.
let RosenbrockJacobian (x : Vector<double>) (J : Matrix<double>) =
J.[0, 0] <- -1.0
J.[0, 1] <- 0.0
J.[1, 0] <- -20.0 * x.[0]
J.[1, 1] <- 10.0
J
Alternatively, the individual gradients can be specified using the
SetGradientFunctions
property. This method takes an array of either FuncT, TResult.
or FuncT1, 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)
};
Dim rosenbrockGradients() As Func(Of Vector(Of Double), Vector(Of Double)) =
{
Function(x) Vector.Create(-1.0, 0.0),
Function(x) Vector.Create(-20.0 * x(0), 10)
}
No code example is currently available or this language may not be supported.
let rosenbrockGradients : Func<Vector<float>, Vector<float>>[] = [|
Func<_,_>(fun x -> Vector.Create(-1.0, 0.0) * 1.0);
Func<_,_>(fun x -> Vector.Create(-20.0 * x.[0], 10.0) * 1.0) |]
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);
Dim dogleg As DoglegSystemSolver = New DoglegSystemSolver()
dogleg.InitialGuess = Vector.Create(0.5, 0.5)
dogleg.TargetFunction = AddressOf Rosenbrock
dogleg.SetGradientFunctions(rosenbrockGradients)
No code example is currently available or this language may not be supported.
let dogleg = DoglegSystemSolver()
dogleg.InitialGuess <- Vector.Create(0.5, 0.5)
dogleg.TargetFunction <- Func<_,_,_>(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]));
dogleg.SetSymbolicTargetFunctions(
Function(x) 1 - x(0),
Function(x) 10 * (x(1) - x(0) * x(0)))
No code example is currently available or this language may not be supported.
dogleg.SetSymbolicTargetFunctions(
(fun (x : Vector<float>) -> 1.0 - x.[0]),
(fun (x : Vector<float>) -> 10.0 * (x.[1] - x.[0] * x.[0])))
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 VectorConvergenceTestT 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 VectorConvergenceTestT 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 VectorConvergenceTestT 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;
dogleg.ValueTest.Tolerance = 0.00001
dogleg.ValueTest.Norm = VectorConvergenceNorm.SumOfAbsoluteValues
dogleg.SolutionTest.Enabled = False
No code example is currently available or this language may not be supported.
dogleg.ValueTest.Tolerance <- 0.00001
dogleg.ValueTest.Norm <- VectorConvergenceNorm.SumOfAbsoluteValues
dogleg.SolutionTest.Enabled <- false
You must specify an initial guess for the solution by setting the InitialGuess property to a suitable
VectorT 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 VectorT 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.
The Status property indicates how the algorithm
terminated. Its possible values and their meaning are listed below.
AlgorithmStatus values
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);
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)
No code example is currently available or this language may not be supported.
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)
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.