Chapter 03.07: 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.