0001 function [ x,k ] = Halley( f,initial,tol,max )
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020 if tol<=0;
0021 error('tolerance must be positive');
0022 elseif max<=0;
0023 error('maximum number of iterations must be positive');
0024 elseif isa(f,'sym');
0025 df1 = diff(f);
0026 df2 = diff(df1);
0027 f = matlabFunction(f);
0028 df1 = matlabFunction(df1);
0029 df2 = matlabFunction(df2);
0030 elseif isa(f,'double');
0031 [n,m]=size(f);
0032 if n~=1 && m~=1;
0033 error(['input must be a symbolic function or vector holding ',...
0034 'polynomial coefficients']);
0035 else
0036 v=poly2sym(f);
0037
0038 df1 = diff(v);
0039 df2 = diff(df1);
0040 f=matlabFunction(v);
0041
0042 df1 = matlabFunction(df1);
0043 df2 = matlabFunction(df2);
0044 end
0045 else
0046 error(['input must be a symbolic function or vector holding ',...
0047 'polynomial coefficients']);
0048 end
0049
0050
0051 x=zeros(max,1);
0052 x(1)=initial;
0053 fx=feval(f,x(1));
0054 df1x=feval(df1,x(1));
0055 df2x=feval(df2,x(1));
0056 k=2;
0057 l=0;
0058 d0=0;
0059
0060 while abs(fx)>tol && k<=max;
0061
0062 denom=2*(df1x)^2-fx*df2x;
0063 if abs(denom)<1e-12;
0064 disp('denominator vanishes');
0065 x=x(1:k); return;
0066 end
0067
0068
0069 x(k)=x(k-1)-(2*fx*df1x)/(denom);
0070
0071 test=find(x(1:k-1)==x(k));
0072 if isempty(test)==0;
0073
0074
0075 disp('algorithm fell into infinite loop');
0076 x=x(1:k); return;
0077 end
0078
0079 d1=abs(x(k)-x(k-1));
0080 if d1>d0;
0081 l=l+1;
0082 if l>20;
0083 disp('distance between succesive iterates is increasing');
0084 x=x(1:k); return;
0085 end
0086 else
0087 l=0;
0088 end
0089 d0=d1;
0090
0091 fx=feval(f,x(k));
0092 df1x=feval(df1,x(k));
0093 df2x=feval(df2,x(k));
0094
0095 k=k+1;
0096 end
0097
0098 if k==max+1;
0099 disp('maximum number of iterations reached');
0100 end
0101 x=x(k-1);
0102
0103 end