0001 function [ x,k ] = inv_quad_it( f,x0,x1,x2,tol,max )
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013 if tol<=0;
0014 error('tolerance must be positive');
0015 elseif max<=0;
0016 error('maximum number of iterations must be positive');
0017 elseif isa(f,'function_handle');
0018
0019 elseif isa(f,'double');
0020 [n,m]=size(f);
0021 if n~=1 && m~=1;
0022 error(['input must be a function handle or vector holding ',...
0023 'polynomial coefficients']);
0024 else
0025 v=poly2sym(f);
0026
0027 f=matlabFunction(v);
0028
0029 end
0030 else
0031 error(['f must be a function handle or vector holding ',...
0032 'polynomial coefficients']);
0033 end
0034
0035
0036 x=zeros(max,1);
0037 x(1)=x0;x(2)=x1;x(3)=x2;
0038 f0=feval(f,x(1));
0039 f1=feval(f,x(2));
0040 f2=feval(f,x(3));
0041 k=4;
0042 l=0;
0043 d0=abs(x2-x1);
0044
0045
0046 while abs(f2)>tol && k<=max;
0047
0048 if abs(f0-f1)<1e-12;
0049 disp('denominator vanishes');
0050 x=x(1:k); return;
0051 elseif abs(f1-f2)<1e-12;
0052 disp('denominator vanishes');
0053 x=x(1:k); return;
0054 elseif abs(f0-f2)<1e-12;
0055 disp('denominator vanishes');
0056 x=x(1:k); return;
0057 end
0058
0059
0060 x(k)=x0*f1*f2/((f0-f1)*(f0-f2))+f0*x1*f2/((f1-f0)*(f1-f2))+...
0061 f0*f1*x2/((f2-f0)*(f2-f1));
0062
0063 d1=abs(x(k)-x(k-1));
0064 if d1>d0;
0065 l=l+1;
0066 if l>max/10;
0067
0068
0069 disp('distance between succesive iterates is increasing');
0070 x=x(1:k); return;
0071 end
0072 else
0073 l=0;
0074 end
0075 d0=d1;
0076
0077 x0 = x1;
0078 x1 = x2;
0079 x2 = x(k);
0080 f0=f1;
0081 f1=f2;
0082 f2=feval(f,x(k));
0083
0084 k=k+1;
0085 end
0086
0087 if k==max+1;
0088 disp('maximum number of iterations reached');
0089 end
0090 x=x(k-1);
0091
0092 end