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

[Thomas Klemm]

(04-05-2014 04:44 AM)Namir Wrote:  I think the Matlab code can be adapted to work on the HP Prime.

This post probably fits better the forum HP Prime.

Kind regards
Thomas

[ttw]

Halley's method does converge faster (as do several others), but the domain of convergence is usually much smaller that that of Newton's method (which can be rather small itself.)

[Namir]

Let me see if I am reading you correctly. Using Newton's method for nonlinear equation with multiple solution may not be (easily) able to converge to the solutions you are looking for? If this is what you mean, than I am in full agreement with you. HP Prime has problems converging to certain solutions of nonlinear equations. I recently wrote a similar Newton's method for 3 eqns on the HP-85 and I came across the same problem. This leads me to think if Newton's (and Halley's) methods are highly affected by the curvature of the nonlinear functions.

Yes? No?

Namir

[ttw]

Both Newton's and Halley's methods (and Chebychev's and all the other higher order methods) strongly depend on the shape of the curves. The topic is usually covered in numerical analysis classes.

[Eddie W. Shore]

When I saw Halley, I immediately thought of Halley's Comet. It turns out that it is the same person who Halley's Comet is named after: Edmund Halley.