The Newton Raphson method is a numerical method that in a number of cases can be used for the solution of systems of simultaneous nonlinear equations.

Using the matrix version of a Taylor series, a matrix function can be approximated as shown below.

F_{x1} = F_{x0} + J(x_{1}-x_{0})

Where x_{0} is a column vector of starting points, F_{x0} is the values of the function at the starting points, x_{1} is a column vector of points at some other location, F_{x1} is the approximate function values at that other location, and J is the Jacobian matrix of the system of equations.

The Jacobian matrix is a matrix of the first derivatives of the functions. The functions that are differentiated go down the matrix as rows and the parameters for which the partial derivatives are calculated go across the matrix as columns. An example Jacobian matrix for a system with three functions and three parameters is shown below.

Setting F_{x1} to 0 and rearranging the above equation gives:

The Newton Raphson algorithm works by starting with guess values for the x_{0}s, then the functions are evaluated at that point to get F_{x0}, then the Jacobian is calculated at that point. Using the above equation, x_{1} values are then calculated. These x_{1} values are then used as x_{0} values for a new iteration and the whole process is repeated until either the equations converge or go unstable.

The code for the Newton Raphson class is shown at the bottom of this post.

The Jacobian, function values, and point values are stored in Matrix fields. These are done at the class level to avoid allocating and deallocating the memory for these structures with each iteration.

A field of the class also holds a derivatives member which is used to calculate the partial derivatives of the Jacobian.

The parameters used in the Jacobian are stored in a parameter array and the functions used in the Jacobian are stored in an array of functions of the form double F().

There are two constructors for the class. The first takes the array of functions to use in calculating the Jacobian, a class which derives from, and an integer value indicating the number of points to use when numerically calculating the partial derivatives for the Jacobian. The second constructor is similar to the first except that the number of points to use when numerically calculating the partial derivatives is defaulted to 3. In the constructor, the array of parameters is extracted from the class derived from ModelBase using the GetSolvedForParameters function. The instance of the Derivatives class is then created, and memory is allocated for the Jacobian, function values, and point values.

The function Iterate performs an iteration of the Newton Raphson method. It first calculates the function values at the current point and stores those values and the current parameter values for use in the iteration calculation. It then calculates the Jacobian matrix using the Derivatives class. The new parameter values are then calculated using the Jacobian, function values, and parameter values. The model parameter values are then set to the new calculated values.

using System; using System.Diagnostics; using NumericalMethods.FourthBlog;namespace NumericalMethods {public class NewtonRaphson{private Matrix _jacobian;private Matrix _functionMatrix;private Matrix _x0;private Derivatives _derivatives;private Parameter[] _parameters;private Func<double>[] _functions;public NewtonRaphson(Parameter[] parameters, Func<double>[] functions, int numberOfDerivativePoints){_parameters = parameters;_functions = functions;int numberOfFunctions = _functions.Length;int numberOfParameters = _parameters.Length;Debug.Assert(numberOfParameters == numberOfFunctions);_derivatives = new Derivatives(numberOfDerivativePoints);_jacobian = new Matrix(numberOfFunctions, numberOfParameters);_functionMatrix = new Matrix(numberOfFunctions, 1);_x0 = new Matrix(numberOfFunctions, 1);}public NewtonRaphson(Parameter[] parameters, Func<double>[] functions): this(parameters, functions, 3){}public void Iterate(){int numberOfFunctions = _functions.Length;int numberOfParameters = _parameters.Length;for (int i = 0; i < numberOfFunctions; i++){_functionMatrix[i, 0] = _functions[i]();_x0[i, 0] = _parameters[i];}for (int i = 0; i < numberOfFunctions; i++){for (int j = 0; j < numberOfParameters; j++){_jacobian[i, j] = _derivatives.ComputePartialDerivative(_functions[i], _parameters[j], 1, _functionMatrix[i, 0]);}}Matrix newXs = _x0 - _jacobian.SolveFor(_functionMatrix);for (int i = 0; i < numberOfFunctions; i++){_parameters[i].Value = newXs[i, 0];}}} }