Solving a system on nonlinear equations using Halley's method

[Namir]

I have known for many years that Halley's method for solving single-variable nonlinear functions converged faster than Newton's method.

Early today I found the following 1985 article that discusses Halley's method for solving a system of nonlinear equations. Developing the matrix version of Halley from the regular version is not trivial. The article find is a gem!

I wrote the following Matlab code (should be stored in file sne_halley.m) for the algorithm mentioned in the cited article.

Code:


function [X,Iter] = sne_halley(N, X, Xtol, Ftol, Maxiter, varargin)
% Function that solves nonlinear simultaneous equations

% 
% Example:
% 
% [X, Iter]=sne_halley(4, [6 6 6 6],1e-7,1e-8,1000,@fx1,@fx2,@fx3,@fx4)
% 
% 
% X =
% 
%      2
%      3
%      4
%      5
% 
% 
% Iter =
% 
%      6
% 

X=X'; % transpose row vector X into a column vector
jacobian = zeros(N,N);
hessian=zeros(N,N);
Dx = zeros(N,1);
Cx = zeros(N,1);
Fx = zeros(N,1);

Iter = 0;
bStop = false;

while ~bStop

  Iter = Iter + 1;
  if Iter > Maxiter
    break
  end 
  
  % calculate Jacobian matrix
  for i=1:N
    % First derivatives are calculated using:
    %
    %  dF(X)/dX(i) = (F(X(i)+h) - F(X(i)-h))/(2h)
    %
    % Where h = 0.01*(1 +abs(X(i))
    %
    Fx(i,1) = -feval(varargin{i}, X);
    for j=1:N
     xi = X(i);
     h = 0.01 * (1 + abs(xi));
     X(i) = xi + h;
     fp = feval(varargin{j}, X);
     X(i) = xi - h;
     fm = feval(varargin{j}, X);
     X(i) = xi;     
     jacobian(j,i) = (fp - fm) / 2 / h; % jacobianian matrix elements
    end
    X(i) = xi;
  end
  
  % calculate Hessian matrix
  for i=1:N
    for j=1:N
      if i==j
        % diagonal element. Calculate second derivative with
        % respect to X(i)
        % Second derivatives are calculated using:
        %
        %  dF2(X)/dX(i)2 = (F(X(i)+h) - 2*F(X(i)) + F(X(i)-h))/h^2
        %
        % Where h = 0.01*(1 +abs(X(i))
        %
        xi = X(i);
        h = 0.01 * (1 + abs(xi));
        X(i) = xi + h;
        fp = feval(varargin{i}, X);
        X(i) = xi - h;
        fm = feval(varargin{i}, X);
        X(i) = xi;
        f0 = feval(varargin{i}, X);
        % calculate second derivative of function i for variable X(i)
        hessian(i,j) = (fp - 2 * f0 + fm)/h^2;
      else
        % non-diagonal element. Calculate second derivative with
        % respect to X(i) and X(j)     
        % Second derivatives are calculated using:
        %
        %  dF2(X)/dX(i)dX(j) = ((F(X(i)+hi,X(j)+hj) - F(X(i)-hi,X(j)+hj))/(2hi) 
        %                      - (F(X(i)+hi,X(j)-hj) - F(X(i)-hi,X(j)-hj))/(2hi) 
        %                      / (2hj)
        %
        % Where hi = 0.01*(1 +abs(X(i)) and hj = 0.01*(1+abs(X(j))
        %
        
        xi = X(i);
        xj = X(j);
        hi = 0.01 * (1 + abs(xi));
        hj = 0.01 * (1 + abs(xj));
        X(j) = xj + hj;
        X(i) = xi + hi;
        fp1 = feval(varargin{i}, X);
        X(i) = xi - hi;
        fm1 = feval(varargin{i}, X);
        X(j) = xj - hj;
        X(i) = xi + hi;
        fp2 = feval(varargin{i}, X);
        X(i) = xi - hi;
        X(j) = xj - hj;
        fm2 = feval(varargin{i}, X);
        X(i) =xi;
        X(j) = xj;
        hessian(i,j) + ((fp1-fm1) - (fp2-fm2))/(4*hi*hj);
      end
    end
  end
  
  % calculate Newton's correction
  Dx = jacobian \ Fx;
  Bx = jacobian \ ((hessian*Dx) .* Dx);
  % Calculate Halley's correction
  Cx = (Dx.*Dx)./(Dx + Bx / 2);
  X = X + Cx;
  if ((norm(X) / realsqrt(N)) <= Xtol) || ((norm(Fx) / realsqrt(N)) <= Ftol)
    bStop = true;
  end
end

The code for test file fx1.m is:

Code:

function y = fx1(x)
y = x(1)*x(2)*x(3)-x(2)^2-15;

The code for test file fx2.m is:

Code:

function y = fx2(x)
y = x(2)^2 - x(1)^3 + x(3)^2 - 17;

The code for test file fx3.m is:

Code:

function y = fx3(x)
y = (x(1)+x(3))/x(2)-2;

Sample run in Matlab is:

[X, ~]=sne_halley(3, [4 2 3],1e-7,1e-8,1000,@fx1,@fx2,@fx3)

You get:

Code:

% 
% X =
% 
%     2.0000    3.0000    4.0000
%

I think the Matlab code can be adapted to work on the HP-39gII, HP Prime, HP-71B with MATH ROM, HP-75C with MATH ROM, and the HP48S/SX/G/G+/Gx/49G/49G+/Gii/50 line of graphing calculators. The code for the nonlinear functions can be consolidated into a single program function that evaluate the different non-linear functions using an IF statement.

Namir