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 Namir Senior Member Posts: 702 Joined: Dec 2013
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
04-05-2014, 07:46 AM
Post: #2
 Thomas Klemm Senior Member Posts: 1,447 Joined: Dec 2013
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
06-06-2014, 01:48 AM
Post: #3
 ttw Member Posts: 198 Joined: Jun 2014
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.)
06-06-2014, 03:02 AM
Post: #4 Namir Senior Member Posts: 702 Joined: Dec 2013
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
06-10-2014, 02:43 AM
Post: #5
 ttw Member Posts: 198 Joined: Jun 2014
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.
06-12-2014, 01:39 PM (This post was last modified: 06-12-2014 01:46 PM by Eddie W. Shore.)
Post: #6 Eddie W. Shore Senior Member Posts: 1,112 Joined: Dec 2013
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.
 « Next Oldest | Next Newest »

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