Post Reply 
Solving a system on nonlinear equations using Halley's method
04-05-2014, 04:44 AM (This post was last modified: 11-09-2015 09:23 AM by Namir.)
Post: #1
Solving a system on nonlinear equations using Halley's method
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
Find all posts by this user
Quote this message in a reply
04-05-2014, 07:46 AM
Post: #2
RE: Solving a system on nonlinear equations using Halley's method
(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
Find all posts by this user
Quote this message in a reply
06-06-2014, 01:48 AM
Post: #3
RE: Solving a system on nonlinear equations using Halley's method
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.)
Find all posts by this user
Quote this message in a reply
06-06-2014, 03:02 AM
Post: #4
RE: Solving a system on nonlinear equations using Halley's method
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
Find all posts by this user
Quote this message in a reply
06-10-2014, 02:43 AM
Post: #5
RE: Solving a system on nonlinear equations using Halley's method
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.
Find all posts by this user
Quote this message in a reply
06-12-2014, 01:39 PM (This post was last modified: 06-12-2014 01:46 PM by Eddie W. Shore.)
Post: #6
RE: Solving a system on nonlinear equations using Halley's method
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.
Visit this user's website Find all posts by this user
Quote this message in a reply
Post Reply 




User(s) browsing this thread: 1 Guest(s)