%%Newton's Method for Root Finding

% Uses basic calculus to find the root of a function
% Idea: Start with an initial guess x0 and repeatedly update
% x(k+1) = x(k) - f(x(k))/f'(x(k))
% where f'(x(k)) is the first derivative of f evaluated at x(k)

% Why this update?
% Taylor's Theorem: A (smooth) function can be approximated near x
% using the first derivative at x
% f(x(k) + p) approx= f(x(k)) + f'(x(k))p ->
% f(x(k+1)) approx= f(x(k)) + f'(x(k))(x(k+1) - x(k))
% Remember, our goal is to find x(k+1) such that f(x(k+1)) = 0.
% So we can set this equation equal to 0 and solve for x(k+1)!
% 0 = f(x(k)) + f'(x(k))(x(k+1) - x(k)) ->
% x(k+1) = x(k) - f(x(k))/f'(x(k))

% Aside: Similar idea that can be applied to many areas. One example is
% for minimizing a function f (i.e. find x such that f(x) is minimized).
% In this case, we'd set x(k+1) = x(k) - c f'(x(k))
% f(x(k+1)) approx= f(x(k)) - c (f'(x(k))^2
% Known as gradient descent, basic algorithm in optimization.
% Applies to problems in machine learning/statistics, economics, operations
% research, engineering, etc. Check out CS/ECE/ISye 524 for more info.
% In MATLAB, if you want to minimize a function, call fminunc.

% Can use calculus to show that the error decreases at a rate that is
% quadratic in the previous error, aka a quadratic rate of convergence.
% Note: this is only guaranteed to work for smooth, continuously
% differentiable functions.

%Example: f(x) = 6 - 2x^2
%True roots are at +-sqrt(3)
f = @(x)6 - 2*x.^2;

%First step: plot to find a good starting point
xx = -3 : 0.01 : 3;
yy = f(xx);
ymin = -7;
ymax = 7;
xmin = -3;
xmax = 3;
plot(xx,yy,'b',[0,0],[ymin,ymax],'k',[xmin,xmax],[0,0],'k')
%Adds a grid to help see exactly where the line is
grid on
axis([xmin xmax ymin ymax])

%Need the first derivative to call Newton's method
fprime = @(x)-4*x;
xstart = 1;
epsilon = 1e-12;
x = newton(f,fprime,xstart,epsilon);

precision = 20;

% display successive approximations from Newton's method
for k = 1 : length(x)
    disp(['x(', num2str(k), ') = ', num2str(x(k),  precision)])
end

x_root = x(length(x));
disp(['root = ', num2str(x_root, precision)])
perc_error=abs(sqrt(3)-x_root)/sqrt(3)*100;
disp(['% error = ' num2str(perc_error, precision) ' %'])

