How to solve system of (non-linear) equations using Jacobian and Newton's Method in Matlab

Question

I am trying to build a function that can solve a system of n-1 (non-linear) equations with n unknowns in Matlab making use of Newton's method.

I did some research online and came to the following procedure:

Use matlab's 'null' function applied to an array of first partial derivatives (Jacobian) to obtain an m-element vector, (call this vector v,) which is to be orthogonal to all the m-1 vectors of the partials. After proceeding along this vector a distance h0, use the 'null' function again, this time to find the m-1 dimensional subspace orthogonal to v from that advanced point. Projecting all the original coordinates to coordinates in this subspace you can perform an m-1-dimensional Newton-Raphson procedure which should converge very rapidly provided your step h0 is not too large. Then convert this back to the original coordinates.

I have the following system of equations:

fun = @(x,y,z)[x.^3+y.^2+z.^2,x.^2-y.^3+sin(z)]

and it's Jacobian:

Jac = @(x,y,z)[3*x^2,    2*y,    2*z; 2*x, -3*y^2, cos(z)]

And an initial guess: x0 = [1,1,1]

In my function I convert the vector to numbers, so I can work with them: x = num2cell(x0)

Then I want to do N iterations:

for numIterations = 1:N

%calculate the function values at the current iteration
F = fun(x{:})

%calculate the jacobian matrix
J = Jac(x{:})

%find direction in which the curve extends:
v = null(J)

xnew = xold + h0 * v
xnew = xnew'

subspace = null(xnew)

nextX = multinewton(f,df,xnew)

Now I am stuck in the part where I use Newton's method to calculate the next coordinates. How do I do this?

The Newton method I want to use to calculate the next coordinates, is the following method:

function [zero,res,niter]=newton(f,df,x0,tol,nmax,varargin)
%NEWTON Find function zeros.
%   ZERO=NEWTON(FUN,DFUN,X0,TOL,NMAX) tries to find the zero ZERO of the 
%   continuous and differentiable function FUN nearest to X0 using the    Newton 
%   method. FUN and its derivative DFUN accept real scalar input x and returns 
%   a real scalar value. If the search fails an error message is displayed.
%   FUN and DFUN can also be inline objects.
%
%   [ZERO,RES,NITER]= NEWTON(FUN,...) returns the value of the residual in ZERO
%   and the iteration number at which ZERO was computed.

x = x0; 
fx = feval(f,x,varargin{:}); 
dfx = feval(df,x,varargin{:});
niter = 0; 
diff = tol+1;
resultVector = [x];
while diff >= tol && niter <= nmax
disp(x)
niter = niter + 1;
diff = - dfx\fx;
x = x + diff;     
diff = abs(diff);
fx = feval(f,x,varargin{:});        
dfx = feval(df,x,varargin{:});
resultVector = [resultVector x];
end
if niter > nmax
fprintf(['newton stopped without converging to the desired tolerance',...
'because the maximum number of iterations was reached\n']);
end
res = fx
zero = x;
disp('Number of iterations is:');
disp(niter);

plot(resultVector);
return

I think I do not completely understand the part 'projecting all the original coordinates to coordinates in this subspace you can perform an m-1-dimensional Newton-Raphson procedure'. What is meant by this? What do I give as input to the newton method?

Many thanks in advance!

Answer 1

I'll answer the question of how one can solve a system of n-1 equations with n unknowns in Matlab by adapting Newton's method. My adaptation is not the one you found through your research -- it's simpler.

The idea of Newton's method is that we linearize the system around some guess point and solve the resulting linear system. The solution becomes our next guess.

If your system is underdetermined (has fewer equations than unknowns), the linear system will be underdetermined too. Matlab's \ operator handles such systems quite well: returns the solution of smallest norm. So, why not just use this, if it works? No changes of code are needed, compared to the case of n equations with n unknowns.

Here's my implementation: calling newtonund(1,1,1) find a root in 20 steps, namely -1.226458955, 1.357877375, -0.031788652. I changed your system of equations slightly, replacing sin(z) by cos(z). It works in the original form too, but the result is (predictably) 0, 0, 0 which is not as interesting.

function newtonund(x,y,z)
tol = 1e-6;
diff = tol+1;
n = 0
nmax = 1000;
disp('  n      x(n)            y(n)             z(n)          |f(x)|');
X = [x;y;z];
Y = F(X);
while diff>=tol & n<=nmax
    changeX = -dF(X)\Y;
    X = X + changeX;
    Y = F(X);
    diff = norm(changeX);
    n = n+1;
    fprintf('%3d %15.9f %15.9f %15.9f %10.5g \n', n, X(1), X(2), X(3), norm(Y));
end
if n>nmax
    disp('Failed to converge');
else
    disp('Root found');
end
end

function Y=F(X)
x=X(1); y=X(2); z=X(3);
Y=[x^3+y^2+z^2 ; x^2-y^3+cos(z)];
end

function J=dF(X)
x=X(1); y=X(2); z=X(3);
J=[3*x^2, 2*y, 2*z; 2*x, -3*y^2, -sin(z)];
end

The line changeX = -dF(X)\Y expresses Newton's method. Normally dF(X) would be a square matrix; in the underdetermined case it is rectangular but the algorithm works anyway.