3.5: System of Nonlinear Equations

A system of nonlinear equations can be solved using a version of Newton’s Method. We illustrate this method for a system of two equations and two unknowns. Suppose that we want to solve

\[f(x, y)=0, \quad g(x, y)=0 \nonumber \]

for the unknowns \(x\) and \(y .\) We want to construct two simultaneous sequences \(x_{0}, x_{1}, x_{2}, \ldots\) and \(y_{0}, y_{1}, y_{2}, \ldots\) that converge to the roots. To construct these sequences, we Taylor series expand \(f\left(x_{n+1}, y_{n+1}\right)\) and \(g\left(x_{n+1}, y_{n+1}\right)\) about the point \(\left(x_{n}, y_{n}\right) .\) Using the partial derivatives \(f_{x}=\partial f / \partial x, f_{y}=\partial f / \partial y\), etc., the twodimensional Taylor series expansions, displaying only the linear terms, are given by

\[\begin{aligned} f\left(x_{n+1}, y_{n+1}\right)=f\left(x_{n}, y_{n}\right)+\left(x_{n+1}-x_{n}\right) f_{x}\left(x_{n}, y_{n}\right) & \\ &+\left(y_{n+1}-y_{n}\right) f_{y}\left(x_{n}, y_{n}\right)+\ldots \end{aligned} \nonumber \]

\(g\left(x_{n+1}, y_{n+1}\right)=g\left(x_{n}, y_{n}\right)+\left(x_{n+1}-x_{n}\right) g_{x}\left(x_{n}, y_{n}\right)\)

\[+\left(y_{n+1}-y_{n}\right) g_{y}\left(x_{n}, y_{n}\right)+\ldots \nonumber \]

To obtain Newton’s method, we take \(f\left(x_{n+1}, y_{n+1}\right)=0, g\left(x_{n+1}, y_{n+1}\right)=0\) and drop higher-order terms above linear. Although one can then find a system of linear equations for \(x_{n+1}\) and \(y_{n+1}\), it is more convenient to define the variables

\[\Delta x_{n}=x_{n+1}-x_{n}, \quad \Delta y_{n}=y_{n+1}-y_{n} . \nonumber \]

The iteration equations will then be given by

\[x_{n+1}=x_{n}+\Delta x_{n}, \quad y_{n+1}=y_{n}+\Delta y_{n} \nonumber \]

and the linear equations to be solved for \(\Delta x_{n}\) and \(\Delta y_{n}\) are given by

\[\left(\begin{array}{ll} f_{x} & f_{y} \\ g_{x} & g_{y} \end{array}\right)\left(\begin{array}{l} \Delta x_{n} \\ \Delta y_{n} \end{array}\right)=\left(\begin{array}{l} -f \\ -g \end{array}\right) \nonumber \]

where \(f, g, f_{x}, f_{y}, g_{x}\), and \(g_{y}\) are all evaluated at the point \(\left(x_{n}, y_{n}\right) .\) The twodimensional case is easily generalized to \(n\) dimensions. The matrix of partial derivatives is called the Jacobian Matrix.

We illustrate Newton’s Method by finding the steady state solution of the Lorenz equations, given by

\[\begin{array}{r} \sigma(y-x)=0 \\ r x-y-x z=0 \\ x y-b z=0 \end{array} \nonumber \]

where \(x, y\), and \(z\) are the unknown variables and \(\sigma, r\), and \(b\) are the known parameters. Here, we have a three-dimensional homogeneous system \(f=0, g=0\), and \(h=0\), with

\[\begin{aligned} &f(x, y, z)=\sigma(y-x) \\ &g(x, y, z)=r x-y-x z \\ &h(x, y, z)=x y-b z . \end{aligned} \nonumber \]

The partial derivatives can be computed to be

\[\begin{array}{lll} f_{x}=-\sigma, & f_{y}=\sigma, & f_{z}=0, \\ g_{x}=r-z, & g_{y}=-1, & g_{z}=-x, \\ h_{x}=y, & h_{y}=x, & h_{z}=-b . \end{array} \nonumber \]

The iteration equation is therefore

\[\left(\begin{array}{ccc} -\sigma & \sigma & 0 \\ r-z_{n} & -1 & -x_{n} \\ y_{n} & x_{n} & -b \end{array}\right)\left(\begin{array}{c} \Delta x_{n} \\ \Delta y_{n} \\ \Delta z_{n} \end{array}\right)=-\left(\begin{array}{c} \sigma\left(y_{n}-x_{n}\right) \\ r x_{n}-y_{n}-x_{n} z_{n} \\ x_{n} y_{n}-b z_{n} \end{array}\right) \nonumber \]

with

\[\begin{aligned} &x_{n+1}=x_{n}+\Delta x_{n} \\ &y_{n+1}=y_{n}+\Delta y_{n} \\ &z_{n+1}=z_{n}+\Delta z_{n} \end{aligned} \nonumber \]

The MATLAB program that solves this system is contained in newton_system.m.