Now that we are familiar with how Newton’s method works, how can we apply it to a system 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 last chapter,

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 (x0 , y0 ) 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, x0 can be represented as a vector

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

We can represent f(x0) with a vector

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

To represent f'(x0), we will have to use what is called a Jacobian matrix. In this example, we will use a Jacobian matrix of f1 and f2. 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}

We can then rearrange the general equation and substitute the matrices and vectors

$latex  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 xn+1-xn and yn+1-yn, and obtain the values for xn and yn. 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)

From here, we can then calculate the value of x1 and y1, 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 solve for y1 and x1 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 a few more iterations, the guess will be so close to the actual root that the difference is negligible and the guess can be considered the actual 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)
				
			

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)
				
			

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)
				
			

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))))
				
			

The equations each would all return 0 if the values of x,y and 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.