Bivariate Newton-Raphson solution for non-linear equation system

Systems of non-linear equations can be solved using iterative optimization. The example presented hereafter employs the Newton-Raphson method to find the intersection of two bi-variate polynomials $g_1(x,y)$ and $g_2(x,y)$.
The provided example is authored by Christopher Hahne and inspired by a lab note from Colorado State University.

last update: 01/10/2020

Hahne, PhD

Newton's method suits to solve a non-linear equation system via first-order Taylor expansion of high-order polynomials $g_1(x,y)$ and $g_2(x,y)$ given by

$$g_1(x,y)=x^2+y^2-1,\quad g_2(x,y)=x^4-y^4+xy$$

This is converted to vector-valued functions $\mathbf f\left(x, y\right)$ containing the above equations via

$$ \mathbf f\left(x, y\right) = \begin{bmatrix} g_1(x,y)\\ g_2(x,y)\end{bmatrix} = \begin{bmatrix} x^2+y^2-1 \\ x^4-y^4+xy \end{bmatrix} $$

which is implemented as

In the bi-variate case, we algebraically formulate the Fréchet-derivative as the Jacobian $\mathbf J(x, y)$ which can be obtained analytically as a square matrix given by

$$ \mathbf{J}(x, y) = \begin{bmatrix} \dfrac{\partial g_1(x,y)}{\partial x} & \dfrac{\partial g_1(x,y)}{\partial y}\\[1em] \dfrac{\partial g_2(x,y)}{\partial x} & \dfrac{\partial g_2(x,y)}{\partial y} \end{bmatrix} = \begin{bmatrix} 2x & 2y \\ 4x^3+y & -4y^3+x \end{bmatrix} $$

For compactness, we let $\mathbf{x}_k \in \mathbb{R}^m$ with $\mathbf{x}_k=\left(x_k^{(1)}, x_k^{(2)}, ... ,x_k^{(m)}\right)^\intercal$ and $k\ge0$ is the iteration index so that a multivariate Newton method writes

$$ \mathbf{x}_{k+1}=\mathbf{x}_k - \big[\mathbf{J}(\mathbf{x}_k)\big]^{-1}\mathbf{f}(\mathbf{x}_k) $$

for which we may plug in $\mathbf x_k=(x_k,y_k)^{\intercal}$ from the above terms giving

$$ \begin{bmatrix} x_{k+1} \\ y_{k+1}\end{bmatrix}= \begin{bmatrix} x_k \\ y_k\end{bmatrix} - \begin{bmatrix} 2x_k & 2y_k \\ 4x_k^3+y_k & -4y_k^3+x_k \end{bmatrix}^{-1} \begin{bmatrix} x_k^2+y_k^2-1 \\ x_k^4-y_k^4+x_k y_k \end{bmatrix} $$

An implementation of the multivariate Newton method is provided hereafter.