Chapter 03.07: Newton-Raphson Method For Solving Simultaneous Nonlinear Equations

Lesson: Newton-Raphson Method For Solving Simultaneous Nonlinear Equations

Learning Objectives

After reading this chapter, you should be able to:

1) derive the Newton-Raphson method formula for simultaneous nonlinear equations,

2) develop the algorithm of the Newton-Raphson method for solving simultaneous nonlinear equations,

3) use the Newton-Raphson method to solve a set of simultaneous nonlinear equations,

4) model a real-life problem that results in a set of simultaneous nonlinear equations.

  

Introduction

Several physical systems result in a mathematical model in terms of simultaneous nonlinear equations. A set of such equations can be written as

$$\begin{split} &f_{1}(x_{1},x_{2},...,x_{n}) = 0\\ &\text{...}\\ &\text{...}\\ &f_{n}(x_{1},x_{2},...,x_{n}) = 0 \end{split}\;\;\;\;\;\;\; (1)$$

The solution to these simultaneous nonlinear equations are values of $x_{1},x_{2},...,x_{n}$ which satisfy all the above $n$ equations. The number of set of solutions to these equations could be none, unique, more than one but finite, or infinite. In this chapter, we use the Newton-Raphson method to solve these equations.

The Newton-Raphson method of solving a single nonlinear equation, $f(x) = 0$, can be derived using first-order Taylor series (first two terms of Taylor series) and are given by

$$f(x_{i + 1}) = f(x_{i}) + f^{\prime}(x_{i})(x_{i + 1} - x_{i})\;\;\;\;\;\;\; (2)$$

where

$$x_{i} = \text{previous estimate of root}$$

$$x_{i + 1} = \text{present estimate of root}$$

Since we are looking for $x_{i + 1}$where $f(x_{i + 1})$ becomes zero, Equation (2) can be re-written as

$$0 = f(x_{i}) + f^{\prime}(x_{i})(x_{i + 1} - x_{i})$$

and then as a recursive formula as

$$x_{i + 1} = x_{i} - \frac{f(x_{i})}{f^{\prime}(x_{i})}\;\;\;\;\;\;\; (3)$$

  

Derivation

Now how do we extend the same to simultaneous nonlinear equations? For sake of simplicity, let us limit the number of nonlinear equations to two as

$$u(x,y) = 0\;\;\;\;\;\;\; (4a)$$

$$v(x,y) = 0\;\;\;\;\;\;\; (4b)$$

The first order Taylor-series for nonlinear equation is (4a) & (4b) are

$$u(x_{i + 1},y_{i + 1}) = u(x_{i},y_{i}) + \left.\frac{\partial u}{\partial x}\right|_{x_{i},y_{i}} (x_{i + 1} - x_{i}) + \left. \frac{\partial u}{\partial y}\right|_{x_{i},y_{i}} (y_{i + 1} - y_{i})\;\;\;\;\;\;\; (5a)$$

$$v(x_{i + 1},y_{i + 1}) = v(x_{i},y_{i}) + \left. \frac{\partial v}{\partial x}\right|_{x_{i},y_{i}} (x_{i + 1} - x_{i}) + \left.\frac{\partial v}{\partial y}\right|_{x_{i},y_{i}}(y_{i + 1} - y_{i})\;\;\;\;\;\;\; (5b)$$

We are looking for $(x_{i + 1},y_{i + 1})$ where $u(x_{i + 1},y_{i + 1})$ and $v(x_{i + 1},y_{i + 1})$ are zero. Hence

$$0 = u(x_{i},y_{i}) + \left.\frac{\partial u}{\partial x}\right|_{x_{i},y_{i}} (x_{i + 1} - x_{i}) + \left.\frac{\partial u}{\partial y}\right|_{x_{i},y_{i}}(y_{i + 1} - y_{i}) \;\;\;\;\;\;\;(6a)$$

$$0 = v(x_{i},y_{i}) + \left.\frac{\partial v}{\partial x}\right|_{x_{i},y_{i}} (x_{i + 1} - x_{i}) + \left.\frac{\partial v}{\partial y}\right|_{x_{i},y_{i}} (y_{i + 1} - y_{i})\;\;\;\;\;\;\; (6b)$$

Writing

$$x_{i + 1} - x_{i} = \Delta x\;\;\;\;\;\;\; (7a)$$

$$y_{i + 1} - y_{i} = \Delta y\;\;\;\;\;\;\; (7b)$$

we get

$$0 = u(x_{i},y_{i}) + \left.\frac{\partial u}{\partial x}\right|_{x_{i},y_{i}} \Delta x + \left.\frac{\partial u}{\partial y}\right|_{x_{i},y_{i}} \Delta y\;\;\;\;\;\;\; (8a)$$

$$0 = v(x_{i},y_{i}) + \left.\frac{\partial v}{\partial x}\right|_{x_{i},y_{i}} \Delta x + \left.\frac{\partial v}{\partial y}\right|_{x_{i},y_{i}} \Delta y\;\;\;\;\;\;\; (8b)$$

Rewriting Equations (8a) and (8b)

$$\left.\frac{\partial u}{\partial x}\right| x_{i},y_{i} \Delta x + \left.\frac{\partial u}{\partial y}\right|_{x_{i},y_{i}} \Delta y = - u(x_{i},y_{i})\;\;\;\;\;\;\; (9a)$$

$$\left.\frac{\partial v}{\partial x}\right| x_{i},y_{i} \Delta x + \left.\frac{\partial v}{\partial y}\right|_{x_{i},y_{i}} \Delta y = - v(x_{i},y_{i})\;\;\;\;\;\;\; (9b)$$

and then in the matrix form

$$\begin{bmatrix} \displaystyle\left.\frac{\partial u}{\partial x}\right|_{x_i,y_i} & \displaystyle\left.\frac{\partial u}{\partial y}\right|_{x_i,y_i}\\ \displaystyle\left.\frac{\partial v}{\partial x}\right|_{x_i,y_i} & \displaystyle\left.\frac{\partial v}{\partial y}\right|_{x_i,y_i} \end{bmatrix} \begin{bmatrix} \Delta x \\ \Delta y \\ \end{bmatrix} = \begin{bmatrix} - u(x_{i},y_{i}) \\ - v(x_{i},y_{i}) \\ \end{bmatrix}\;\;\;\;\;\;\; (10)$$

Solving Equation (10) would give us $\Delta x$ and $\Delta y$. Since the previous estimate of the root is $(x_{i},y_{i})$, one can find from Equation (7)

$$x_{i + 1} = x_{i} + \Delta x\;\;\;\;\;\;\; (11a)$$

$$y_{i + 1} = y_{i} + \Delta y\;\;\;\;\;\;\; (11b)$$

This process is repeated till one obtains the root of the equation within a prespecified tolerance,$\varepsilon_{s}$ such that

$$\left| \varepsilon_{a} \right|_{x} = \left| \frac{x_{i + 1} - x_{i}}{x_{i + 1}} \right| \times 100 < \varepsilon_{s}$$

$$\left| \varepsilon_{a} \right|_{y} = \left| \frac{y_{i + 1} - y_{i}}{y_{i + 1}} \right| \times 100 < \varepsilon_{s}$$

  

Example 1

Find the roots of the simultaneous nonlinear equations

$$\begin{matrix} x^{2} + y^{2} = 50 \\ x - y = 10 \\ \end{matrix}$$

Use an initial guess of $(x,y) = (2, - 4)$. Conduct two iterations.

Solution

First put the equations in the form

$$u(x,y) = 0$$

$$v(x,y) = 0$$

to give

$$u(x,y) = x^{2} + y^{2} - 50 = 0$$ $$v(x,y) = x - y - 10 = 0$$

Now

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

$$\frac{\partial u}{\partial y} = 2y$$

$$\frac{\partial v}{\partial x} = 1$$

$$\frac{\partial v}{\partial y} = - 1$$

Hence from Equation (10)

$$\begin{bmatrix} \displaystyle\left.\frac{\partial u}{\partial x}\right|_{x_i,y_i} & \displaystyle\left.\frac{\partial u}{\partial y}\right|_{x_i,y_i}\\ \displaystyle\left.\frac{\partial v}{\partial x}\right|_{x_i,y_i} & \displaystyle\left.\frac{\partial v}{\partial y}\right|_{x_i,y_i} \end{bmatrix} \begin{bmatrix} \Delta x \\ \vdots \\ \Delta y \\ \end{bmatrix} = \begin{bmatrix} - u(x_{i},y_{i}) \\ \vdots \\ - v(x_{i},y_{i}) \\ \end{bmatrix}$$

$$\begin{bmatrix} 2x_{i} & 2y_{i} \\ 1 & - 1 \\ \end{bmatrix}\begin{bmatrix} \Delta x \\ \Delta y \\ \end{bmatrix} = \begin{bmatrix} - {x_{i}}^{2} - y_{i}^{2} + 50 \\ - x_{i} + y_{i} + 10 \\ \end{bmatrix}$$

  

Iteration 1

The initial guess

$$(x_{i},y_{i}) = (2, - 4)$$

Hence

$$\begin{bmatrix} 2(2) & 2( - 4) \\ 1 & - 1 \\ \end{bmatrix}\begin{bmatrix} \Delta x \\ \Delta y \\ \end{bmatrix} = \begin{bmatrix} - (2)^{2} - ( - 4)^{2} + 50 \\ - (2) + ( - 4) + 10 \\ \end{bmatrix}$$

$$\begin{bmatrix} 4 & - 8 \\ 1 & - 1 \\ \end{bmatrix}\begin{bmatrix} \Delta x \\ \Delta y \\ \end{bmatrix} = \begin{bmatrix} 30 \\ 4 \\ \end{bmatrix}$$

Solving the equations by any method of your choice, we get

$$\begin{matrix} \Delta x = 0.5000 \\ \Delta y = - 3.500 \\ \end{matrix}$$

Since

$$\begin{matrix} \Delta x = x_{2} - x_{1} = 0.5000 \\ \Delta y = y_{2} - y_{1} = - 3.500 \\ \end{matrix}$$

we get

$$\begin{split} x_{2}&= x_{1} + 0.5000 \\ &= 2 + 0.5000 \\ &= 2.5000 \end{split}$$

$$\begin{split} y_{2}&= y_{1} + ( - 3.500) \\ &= - 4 + ( - 3.500) \\ &= - 7.500 \end{split}$$

The absolute relative approximate errors at the end of the first iteration are

$$\begin{split} \left| \varepsilon_{a} \right|_{x}&= \left| \frac{x_{2} - x_{1}}{x_{2}} \right| \times 100 \\ &= \left| \frac{2.500 - 2.000}{2.500} \right| \times 100 \\ &= 20.00\% \end{split}$$

$$\begin{split} \left| \varepsilon_{a} \right|_{y} &= \left| \frac{y_{2} - y_{1}}{y_{2}} \right| \times 100\\ &= \left| \frac{- 7.500 - ( - 4.000)}{- 7.500} \right| \times 100 \end{split}$$

  

Iteration 2

The estimate of the root at the end of iteration#1 is

$$(x_{2},y_{2}) = (2.500 - 7.500)$$

Hence

$$\begin{bmatrix} 2(2.500) & 2( - 7.500) \\ 1 & - 1 \\ \end{bmatrix}\begin{bmatrix} \Delta x \\ \Delta y \\ \end{bmatrix} = \begin{bmatrix} - (2.500)^{2} - ( - 7.500)^{2} + 50 \\ - (2.500) + ( - 7.500) + 10 \\ \end{bmatrix}$$

$$\begin{bmatrix} 5 & - 15 \\ 1 & - 1 \\ \end{bmatrix}\begin{bmatrix} \Delta x \\ \Delta y \\ \end{bmatrix} = \begin{bmatrix} - 12.50 \\ 0 \\ \end{bmatrix}$$

Solving the above equations by any method of your choice gives

$$\Delta x = 1.250$$

$$\Delta y = 1.250$$

Since

$$\Delta x = x_{3} - x_{2} = 1.250$$

$$\Delta y = y_3 - y_2 = 1.250$$

giving

$$\begin{split} x_3 &= x_2 +1.250\\ &=2.500 + 1.250\\ &=3.750 \end{split}$$

$$\begin{split} y_{3} &= y_{2} + 1.250\\ &=-7.500 + 1.250\\ &=-6.250 \end{split}$$

The absolute relative approximate error at the end of the second iteration is

$$\begin{split} \left| \epsilon_{a} \right|_{x} &= \left| \frac{x_{3} - x_{2}}{x_{3}} \right| \times 100\\ &= \left| \frac{3.750 - 2.5}{3.750} \right| \times 100\\ &= 33.33\% \end{split}$$

$$\begin{split} \left| \epsilon_{a} \right|_{y} &= \left| \frac{y_{3} - y_{2}}{y_{3}} \right| \times 100\\ &= \left| \frac{-6.250 - (-7.500)}{-6.250}\right| \times 100\\ &= 20.00\% \end{split}$$

Although not asked in the example problem statement, the estimated values of the root and the absolute relative approximate errors are given below in Table 1.

Table 1: Estimate of the root and absolute relative approximate error.

Iteration number, $i$	$x_{i}$	$y_{i}$	$\left| \epsilon_{a} \right|_{x}\%$	$\left| \epsilon_{a} \right|_{y}\%$
$1$	$2.500$	$-7.500$	$20.00$	$46.67$
$2$	$3.750$	$-6.250$	$33.33$	$20.00$
$3$	$4.375$	$-5.625$	$14.29$	$11.11$
$4$	$4.688$	$-5.312$	$6.667$	$5.882$
$5$	$4.844$	$-5.156$	$3.226$	$3.030$
$6$	$4.922$	$-5.078$	$1.587$	$1.538$

The exact solution to which the above scheme is converging to is $(x,y) = (5, - 5)$

  

Example 2

A $5 \text{ m}$ long gutter is made from a flat sheet of aluminum which is $5\ \text{m} \times 0.21\ \text{m}$. The shape of the gutter cross-section is shown in Figure 1 and is made by bending the sheet at two locations at an angle $\theta$ (Figure 2). What are the values of $s$ and $\theta$, that will maximize the volume capacity of the gutter so that it drains water quickly during a heavy rainfall?

Figure 1 Sheet metal bent to form a gutter

Figure 2: Labeling the parameters and points of the gutter

Solution

One needs to maximize the cross-sectional area of the gutter. The cross-sectional area $G$ of the gutter is given by

$$\begin{split} G &= \text{Area of trapezoid}\ \text{ABCD}\\ &= \frac{1}{2}(\text{AB} + \text{CD})(\text{FC})\;\;\;\;\;\;\; (E2.1)\end{split}$$

where

$$\text{AB}\ = \ L - 2\ s\;\;\;\;\;\;\; (E2.2)$$

$$\begin{split} CD &= AB + EA + BF\\ &= (L - 2s) + s\cos(\theta) + s\cos(\theta)\\ &= (L - 2s)\ + \ 2s\cos(\theta)\;\;\;\;\;\;\; (E2.3) \end{split}$$

$$\begin{split} FC &= BC \sin(\theta)\\ &= s\sin(\theta)\;\;\;\;\;\;\; (E2.4) \end{split}$$

hence

$$G = \frac{1}{2}(AB\ + CD) (DE)$$

$$\begin{split} G(s,\theta) &= \frac{1}{2}\left\lbrack (L - 2s) + (L - 2s) + 2s\cos(\theta) \right\rbrack s \sin(\theta)\\ &= \left\lbrack L - 2s + s\cos(\theta) \right\rbrack\ s\sin(\theta)\;\;\;\;\;\;\; (E2.5) \end{split}$$

For example, for $\theta = 0{^\circ}$, the cross-sectional area of the gutter becomes

$$G = 0$$

as it represents a flat sheet, for $\theta = 180{^\circ}$, the cross-sectional area of the gutter becomes

$$G = 0$$

as it represents an overlapped sheet, and for $\theta = 90{^\circ}$, it represents a rectangular cross-sectional area with

$$G = (L - 2s) s$$

For the given $L = 0.21\ \text{m,}$

$$G(s,\theta) = (0.21 - 2s + s\cos\theta)s\sin\theta$$

Then

$$\frac{\partial G}{\partial s} = - s^{2}\sin^{2}\theta + (0.21 - 2s + s\cos\theta)s\cos\theta$$

$$\frac{\partial G}{\partial\theta} = ( - 2 + \cos\theta)s\sin\theta + (0.21 + 2s + s\cos\theta)\sin\theta$$

By solving the equations

$$f_{1}(s,\theta) = - s\sin^{2}\theta + (0.21 - 2s + s\cos\theta)s\cos\theta = 0$$

$$f_{2}(s,\theta) = ( - 2 + \cos\theta)s\sin\theta + (0.21 - 2s + s\cos\theta)\sin\theta = 0$$

we can find the local minimas and maximas of $G(s,\theta)$. One of the local maximas may also be the absolute maximum. The values of $s$ and $\theta$ that correspond to the maximum of $G(s,\theta)$ are what we are looking for. Can you solve the equations to find the corresponding values of $s$ and $\theta$ for maximum value of $G$? Intuitively, what do you think would be these values of $s$ and $\theta$?

  

Appendix A: General matrix form of solving simultaneous nonlinear equations.

The general system of equations is given by

$$\begin{split} &f_{1}(x_{1},x_{2},...,x_{n}) = 0\\ &f_{2}(x_{1},x_{2},...,x_{n}) = 0\\ &\text{...}\\ &\text{...}\\ &\text{...}\\ &f_n(x_1,x_2,...,x_n) = 0\;\;\;\;\;\;\; (A.1)\end{split}$$

We can rewrite in a form as

$$F(\overset{\land}{x}) = \widehat{0}\;\;\;\;\;\;\; (A.2)$$

where

$$F(\overset{\land}{x}) = \begin{bmatrix} f_{1}(\overset{\land}{x}) \\ f_{2}(\overset{\land}{x}) \\ . \\ . \\ f_{n}(\overset{\land}{x}) \\ \end{bmatrix}\;\;\;\;\;\;\;(A.3)$$

$$\overset{\land}{0} = \begin{bmatrix} 0 \\ 0 \\ . \\ . \\ 0 \\ \end{bmatrix}\;\;\;\;\;\;\; (A.4)$$

$$\overset{\land}{x} = \begin{bmatrix} x_{1} \\ x_{2} \\ . \\ . \\ x_{n} \\ \end{bmatrix}^{T}\;\;\;\;\;\;\; (A.5)$$

The Jacobian of the system of equations by using Newton-Raphson method then is

$$J(x_{1},x_{2},...,x_{n}) = \begin{bmatrix} &\displaystyle\frac{\partial f_{1}}{\partial x_{1}}\ldots\frac{\partial f_{1}}{\partial x_{n}} \\ &.\ldots. \\ &.\ldots. \\ &.\ldots. \\ &\displaystyle\frac{\partial f_{n}}{\partial x_{1}}\ldots\frac{\partial f_{n}}{\partial x_{n}} \\ \end{bmatrix}\;\;\;\;\;\;\; (A.6)$$

Using the Jacobian for the Newton-Raphson method guess

$$\lbrack J\rbrack\Delta\overset{\land}{x} = - F(\overset{\land}{x})\;\;\;\;\;\;\; (A.7)$$

where

$$\Delta\overset{\land}{x} = {\overset{\land}{x}}_{\text{new}} - {\overset{\land}{x}}_{\text{old}}\;\;\;\;\;\;\; (A.8)$$

Hence

$$\Delta\overset{\land}{x} = - \lbrack J\rbrack^{- 1}F(\overset{\land}{x})$$

$${\overset{\land}{x}}_{\text{new}} - {\overset{\land}{x}}_{\text{old}} = - \lbrack J\rbrack^{- 1}F(\overset{\land}{x})$$

$${\overset{\land}{x}}_{\text{new}} = {\overset{\land}{x}}_{\text{old}} - \lbrack J\rbrack^{- 1}F(\overset{\land}{x})\;\;\;\;\;\;\; (A.10)$$

Evaluating the inverse of $\lbrack J\rbrack$ in Equation (A.10) is computationally more intensive than solving Equation (A.7) for $\Delta\overset{\land}{x}$and calculating ${\overset{\land}{x}}_{\text{new}}$ as

$${\overset{\land}{x}}_{\text{new}} = {\overset{\land}{x}}_{\text{old}} + \Delta\overset{\land}{x}\;\;\;\;\;\;\; (A.11)$$

Problem Set

(1).  Solve the following set of nonlinear equations using Newton-Raphson method.

$${xy = 2}$$

$${x^{2} + y = 5}$$ Start with an initial guess of $(3,4)$.

(2).  Write the drawbacks of Newton-Raphson method and how you would address them.

(3).  The following simultaneous nonlinear equations are found while modeling two continuous non-alcoholic stirred tank reactions

$${f_{1} = (1 - R)\left\lbrack \frac{D}{10(1 + \beta_{1})} - \varphi_{1} \right\rbrack\exp\left( \frac{10\varphi_{1}}{1 + 10\varphi_{1}/\gamma} \right) - \varphi_{1}}$$ $${f_{2} = \varphi_{1} - (1 + \beta_{2})\varphi_{2} + (1 - R) \times \left\lbrack \frac{D}{10} - \beta_{1}\varphi_{1} - (1 + \beta_{2})\varphi_{2} \right\rbrack\exp\left( \frac{10\varphi_{2}}{1 + 10\varphi_{2}/\gamma} \right)}$$

Given $\gamma = 1000,\ D = 22,\ \beta_{1} = 2,\ \beta_{2} = 2,\ R = 0.935$, find a suitable solution for $(\varphi_{1},\varphi_{2})$, given that $\varphi_{i} \in \left\lbrack 0,1 \right\rbrack,i = 1,2$

Source: Michael J. Hirsch, Panos M. Pardalos, Mauricio G.C Resende, Solving systems of nonlinear equations with continuous GRASP, Nonlinear Analysis: Real World Applications, Volume 10, Issue 4, August 2009, Pages 2000-2006, ISSN 1468-1218, http://dx.doi.org/10.1016/j.nonrwa.2008.03.006

(4).  The Lorenz equations are given by

$$\frac{dx}{dt} = \sigma(y - x)$$

$$\frac{{dy}}{{dt}} = x(p - z) - y$$

$$\frac{dz}{dt} = xy - \beta{z}$$

Assuming $\sigma = 3,\ p = 26,\ \beta = 1$, the steady state equations are written as

$${0 = 3(y - x)}$$ $${0 = x(26 - z) - y}$$ $${0 = xy - z}$$

Write the Jacobian for the above system of equations.

(5).  Use MATLAB or equivalent program to find the solution to the problem stated in Example 2.