Newton Method on matrix nonlinear systems

Introduction

In a previous article of the algorithmic chef (a.k.a me !), I presented Newton’s method. In this one, I’ll demonstrate it’s applicability on nonlinear systems, via matrix nonlinear systems.

System of non-linear equations

Newton’s method could also be applied to system of nonlinear equations. Now, consider

(1) $\begin{equation*}\begin{split}f_1(x_1,x_2 \ldots x_n) & = 0\\f_2(x_1,x_2 \ldots x_n) & = 0\\f_3(x_1,x_2 \ldots x_n) & = 0\\\vdots \\f_n(x_1,x_2 \ldots x_n) & = 0\end{split}\end{equation*}$

Alternatively, and in a more compact form, we can say

(2) $\begin{equation*} \mathbf{f}(\mathbf{x}) = \mathbf{0} \end{equation*}$

where

(3) $\begin{equation*} \mathbf{x}= \left( \begin{array}{c} x_1\\ x_2\\ x_3\\ \vdots\\ x_n \end{array} \right) \end{equation*}$

(4) $\begin{equation*} \mathbf{f}(\mathbf{x})= \left( \begin{array}{c} f_1(x_1,x_2 \ldots x_n)\\ f_2(x_1,x_2 \ldots x_n)\\ f_3(x_1,x_2 \ldots x_n)\\ \vdots\\ f_n(x_1,x_2 \ldots x_n) \end{array} \right) \end{equation*}$

and $\mathbf{0}$ is an all-zero vector.

Newton’s method in Action

Newton method in this case would iterate as follows

(5) $\begin{equation*} \mathbf{x}_k = \mathbf{x}_{k-1} - \mathbf{D}^{-1}(x_{k-1}) \mathbf{f}(\mathbf{x}_{k-1} ) \end{equation*}$

where

(6) $\begin{equation*} \mathbf{D}(\mathbf{x}) = \begin{bmatrix} \frac{\partial f_1}{\partial x_1}(\mathbf{x}) & \frac{\partial f_1}{\partial x_2}(\mathbf{x}) & \frac{\partial f_1}{\partial x_3}(\mathbf{x}) & \dots & \frac{\partial f_1}{\partial x_n}(\mathbf{x}) \\ \frac{\partial f_2}{\partial x_1}(\mathbf{x}) & \frac{\partial f_2}{\partial x_2}(\mathbf{x}) & \frac{\partial f_2}{\partial x_3}(\mathbf{x}) & \dots & \frac{\partial f_2}{\partial x_n}(\mathbf{x}) \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \frac{\partial f_n}{\partial x_1}(\mathbf{x}) & \frac{\partial f_n}{\partial x_2}(\mathbf{x}) & \frac{\partial f_n}{\partial x_3}(\mathbf{x}) & \dots & \frac{\partial f_n}{\partial x_n}(\mathbf{x}) \ \end{bmatrix} \end{equation*}$

Working a step-by-step Example

To avoid confusion of notation, and align the problem with the statements above we note

(7) $\begin{equation*} \begin{split} x_1 &= x \ x_2 &= y \ x_3 &= z \end{split} \end{equation*}$

We have the following system (equation (1) for $n = 3$ )

(8) $\begin{equation*} \begin{split} f_1(x_1,x_2,x_3) & = 16x_1^4 + 16x_2^4 + x_3^4 - 16 = 0\\ f_2(x_1,x_2,x_3) & = x_1^2 + x_2^2 + x_3^2 - 3= 0\\ f_3(x_1,x_2,x_3) & = x_1^3 - x_2 = 0 \end{split} \end{equation*}$

Let’s compute the gradient matrix $\mathbf{D}(\mathbf{x})$ , which in this case is $3 \times 3$

(9) $\begin{equation*} \mathbf{D}(\mathbf{x}) = \begin{bmatrix} \frac{\partial f_1}{\partial x_1}(\mathbf{x}) & \frac{\partial f_1}{\partial x_2}(\mathbf{x}) & \frac{\partial f_1}{\partial x_3}(\mathbf{x}) \\ \frac{\partial f_2}{\partial x_1}(\mathbf{x}) & \frac{\partial f_2}{\partial x_2}(\mathbf{x}) & \frac{\partial f_2}{\partial x_3}(\mathbf{x}) \\ \frac{\partial f_3}{\partial x_1}(\mathbf{x}) & \frac{\partial f_3}{\partial x_2}(\mathbf{x}) & \frac{\partial f_3}{\partial x_3}(\mathbf{x}) \end{bmatrix} \end{equation*}$

The $3 \times 3 = 9$ above derivatives are easily computed as follows

(10) $\begin{equation*} \begin{split} \frac{\partial f_1}{\partial x_1}(\mathbf{x}) &= \frac{\partial}{\partial x_1} \big( 16x_1^4 + 16x_2^4 + x_3^4 - 16 \big) = 64x_1^3 \\ \frac{\partial f_1}{\partial x_2}(\mathbf{x}) &= \frac{\partial}{\partial x_2} \big( 16x_1^4 + 16x_2^4 + x_3^4 - 16 \big) = 64x_2^3 \\ \frac{\partial f_1}{\partial x_3}(\mathbf{x}) &= \frac{\partial}{\partial x_3} \big( 16x_1^4 + 16x_2^4 + x_3^4 - 16 \big) = 4x_3^3 \\ \frac{\partial f_2}{\partial x_1}(\mathbf{x}) &= \frac{\partial}{\partial x_1} \big( x_1^2 + x_2^2 + x_3^2 -3 \big) = 2x_1 \\ \frac{\partial f_2}{\partial x_2}(\mathbf{x}) &= \frac{\partial}{\partial x_2} \big( x_1^2 + x_2^2 + x_3^2 -3 \big) = 2x_2 \\ \frac{\partial f_2}{\partial x_3}(\mathbf{x}) &= \frac{\partial}{\partial x_3} \big( x_1^2 + x_2^2 + x_3^2 -3 \big) = 2x_3 \\ \frac{\partial f_3}{\partial x_1}(\mathbf{x}) &= \frac{\partial}{\partial x_1} \big( x_1^3 - x_2 \big) = 3x_1^2 \\ \frac{\partial f_3}{\partial x_2}(\mathbf{x}) &= \frac{\partial}{\partial x_2} \big( x_1^3 - x_2 \big) = -1 \\ \frac{\partial f_3}{\partial x_3}(\mathbf{x}) &= \frac{\partial}{\partial x_3} \big( x_1^3 - x_2 \big) = 0 \ \end{split} \end{equation*}$

Plugging (10) in (9) we get

(11) $\begin{equation*} \mathbf{D}(\mathbf{x}) = \begin{bmatrix} 64x_1^3 & 64x_2^3 & 4x_3^3 \\ 2x_1 & 2x_2 & 2x_3\\ 3x_1^2 & -1 & 0 \end{bmatrix} \end{equation*}$

Therefore, the Newton iterations are done as follows (substitute (11) in (5) and use $x,y,z$ instead of $x_1,x_2,x_3$ )

(12) $\begin{equation*} \begin{bmatrix} x_{k} \\ y_{k}\\ z_{k} \end{bmatrix} = \begin{bmatrix} x_{k-1} \\ y_{k-1}\\ z_{k-1} \end{bmatrix} - \begin{bmatrix} 64x_{k-1}^3 & 64y_{k-1}^3 & 4z_{k-1}^3 \\ 2x_{k-1} & 2y_{k-1}& 2z_{k-1}\\ 3x_{k-1}^2 & -1 & 0 \end{bmatrix}^{-1} \begin{bmatrix} f_1(x_{k-1} )\\ f_2(y_{k-1})\\ f_3(z_{k-1}) \end{bmatrix} \end{equation*}$

where the initial guess is (this is given in the requirements)

(13) $\begin{equation*} \begin{bmatrix} x_{0} \\ y_{0}\\ z_{0} \end{bmatrix} = \begin{bmatrix} 1 \\ 1\\ 1 \end{bmatrix} \end{equation*}$

It is now easy to iterate, let’s compute the values at first iteration for example

(14) $\begin{equation*} \begin{bmatrix} x_{1} \\ y_{1}\\ z_{1} \end{bmatrix} = \begin{bmatrix} x_{0} \\ y_{0}\\ z_{0} \end{bmatrix} - \begin{bmatrix} 64x_{0}^3 & 64y_{0}^3 & 4z_{0}^3 \\ 2x_{0} & 2y_{0}& 2z_{0}\\ 3x_{0}^2 & -1 & 0 \end{bmatrix}^{-1} \begin{bmatrix} f( x_{0} )\\ f( y_{0} )\\ f( z_{0}) \end{bmatrix} \end{equation*}$

Using (13), we have

(15) $\begin{equation*} \begin{bmatrix} x_{1} \\ y_{1}\\ z_{1} \end{bmatrix} = \begin{bmatrix} 1 \\ 1\\ 1 \end{bmatrix} - \begin{bmatrix} 64 & 64 & 4 \\ 2& 2& 2\\ 3 & -1 & 0 \end{bmatrix}^{-1} \begin{bmatrix} 17 \\ 0\\ 0 \end{bmatrix} = \begin{bmatrix} 0.9292 \\ 0.7875 \\ 1.2833 \end{bmatrix} \end{equation*}$

and so on.. Using Newton iteration, refer to the figure below

Convergence behavior of Newton method for the given non linear system

Furthermore, we see in the above, which is obtained by running the MATLAB code, which is presented in the next section, the convergence of all methods. Also, note here that subroutine is attained by using MATLAB’s fsolve function to compare it with the Newton method

MATLAB step-by-step explanations

Now, it’s time to implement the matrix form of Newton method, enabling us to solve non-linear. Firstly, let’s define the system we have in equation (8) in a separate MATLAB function calling it root3d.m

function F = root2d(x,y,z)

F(1) = 16*x(1)^4 + 16*x(2)^4+ x(3)^4 - 16;
F(2) = x(1)^2 + x(2)^2 + x(3)^2 - 3;
F(3) = x(1)^3 - x(2);

Following the above, we will be using the above function in our main function to get the exact roots of the system (referred to as subroutine) and comparing Now, in the main function, we can easily implement our Newton method as such

Niter = 1e1;
% values of x,y,z at initial guess
v(:,1) = [1;1;1]; 
%functions 
f =@(x1,x2,x3)[16*x1^4 + 16*x2^4+ x3^4 - 16;x1^2 + x2^2 + x3^2 - 3;x1^3 - x2];

% Newton iterations
for n = 2:Niter
   x = v(1,n-1);
   y = v(2,n-1);
   z = v(3,n-1);
   
   D = [64*x^3  64*y^3   4*z^3;...
   2*x  2*y  2*z;...
   3*x^2  -1  0];
v(:,n) = v(:,n-1) - inv(D)*f(x,y,z);
end

In order to compare our iterative approach, we shall compare it to the solution given to us by MATLAB’s fsolve function as

%Matlab subroutine
fun = @root3d;
x0 = [1,1,1];
x_subroutine = fsolve(fun,x0);

where in the above we initialized fsolve to all-ones. Next, we plot

figure
plot(v(1,:),'r','Linewidth',1)
hold on
plot(v(2,:),'m','Linewidth',1)
plot(v(3,:),'k','Linewidth',1)
plot(x_subroutine(1)*ones(Niter,1),'--*r','Linewidth',2)
plot(x_subroutine(2)*ones(Niter,1),'--*m','Linewidth',2)
plot(x_subroutine(3)*ones(Niter,1),'--*k','Linewidth',2)
xlabel('Iteration number ','Interpreter','latex')
ylabel('x','Interpreter','latex')
legend('','','',' subroutine',' subroutine',' subroutine','Interpreter','latex')
grid on
grid minor

Finally, our full main function looks like this


Niter = 1e1;
% values of x,y,z at initial guess
v(:,1) = [1;1;1]; 
%functions 
f =@(x1,x2,x3)[16*x1^4 + 16*x2^4+ x3^4 - 16;x1^2 + x2^2 + x3^2 - 3;x1^3 - x2];

% Newton iterations
for n = 2:Niter
   x = v(1,n-1);
   y = v(2,n-1);
   z = v(3,n-1);
   
   D = [64*x^3  64*y^3   4*z^3;...
   2*x  2*y  2*z;...
   3*x^2  -1  0];
v(:,n) = v(:,n-1) - inv(D)*f(x,y,z);
end
%Matlab subroutine
fun = @root3d;
x0 = [1,1,1];
x_subroutine = fsolve(fun,x0);
figure
plot(v(1,:),'r','Linewidth',1)
hold on
plot(v(2,:),'m','Linewidth',1)
plot(v(3,:),'k','Linewidth',1)
plot(x_subroutine(1)*ones(Niter,1),'--*r','Linewidth',2)
plot(x_subroutine(2)*ones(Niter,1),'--*m','Linewidth',2)
plot(x_subroutine(3)*ones(Niter,1),'--*k','Linewidth',2)
xlabel('Iteration number ','Interpreter','latex')
ylabel('x','Interpreter','latex')
legend('','','',' subroutine',' subroutine',' subroutine','Interpreter','latex')
grid on
grid minor

The above plots the iterations given by the Newton method for the given non-linear system in matrix form, so as to get the figure below

matrix newton method — Convergence behavior of Newton method for the given non linear system

Finally, Subscribe to my channel to support us ! We are very close to 100K subscribers <3 w/ love. The algorithmic chef – Ahmad Bazzi. Click here to subscribe. Check my other articles here.

M	T	W	T	F	S	S
1	2	3	4	5	6	7
8	9	10	11	12	13	14
15	16	17	18	19	20	21
22	23	24	25	26	27	28
29	30

Newton Method on matrix nonlinear systems

Introduction

System of non-linear equations

Newton’s method in Action

Working a step-by-step Example

MATLAB step-by-step explanations

You may also like...

Leave a Reply Cancel reply

Latest Vid: Tensor RT 8.2

Lecture Notes

Recent Comments

Archives

Categories

Meta

Pages

Ahmad Bazzi

Introduction

System of non-linear equations

Newton’s method in Action

Working a step-by-step Example

MATLAB step-by-step explanations

You may also like...

NVIDIA TensorRT

Stocks Programming

Cone Programming on CVXOPT in Python

Leave a Reply Cancel reply

Latest Vid: Tensor RT 8.2

Lecture Notes

Recent Comments

Archives

Categories

Meta

Pages

Tags

Ahmad Bazzi