So I have this example Newton's method for root finding with quadratic convergence below. It takes a function f, the derivative of f df, initial guess g, and tolerance tol. It outputs the # iterations it to reach nth root, the root estimate r, and the error approximation err.
function [it, r, err] = QuadraticN(f, df, g , tol)
num_it = 20;
it_max = num_it + 1;
x(1)=(g + num_it)/2; %set x at n = 1
n = 1;
r = 0; % root
it = 0; % iteration #
err =0; % error
% n = 1
r(1) = x(1)-f(x(1))/df(x(1)); % root at n = 1
err(1) = abs(x(1)-g); % error at n = 1
% n > 1
while (min(abs(f(x(n))))) && (abs(x(n)-it_max)>tol) && (it < num_it)
x(n+1)=x(n)-f(x(n))/df(x(n)); % quadratic method
it_max=x(n); % reset counter
r= x(n); % set root to current at n
err = f(x(n)); % set error to current
it = n; % keep track of iterations
n=n+1; % increment to next n
end
end
Given this example, I am trying to implement a second version of the method that uses cubic convergence as opposed to quadratic and this is what I have (this is just the part of the function that would be changed to cubic, it takes the same inputs with the addition of ddf = the second derivative):
r(1) = x(1)-df(x(1))/ddf(x(1)) + (sqrt((df(x(1)))^2 -
2*f(x(1))*ddf(x(1)))) / ddf(x(1));
err(1) = abs(x(1)-g); % error at n = 1
while (min(abs(ddf(x(n))))) && (abs(x(n)-it_max)>tol) && (it < num_it)
cube = ((sqrt((df(x(n)))^2 - 2*f(x(n))*ddf(x(n)))) / ddf(x(n)));
first = x(n) - df(x(n))/ddf(x(n));
first = real(first); % is this right?
x(n+1) = first + (cube);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% FIX THIS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% check when to do plus or minus
if( first > 0)
x(n+1)=x(n)-df(x(n))/ddf(x(n)) + cube;
else
x(n+1)=x(n)-df(x(n))/ddf(x(n)) - cube;
% end
I'm not sure if I am checking when to use +/- cube correctly. I know it will depend on the multiplicity but I'm not sure if using first (Newton's basic method). It is outputting some of the roots correctly, but for instance for the equation -x^3 + 8 with initial guess 2 I am getting -9.85985274E-01. Now I don't know if this is part of the problem above or if converting to real before calculating is. I am also getting negative errors. I am checking in the same way as I did with quadratic but I am not sure if this is right.