Solving Non-linear System of Equations with Newton's Method

Introduction to Data Science Finding Roots of System of Equations Solving Non-linear System of Equations with Newton’s Method

As systems of nonlinear equations becomes more complex, it will be harder to obtain solutions by substitutions or graphing. Newton’s method allows us to quickly solve for the system. It is also the method that we will implement in solving using NMdev.

What is Newton’s Method?

The newton’s method, also known as the Newton–Raphson method, is a root-finding algorithm which produces successively closer approximations to the actual roots of a real-valued function. The basic concept is to start with an initial guess which is close to the true root, $x_{0} $, then to approximate the function by its tangents, and then compute the x-intercept of this tangent line. Suppose we need to find the roots this function $f $, where the equation of the graph is

$f(x)=x^3 +3x^2+1 $

From the graph, a good close initial guess would be $-4 $. With that, we can draw the tangent line and find the x-intercept of the tangent. This x-intercept will typically be a better approximation to the original function’s root than the first guess. By Taylor series expansion,

$f(x)=f(a)+f^{\prime}(a)(x-a)+\frac{f^{\prime\prime}(a)(x-a)^{2}}{2}+ $ . . .

We can substitute in $x $ as the solution and $a $ as the guess, $x_{0} $.

$f(x)=f(x_{0})+f^{\prime}(x_{0})(x-x_{0})+\frac{f^{\prime\prime}(x_{0})(x-x_{0})^{2}}{2}+ $ . . .

If $x_{0} $ is sufficiently close to $x $, $(x_{0}-x) $ is small, and higher powers of $(x_{0}-x) $ can be ignored. We can then rewrite the equation, substituting $x $ as $x_{1} $ which is the value of $x $ that is comparatively to the solution, and let $f(x)=0 $ since we are finding the roots. We will come to this equation

$0=f(x_{0})+f'(x_{0})(x_{1}-x_{0}) $

Solving for $x_{1} $ gives

$f'(x_{0})(x_{1})=f'(x_{0})(x_{0})-f(x_{0}) $

$x_{1}=x_{0}-\frac{f(x_{0})}{f'(x_{0})} $

The equation can be rewritten in a general form as a iterative equation,

$x_{n+1}=x_{n}-\frac{f(x_{n})}{f'(x_{n})} $

where $x_{n+1} $ is a value that is more accurate than $x_{n} $, provided this initial guess is close enough to the unknown zero, and that $f'(x_{0})\neq0 $.Using this equation of the graph, we are now able to deduce the value of $x_{1} $

$x_{1}=-4-f(-4)/(f'(-4)) $

$f(-4)=(-4)^3+3(-4)^2+1=-15 $

$f'(-4)=3(-4)^2+6(-4)+1=25 $

$x_{1}=-4+15/25=-3.4 $

The idea behind Newton’s system is that each x-intercept that we find with our tangent will be closer to the root than our initial guess. Graphically, when we plot the values and the tangent of $x_{0} $ on the graph, we will find that the new x-intercept, $x_{1} $, of this tangent is closer to the root than the initial value of $x $ that we guessed.

We can then continue to repeat this process iteratively until we arrive at a value that is precise enough. Continuing with the process using $x_{1} $, we can obtain the value $x_{2} $.

$x_{1}=-3.4-f(-3.4)/(f'(-3.4)) $

$f(-3.4)=(-3.4)^3+3(-3.4)^2+1=-3.624 $

$f'(-3.4)=3(-3.4)^2+6(-3.4)+1=15.28 $

$x_{1}=-3.4+3.624/15.28=-3.1628 $

Graphically, we can show that $x_{2} $, is closer to $x_{1} $ by plotting the tangent at $x=x_{1} $.

We will keep iterating Newton’s method until the precision threshold is satisfied or the maximum number of iterations is reached. The difference of $x $, to the actual root by the time the precision threshold is satisfied should be negligible, and hence with Newton’s root method, we are able to find the roots of a nonlinear equation.

Now, we will look at how to apply Newton’s method to solve a system of non-linear equations, as well as how use the NewtonSystemRoot class in NMdev to implement a solver for systems of non-linear equations.

Solving non-linear system of equations using Newton’s method

Recall the general solution of the Newton’s method that we had deduced,

$x_{n+1}=x_{n}-\frac{f(x_{n})}{f'(x_{n})} $

We are able to use the same equation for a system of non-linear equations, but instead of using $f(x) $ and $f'(x) $, we can substitute them as matrices. Suppose we are finding the solution to a system of 2 equations and 2 unknowns,

$f_{1}:3x + y^2-12=0 $
$f_{2}:x^2 +y-4=0 $

We begin by choosing an initial guess $(x_{0},y_{0}) $ near the intersection of both equations.

Looking at the graph, there are two intersections. Let us guess the lower intersection first. A reasonable guess would be $(2,-2) $. The initial guess, $x_{0} $ can be represented as a vector

$x_{0}=\left(\begin{array}{c}2\\-2\end{array}\right) $

We can represent $f(x_{0}) $ with a vector

$f(x_{0})=\left(\begin{array}{c}f_{1}(2,-2)\\f_{2}(2,-2)\end{array}\right) $

To represent $f'(x_{0}) $, we will have to use what is called a Jacobian matrix. The Jacobian matrix is the matrix of all the equations’ first-order partial derivatives. It acts as a gradient for a system of multivariate equations. A Jacobian matrix for 2 equations and 2 unknowns will look like this:

$\begin{bmatrix}\frac{\partial f_{1}}{\partial x} & \frac{\partial f_{1}}{\partial y} \\\frac{\partial f_{2}}{\partial x} & \frac{\partial f_{2}}{\partial y} \end{bmatrix} $

In the case of our example, we can find the Jacobian matrix by partially differentiating the 2 equations.

$\frac{\partial f_{1}}{\partial x}=3 $

$\frac{\partial f_{1}}{\partial y}=2y $

$\frac{\partial f_{2}}{\partial x}=2x $

$\frac{\partial f_{2}}{\partial y}=1 $

The Jacobian matrix will for our example will thus be

$J(f_{1},f_{2})=\begin{bmatrix}3&2y\\2x&1\end{bmatrix} $

$J(f_{1},f_{2})=\begin{bmatrix}3&-4\\4&1\end{bmatrix} $

With the Jacobian matrix, we can then rearrange the general equation using the matrices and vectors

$x_{n+1}=x_{n}-\frac{f(x_{n})}{f'(x_{n})} $

$\left(\begin{array}{c}x_{n+1}\\y_{n+1}\end{array}\right)=\left(\begin{array}{c}x_{n}\\y_{n}\end{array}\right)-\frac{\left(\begin{array}{c}f(x_{1},y_{1})\\ f(x_{2},y_{2})\end{array}\right)}{J^{-1}(f_{1},f_{2})} $

Since $J^{-1} $ is harder to caculate, we will just bring it over, and we will have our general equation for solving non-linear systems of equations:

$J(f_{1},f_{2})\left(\begin{array}{c}x_{n+1}-x_{n}\\y_{n+1}-y_{n}\end{array}\right)=-\left(\begin{array}{c}f(x_{1},y_{1})\\ f(x_{2},y_{2})\end{array}\right) $

From here, we will be able to solve for $x_{n+1}-x_{n} $ and $y_{n+1}-y_{n} $, and obtain the values for $x_{n} $ and $y_{n} $. Using our example, we will obtain this equation

$\begin{bmatrix}3&-4\\4&1\end{bmatrix}\left(\begin{array}{c}x_{1}-2\\y_{1}+2\end{array}\right)=-\left(\begin{array}{c}f_{1}(2,-2)\\f_{2}(2,-2)\end{array}\right) $

$f_{1}(2,-2)=3(2)+4-12=-2 $

$f_{2}(2,-2)=4-2-4=-2 $

$\begin{bmatrix}3&-4\\4&1\end{bmatrix}\left(\begin{array}{c}x_{1}-2\\y_{1}+2\end{array}\right)=-\left(\begin{array}{c}-2\\-2\end{array}\right) $

Now, we can then calculate the value of $x_{1} $ and $y_{1} $, which are closer estimates to the intersections of the 2 functions than our previous guess.

$\begin{bmatrix}3(x_{1}-2)-4(y_{1}+2)\\4(x_{1}-2)+y_{1}+2\end{bmatrix}=\left(\begin{array}{c}2\\2\end{array}\right) $

$\begin{bmatrix}3x_{1}-6-4y_{1}-8\\4x_{1}-8+y_{1}+2\end{bmatrix}=\left(\begin{array}{c}2\\2\end{array}\right) $

We can then solve for $x_{1} $ and $y_{1} $ using Gaussian Elimination.

$\begin{bmatrix}3x_{1}-4y_{1}-14\\4x_{1}+y_{1}-6\end{bmatrix}=\left(\begin{array}{c}2\\2\end{array}\right) $

$\begin{bmatrix}3x_{1}-4y_{1}\\4x_{1}+y_{1}\end{bmatrix}=\left(\begin{array}{c}16\\8\end{array}\right) $

$\begin{bmatrix}3x_{1}-4y_{1}\\4x_{1}+y_{1}+3/4x_{1}-y_{1}\end{bmatrix}=\left(\begin{array}{c}16\\12\end{array}\right) $

$\begin{bmatrix}-4y_{1}+3x_{1}\\0+19/4x_{1}\end{bmatrix}=\left(\begin{array}{c}16\\12\end{array}\right) $

$x_{1}=2.5263 $

$y_{1}=-2.1053 $

Hence, after one iteration, we can see from the graph that our new guess $(2.5263,-2.1053) $ is already much closer towards the actual root than the original guess of $(2,-2) $, and after more iterations, the difference between the real answer and the guess is negligible and the guess can be considered the solution.

In NM Dev, we can use the class NewtonSystemRoot to implement newton’s method to solve for a system of equations. We will apply the same system

$f_{1}:3x + y^2-12=0 $

$f_{2}:x^2 +y-4=0 $

to demonstrate as an example.

				
					%use s2
//First let us create the two functions f1 and f2
val f1:BivariateRealFunction = object : AbstractBivariateRealFunction() {
    //Using AbstractBivariateRealFunction as we are only using 2 variables
    public override fun evaluate(x: Double, y: Double): Double {
        //represents 3x+y^2-12=0
        //if the values are correct, it returns 0
        return 3 * x + y * y - 12
}}

val f2:BivariateRealFunction = object : AbstractBivariateRealFunction() {
    //Using AbstractBivariateRealFunction as we are only using 2 variables
    public override fun evaluate(x: Double, y: Double): Double {
        //represents x^2+y-4=0
        //if the values are correct, it returns 0
        return x * x + y - 4
}}
// public NewtonSystemRoot(double accuracy, int maxIter) takes in 2 parameters:
// @param accuracy the terminating accuracy
// @param maxIter the maximum number of iterations

// NewtonSystemRoot has a function solve(RealVectorFunction f,Vector guess)
// @param f is a multivariate function
// @param guess is an initial guess of the root
// It returns an approximate root
// and throws NoRootFoundException when the search fails to find a root

//creating a NewtonSystemRoot solver with tolerance of 1e-8
//and max iterations of 1
val solver:NewtonSystemRoot = NewtonSystemRoot(1e-8, 1)
//This is our initial guess
val initial:Vector = DenseVector(arrayOf(2.0,-2.0))
//We use the solve function to return the root and print the output
val root:Vector = solver.solve(arrayOf(f1,f2),initial)
print(root)

output: [2.526316, -2.105263]

The output is the same as what we had gotten by performing the Newton’s method manually. In the above code, we set max iterations to 1 to compare with our workings, which meant that we did one iteration of the Newton’s method. However, to get to more accurate results, we can set max iterations to a higher number, so that our answer will be even closer to the actual roots. In this case, we set our max iterations to 10 while using the same system of equations.

				
					//creating a NewtonSystemRoot solver with tolerance of 1e-8
//and max iterations of 10
val solver:NewtonSystemRoot = NewtonSystemRoot(1e-8, 10)
//This is our initial guess
val initial:Vector = DenseVector(arrayOf(2.0,-2.0))
//We use the solve function to return the root and print the output
val root:Vector = solver.solve(arrayOf(f1,f2),initial)
print(root)

output: [2.477352, -2.137275]

After 10 iterations, this answer can now be said to be the solution of this system, as the difference between the actual solution and these values are so insignificant that it can be considered negligible.

You can use the same concept to apply to a system of non-linear equations with 3 or more unknowns. Let’s use this system as an example,

$x^2+y^3-z=6 $

$2x+9y-z=17 $

$x^4 + 5y + 6z = 29 $

Let us use AbstractRealVectorFunction instead to represent the system, to make it easier to read

				
					%use s2
val G:RealVectorFunction = object : AbstractRealVectorFunction(3, 3) {
    public override fun evaluate(v:Vector) : Vector {
        val x:Double = v.get(1)
        val y:Double = v.get(2)
        val z:Double = v.get(3)
        //first equation of the system
        //x^2 + y^3 − z = 6
        val g1:Double = Math.pow(x, 2.0) + Math.pow(y, 3.0) - z - 6.0
        //second equation of the system
        //2x + 9y − z = 17
        val g2:Double = 2.0 * x + 9.0 * y - z - 17.0
        //third equation of the system
        //x^4 + 5y + 6z = 29
        val g3:Double = Math.pow(x, 4.0) + 5.0 * y + 6.0 * z - 29.0
        
        val g:Vector = DenseVector(arrayOf(g1, g2, g3))
        return g
    }
}
//creating a NewtonSystemRoot solver with tolerance of 1e-8
//and max iterations of 20
val solver:NewtonSystemRoot = NewtonSystemRoot(1e-8, 20)
//This is our initial guess
val initial:Vector = DenseVector(arrayOf(0.0,0.0,0.0))
//We use the solve function to return the root and print the output
val root:Vector = solver.solve(G,initial)
print(root)

output:[-4.580190, -4.307120, -64.924460]

The solution of

$x=-4.580190 $

$y=-4.307120 $

$z=-64.924460 $

can be proven to be correct by substituting the values back into function G.

				
					print(G.evaluate(DenseVector(arrayOf(-4.580190, -4.307120, -64.924460))))

output:[-0.000000, 0.000000, 0.000016]

The equations each would all return 0 if the values of $x,y and $latex z $ exactly equate to the root of the system of non-linear equations. Though the last equation did not return an exact zero, the output is still close enough for our solution of $x, y $ and $z $ to be considered the roots to the system.

Previous Topic

Back to Lesson

Next Lesson

Solving Non-linear System of Equations with Newton’s Method

Did you know?