## Nonlinear Systems of Equations: Newton-Raphson Method

### Newton-Raphson Method

The Newton-Raphson method is the method of choice for solving nonlinear systems of equations. Many engineering software packages (especially finite element analysis software) that solve nonlinear systems of equations use the Newton-Raphson method. The derivation of the method for nonlinear systems is very similar to the one-dimensional version in the root finding section. Assume a nonlinear system of equations of the form:

If the components of one iteration are known as: , then, the Taylor expansion of the first equation around these components is given by:

Applying the Taylor expansion in the same manner for , we obtained the following system of linear equations with the unknowns being the components of the vector :

By setting the left hand side to zero (which is the desired value for the functions , then, the system can be written as:

Setting , the above equation can be written in matrix form as follows:

where is an matrix, is a vector of components and is an -dimensional vector with the components . If is invertible, then, the above system can be solved as follows:

#### Example

Use the Newton-Raphson method with to find the solution to the following nonlinear system of equations:

##### Solution

In addition to requiring an initial guess, the Newton-Raphson method requires evaluating the derivatives of the functions and . If , then it has the following form:

Assuming an initial guess of and , then the vector and the matrix have components:

The components of the vector can be computed as follows:

Therefore, the new estimates for and are:

The approximate relative error is given by:

For the second iteration the vector and the matrix have components:

The components of the vector can be computed as follows:

Therefore, the new estimates for and are:

The approximate relative error is given by:

For the third iteration the vector and the matrix have components:

The components of the vector can be computed as follows:

Therefore, the new estimates for and are:

The approximate relative error is given by:

Finally, for the fourth iteration the vector and the matrix have components:

The components of the vector can be computed as follows:

Therefore, the new estimates for and are:

The approximate relative error is given by:

Therefore, convergence is achieved after 4 iterations which is much faster than the 9 iterations in the fixed-point iteration method. The following is the Microsoft Excel table showing the values generated in every iteration:

The Newton-Raphson method can be implemented directly in Mathematica using the “FindRoot” built-in function:

View Mathematica CodeFindRoot[{x1^2 + x1*x2 - 10 == 0, x2 + 3*x1*x2^2 == 57}, {{x1, 1.5}, {x2, 3.5}}]

# Using Sympy import sympy as sp x1, x2 = sp.symbols('x1 x2') sol = sp.nsolve((x1**2 + x1*x2 - 10, x2 + 3*x1*x2**2 - 57), (x1,x2), (1.5,3.5)) print("Sympy: ",sol) # Using Scipy's root from scipy.optimize import root def equations(x): x1, x2 = x return [x1**2 + x1*x2 - 10, x2 + 3*x1*x2**2 - 57] sol = root(equations, [1.5,3.5]) print("Scipy: ",sol.x)

MATLAB code for the Newton-Raphson example

MATLAB file 1 (nrsystem.m)

MATLAB file 2 (f.m)

MATLAB file 3 (K.m)

The following “While” loop in Mathematica provides an iterative procedure that implements the Newton-Raphson Method in Mathematica and produces the table shown:

View Mathematica Codex = {x1, x2}; f1 = x1^2 + x1*x2 - 10; f2 = x2 + 3 x1*x2^2 - 57; f = {f1, f2}; K = Table[D[f[[i]], x[[j]]], {i, 1, 2}, {j, 1, 2}]; K // MatrixForm xtable = {{1.5, 3.5}}; xrule = {{x1 -> xtable[[1, 1]], x2 -> xtable[[1, 2]]}}; Ertable = {1} NMax = 100; i = 1; eps = 0.000001 While[And[Ertable[[i]] > eps, i <= NMax], Delta = -(Inverse[K].f) /. xrule[[i]]; xtable = Append[xtable, Delta + xtable[[Length[xtable]]]]; xrule = Append[ xrule, {x1 -> xtable[[Length[xtable], 1]], x2 -> xtable[[Length[xtable], 2]]}]; er = Norm[Delta]/Norm[xtable[[i+1]]]; Ertable = Append[Ertable, er]; i++] Title = {"Iteration", "x1_in", "x2_in", "x1_out", "x2_out", "er"}; T2 = Table[{i, xtable[[i, 1]], xtable[[i, 2]], xtable[[i + 1, 1]], xtable[[i + 1, 2]], Ertable[[i + 1]]}, {i, 1, Length[xtable] - 1}]; T2 = Prepend[T2, Title]; T2 // MatrixForm

import numpy as np import sympy as sp import pandas as pd x1, x2 = sp.symbols('x1 x2') x = [x1, x2] f1 = x1**2 + x1*x2 - 10 f2 = x2 + 3*x1*x2**2 - 57 f = sp.Matrix([f1, f2]) K = sp.Matrix([[f[i].diff(x[j]) for j in range(2)] for i in range(2)]) xtable = [[1.5, 3.5]] Ertable = [1] NMax = 100 i = 0 eps = 0.000001 Delta = -(K.inv()*f) while Ertable[i] > eps and i <= NMax: xtable.append(np.add(list(Delta.subs({x1:xtable[-1][0], x2:xtable[-1][1]})), xtable[-1])) er = np.linalg.norm(np.array(Delta.subs({x1:xtable[-2][0], x2:xtable[-2][1]})).astype(np.float64))/np.linalg.norm(np.array(xtable[i+1]).astype(np.float64)) Ertable.append(er) i+=1 T2 = [[i + 1, xtable[i][0], xtable[i][1], xtable[i + 1][0], xtable[i + 1][1], Ertable[i + 1]] for i in range(len(xtable) - 1)] T2 = pd.DataFrame(T2, columns=["Iteration", "x1_in", "x2_in", "x1_out", "x2_out", "er"]) display(T2)