%Iterative methods to find solutions to equations:
% 1. Fixed-point iteration
% 2. Newton's method
% 3. Secant method
%
%% FIXED POINT ITERATION: EX 1
clc

fprintf(' We can attempt to solve  x = g(x)  similar to Gauss-Seidel from before. \n')
fprintf('   Step 1: make a guess for x \n')
fprintf('   Step 2: check if  g(x) = x  (probably won''t happen) \n')
fprintf('   Step 3: replace guess by  g(x) \n')
fprintf('   Step 4: repeat \n')
fprintf('      (maybe calculate error = |x-g(x)| to see if you are getting closer \n')
fprintf('       before you repeat too many times...) \n')


%%
fprintf('EX 1: Solve  x = sin(x) + cos(x) \n')

clc

% We will use an INLINE FUNCTION so we don't have to worry about mistyping
% (sin(x) + cos(x)) all of the time.  The format for inline functions is:
%   <name> = inline('<function>', '<variable>', '<variable>', ...)
%  NOTE:  <function> and <variables>, must be in quotes (so that matlab 
%   saves them as text instead of evaluating them to be numbers right away)

g = inline('sin(x) + cos(x)', 'x')

fprintf(' Initial guess: \n')
x = 1

%%
fprintf(' How close are we? \n')
error = abs(x - g(x))

%%
% Evaluate this cell 4 or 5 times so that the error is small.
%
fprintf(' Replace x by g(x): \n')
x = g(x)

error = abs(x - g(x))

%%
fprintf(' Go back to the previous cell and click "Evaluate cell" until the error is small.\n\n')

disp('      x     sin(x)+cos(x)')
disp([x, sin(x)+cos(x)])




%% FIXED POINT ITERATION: EX 2 (while loop)

clc

fprintf('Time to introduce a loop!\n')
fprintf('EX 2: Solve  x = cos(x) + 2*sin(x)\n')

g = inline('cos(x) + 2*sin(x)', 'x')
x = -1

%%
% We will iterate until the error is small.  Just to be safe, we'll also
% count how many times we've iterated and stop if we iterate too much.
%

error = abs(x - g(x))

iteration = 0

% In the command below, (error > .001) will be either TRUE or FALSE
%  also (iteration < 11) is either TRUE or FALSE
% The operation "&&" is called a "logical and": it tells matlab to check
%  that the conditions on the left and on the right are both TRUE

while ((error > .001) && (iteration < 100))
    iteration = iteration + 1
    x = g(x)
    error = abs(x - g(x))
end




%% FIXED POINT ITERATION: EX 3 (importance of initial guess)
%
clc
fprintf('EX 3: Solve sin(x)=1/3\n')

fprintf('  To solve sin(x) = 1/3 using iteration, we need to introduce an x:\n')
fprintf('       x + sin(x) = x + 1/3\n')
fprintf(' x + sin(x) - 1/3 = x\n--------\n')

fprintf('The initial guess has two effects:\n')
fprintf(' 1. How fast we converge to an answer.\n')
fprintf(' 2. What answer we converge to.\n')

fprintf('----\n We will investigate both of these for the equation x = x+sin(x)-1/3\n');

g = inline('x+sin(x)-1/3', 'x')

%%

fprintf('----------\nInitial guess:')
x = -1
iteration=0;
while ((iteration < 1000) && (abs(x-g(x)) > .00001))
    x = g(x); 
    iteration = iteration + 1;
end
disp('Result:')
disp('      x       error   iteration')
disp([x, abs(x-g(x)), iteration])

%%

fprintf('--------\nInitial guess:')
x = 1
iteration=0;
while ((iteration < 1000) && (abs(x-g(x)) > .00001))
    x = g(x); 
    iteration = iteration + 1;
end
disp('Result:')
disp('      x       error   iteration')
disp([x, abs(x-g(x)), iteration])

%%

fprintf('--------\n There''s lots of solutions to sin(x)=1/3: \n')
fprintf('   x = -3.481429563...\n')
fprintf('   x =  0.339836909...\n')
fprintf('   x =  2.801755744...\n')
fprintf('  (see plot) \n')

% To show the three solutions to x=g(x), we will graph y=x and y=g(x)
% together.  We use the matlab command "fplot" which is designed for
% plotting functions.  To use fplot we type:
%  fplot( <function pointer>, [<x limits>, <y limits>] )
%
% Don't worry about <function pointer> now.  Basically, it means either
%  "@<function name>" where <function name> is the name of a matlab function
%                                 (or an m-file with a function you wrote)
% or else 
%  "@(x) <function>"  where <function> is something like x^2 or x*sin(x)
% or else
%  "'<function>'"     where '<function>' is like for inline functions
%
% (You can pretend "pointer" means "its name starts with @ or it is in ''.")
%
% EXAMPLES: 
%   fplot(@cos, [-5,5])              plots y=cos(x) from x=-5 to x=5
%   fplot(@(x) x^2+1, [-2,2,-1,1])   plots y=x^2+1  from x=-2 to x=2
%                                                    and y=-1 to y=1
%   fplot('tan(x)', [-4,5,-3,2])     plots y=tan(x) from x=-4 to x=5
%                                                    and y=-3 to y=2
%
% To plot two things together we use the function pointer:
%    @(x) [<function 1>, <function 2>]
% we could also use 
%    '[<function 1>, <function 2>]'  but inline functions wouldn't work

figure                               % open a figure window
fplot(@(x) [g(x), x], [-10, 10])     % plot x and g(x)
title('y=x and y=x+sin(x)-1/3')      % add a title (PRETTY!)

%%
fprintf('-----\n Graph where initial values between -10 and 10 go!\n')

guess = -10:.01:10;                    % this will be our initial guesses
x = guess;

iteration_count = zeros(1,length(x));  % this will keep track of how quickly
                                       % the error gets small for an initial
                                       % guess

for (iteration = 0:1000)             % iterate 1 time for each arabian night
    x = g(x);
    error = abs(x - g(x));
    iteration_count = iteration_count + (error > .000001); 
end

figure
subplot(3,1,1)       % Divide the plot window into a 3x1 matrix of plots
                     % and make the first plot active
plot(guess,x)   % Plots the vector -10:.01:10 vs. x
title('target of convergence')    % Add a title

subplot(3,1,2)       % Keeps the plot window as a 3x1 matrix of plots
                     % but makes the second plot active
plot(guess,iteration_count)  % Plot the number of iterations
title('number of iterations')     % Add a title

subplot(3,1,3)
plot(guess,[guess', (g(guess))'])
title('solutions')

fprintf(' Note that iteration is only giving us half of the solutions...\n-------\n')




%% FIXED POINT ITERATION: EX 4 (importance of how you write your function)
%
clc
fprintf('EX 4: Solve x^2-3x+1=0\n')

fprintf(' ---> Attempt 1:  use x = (x^2+1)/3  to iterate')
g = inline('(x.^2+1)/3', 'x')

%%
fprintf('--------\nInitial guess:')
x = 1
iteration=0;
while ((iteration < 1000) && (abs(x-g(x)) > .00001))
    x = g(x); 
    iteration = iteration + 1;
end
disp('Result:')
disp('      x       error   iteration')
disp([x, abs(x-g(x)), iteration])

fprintf('That''s one of the roots.. but there should be two...\n')
%%
fprintf('--------\nInitial guess:')
x = -4
iteration=0;
while ((iteration < 1000) && (abs(x-g(x)) > .00001))
    x = g(x); 
    iteration = iteration + 1;
end
disp('Result:')
disp('   x    error   iteration')
disp([x, abs(x-g(x)), iteration])

fprintf('This initial guess led to an explosion...\n')
%%
fprintf('--------\nInitial guess:')
x = 4
iteration=0;
while ((iteration < 1000) && (abs(x-g(x)) > .00001))
    x = g(x); 
    iteration = iteration + 1;
end
disp('Result:')
disp('   x    error   iteration')
disp([x, abs(x-g(x)), iteration])

fprintf('This initial guess led to an explosion...\n')
%%
fprintf('-----\n Graph where initial values between -4 and 4 go!\n')

guess = -4:.001:4;                    % this will be our initial guesses
x = guess;

iteration_count = zeros(1,length(x));  % this will keep track of how quickly
                                       % the error gets small for an initial
                                       % guess

for (iteration = 0:1000)             % iterate 1 time for each arabian night
    x = g(x);
    error = abs(x - g(x));
    iteration_count = iteration_count + (error > .000001); 
end

figure
subplot(3,1,1)       % Divide the plot window into a 3x1 matrix of plots
                     % and make the first plot active
plot(guess,x)   % Plots the vector -10:.01:10 vs. x
title('target of convergence')    % Add a title

subplot(3,1,2)       % Keeps the plot window as a 3x1 matrix of plots
                     % but makes the second plot active
plot(guess,iteration_count)  % Plot the number of iterations
title('number of iterations')     % Add a title

subplot(3,1,3)
plot(guess,[guess', (g(guess))'])
title('solutions')

fprintf(' The other solution is around 2.5.... but no initial guess sees it\n')
fprintf('   instead the other inital guesses explode!\n')

%%


fprintf('------------\n Solve x^2 - 3x + 1 = 0  \n')
fprintf(' ---> Attempt 2:  use x = (3x - 1)/x  to iterate')
g = inline('3 - 1./x', 'x')

%%
fprintf('--------\nInitial guess:')
x = 1
iteration=0;
while ((iteration < 1000) && (abs(x-g(x)) > .00001))
    x = g(x); 
    iteration = iteration + 1;
end
disp('Result:')
disp('      x       error   iteration')
disp([x, abs(x-g(x)), iteration])

fprintf('Now we''re getting the other root from this initial guess...\n')
%%
fprintf('--------\nInitial guess:')
x = -1
iteration=0;
while ((iteration < 1000) && (abs(x-g(x)) > .00001))
    x = g(x); 
    iteration = iteration + 1;
end
disp('Result:')
disp('   x    error   iteration')
disp([x, abs(x-g(x)), iteration])

fprintf('This initial guess also gets the other root...\n')
%%
fprintf('--------\nInitial guess:')
x = -10
iteration=0;
while ((iteration < 1000) && (abs(x-g(x)) > .00001))
    x = g(x); 
    iteration = iteration + 1;
end
disp('Result:')
disp('   x    error   iteration')
disp([x, abs(x-g(x)), iteration])

fprintf('This initial guess still gets the other root...\n')
%%
fprintf('-----\n Graph where initial values between -4 and 4 go!\n')

guess = -4:.001:4;                    % this will be our initial guesses
x = guess;

iteration_count = zeros(1,length(x));  % this will keep track of how quickly
                                       % the error gets small for an initial
                                       % guess

for (iteration = 0:1000)             % iterate 1 time for each arabian night
    x = g(x);
    error = abs(x - g(x));
    iteration_count = iteration_count + (error > .000001); 
end

figure
subplot(2,1,1)       % Divide the plot window into a 2x1 matrix of plots
                     % and make the first plot active
plot(guess,x)   % Plots the vector -10:.01:10 vs. x
title('target of convergence')    % Add a title

subplot(2,1,2)       % Keeps the plot window as a 2x1 matrix of plots
                     % but makes the second plot active
plot(guess,iteration_count)  % Plot the number of iterations
title('number of iterations')     % Add a title

fprintf('  The initial guesses all converge, but none of them sees the earlier root!\n')
disp('-----------------------------------------------------')
disp('PROBLEM: the method of iteration pushes away from points where the derivative')
disp(' of g is bigger than 1.  So if g has a big derivative at the point you''re looking')
disp(' for, then iteration will never find it.')
disp('CONCLUSION:  we need a better method -- maybe our method should involve derivatives')
disp(' since they are what''s screwing us up right now....')







%% NEWTON'S METHOD: Warm-Up 
%
clc
disp('NEWTON''S METHOD:  Warm-Up')
disp('')
disp(' To get ready for Newton''s method, let''s do something a little silly...')
disp(' We will solve  0 = 3x+6  using an iterative method.')
disp('   Attempt      x = x+3x+6   and iterate')

g = inline('x+3*x+6', 'x')
%%
fprintf('--------\nInitial guess:')
x = 1
iteration=0;
while ((iteration < 1000) && (abs(x-g(x)) > .00001))
    x = g(x); 
    iteration = iteration + 1;
end
disp('Result:')
disp('   x    error   iteration')
disp([x, abs(x-g(x)), iteration])

disp(' This is never going to work... the derivative of g is too big.')
%%
disp('------------')
disp('We need to fix the derivative of g....')
disp('g = x + f(x), so we wan''t to make the derivative of f(x) = -1')
disp('     0 = 3x+6   can be changed to')
disp('0/(-3) = (3x+6)/(-3)   (divide by negative of old derivative)')
disp('     0 = -(x+2)')
disp('  Attempt    x = x-(x+2)  and iterate')

g = inline('x-(x+2)', 'x')

%%
fprintf('--------\nInitial guess:')
x = 1
iteration=0;
while ((iteration < 1000) && (abs(x-g(x)) > .00001))
    x = g(x); 
    iteration = iteration + 1;
end
disp('Result:')
disp('   x    error   iteration')
disp([x, abs(x-g(x)), iteration])

disp('------------------------------')
disp('What we just did is newton''s method.  If you aren''t dealing with a linear')
disp('function, then the derivative will change as you move around, so it must be')
disp('recalculated at every step... but otherwise it is the same as this.')

%% NEWTON'S METHOD: EX 1
clc

disp('EX 1:   Solve 0 = x^2-3x+1')

f = inline('x.^2 - 3*x + 1', 'x')    % Our function of interest
df = inline('2*x - 3', 'x')          % The derivative of f

%%
fprintf('--------\nInitial guess:')
x = 1
iteration=0;
while ((iteration < 1000) && (abs(f(x)) > .00001))
    x = x - f(x)./df(x); 
    iteration = iteration + 1;
end
disp('Result:')
disp('       x      error   iteration')
disp([x, abs(x-g(x)), iteration])

%%
fprintf('--------\nInitial guess:')
x = 3
iteration=0;
while ((iteration < 1000) && (abs(f(x)) > .00001))
    x = x - f(x)./df(x); 
    iteration = iteration + 1;
end
disp('Result:')
disp('       x      error   iteration')
disp([x, abs(x-g(x)), iteration])

disp('-----------------------')
disp('NOW WE CAN FIND BOTH OF THE ROOTS!!!!  YAY!')



%% NEWTON'S METHOD : EX 2
%
clc
disp('EX 2:   Solve x^6 + 2x^5 - 3x^4 - 10x^3 - 6x^2 + 8x + 8')

f = inline('x.^6 + 2*x.^5 - 3*x.^4 - 10*x.^3 - 6*x.^2 + 8*x + 8', 'x')
df = inline('6*x.^5 + 10*x.^4 - 12*x.^3 - 30*x.^2 - 12*x + 8', 'x')

%%
%
fprintf('-----\n Graph where initial values between -3 and 3 go!\n')

guess = -3:.01:3;                    % this will be our initial guesses
x = guess;

iteration_count = zeros(1,length(x));  % this will keep track of how quickly
                                       % the error gets small for an initial
                                       % guess

for (iteration = 0:1000)             % iterate 1 time for each arabian night
    x = x-(f(x)./df(x));
    error = abs(f(x));
    iteration_count = iteration_count + (error > .000001); 
end

figure               % Open a new figure window

subplot(2,1,1)       % Divide the plot window into a 2x1 matrix of plots
                     % and make the first plot active                
plot(guess,x)        % Plots x
title('NEWTON''S METHOD ANALYSIS') % Add a title
ylabel('target of convergence')    % Add a y label
xlabel('initial guess')

subplot(2,1,2)       % Keeps the plot window as a 2x1 matrix of plots
                     % but makes the second plot active
plot(guess,iteration_count)  % Plot the number of iterations
ylabel('number of iterations')     % Add a y label
xlabel('initial guess')


disp('--------------------------')
disp('We found all of the solutions... note that it is a bit random which solution')
disp('any particular initial guess will lead to, though... In fact, there is some')
disp('kind of weird fractal stuff going on with which solution you go to.')





%% SECANT METHOD: IDEA
clc

fprintf('The problem with Newton''s method is that it requires you to compute the \n')
fprintf('derivative of your function.... This might be hard.  It also may be \n')
fprintf('impossible -- for example, if you don''t actually have a function, but instead\n')
fprintf('just a bunch of sample points.\n\n')

fprintf('--> To remedy this, we can use the SECANT METHOD.\n')
fprintf(' The idea behind the Secant Method is to just estimate the derivative, instead\n')
fprintf(' of trying to compute it\n\n')

fprintf('--> We can estimate df(x) by computing (f(x+h)-f(x))/h for h small\n')

%%
fprintf('---------\n Check the accuracy of sample derivative:  f(x) = 1+sin(x)cos(x)')

f = inline('.5+sin(x)*cos(x)', 'x')
df = inline('(cos(x))^2 - (sin(x))^2', 'x')
estimated_df = inline('((.5+sin(x+h).*cos(x+h))-(.5+sin(x).*cos(x)))./h', 'x', 'h')
     %% NOTE: We cannot use f(x) in the definition of estimated_df, since
     %% inline functions cannot refer to each other.

fprintf(' Estimate derivative at x=.5:')

actual_derivative = df(.5)

figure
plot(.001:.001:.5, estimated_df(.5,.001:.001:.5))
title('estimated derivative at x=.5 (actual value = .5403)')
ylabel('estimated derivative')
xlabel('value of h')



%% SECANT METHOD: EX 1
clc
fprintf('SECANT METHOD:  EX 1 -- solve  .5+sin(x)*cos(x)=0\n')

% In this example, we are going to define the function f differently...
% instead of defining it as an inline function, we will define it as an
% "anonymous function".  The reason we will do this is that matlab handles
% anonymous functions better than inline ones -- specifically, anonymous
% functions can refer to each-other, but inline functions cannot.  The
% format for anonymous functions is:
% 
% <function name> = @(<variables>) (function)
%
% EX:  f = @(x,y) (x*sin(y+1))
% EX:  g = @(x,y) (x^y - x)
% EX:  diff = @(x,y) (f(x,y) - g(x,y))
% 
% ... some day I'll modify all of the MATH 210 matlab scripts to use anonymous
% functions instead of inline functions...

f = @(x) (.5+sin(x)*cos(x))
estimated_df = @(x) ( (f(x+.01) - f(x))/.01 )

%%
fprintf('initial guess:')
x = 1

%%
fprintf('iterating the secant method...')
iteration = 0;
while (abs(f(x)) > .000001)
    x = x-f(x)/estimated_df(x);
    iteration = iteration + 1;
end
fprintf('done!\n')

disp('    result     error   iterations')
disp([x f(x) iteration])

%%
disp('--> Good things about secant method:  ')
disp(' 1. converges more often than fixed point iteration method')
disp(' 2. don''t need to compute derivatives')
disp('--> Bad things about secand method:  ')
disp(' 1. doesn''t converge as fast as Newton''s method')




%% SECANT METHOD: EX 2
clc
fprintf('Compare the secant method to Newton''s method!')

f = @(x) (x.^6 + 2*x.^5 - 3*x.^4 - 10*x.^3 - 6*x.^2 + 8*x + 8)
estimated_df = @(x) ( (f(x+.01) - f(x))/.01 )

%%


fprintf('-----\n Graph where initial values between -3 and 3 go!\n')

guess = -3:.01:3;                    % this will be our initial guesses
x = guess;

iteration_count = zeros(1,length(x));  % this will keep track of how quickly
                                       % the error gets small for an initial
                                       % guess

for (iteration = 0:1000)             % iterate 1 time for each arabian night
    x = x-(f(x)./estimated_df(x));
    error = abs(f(x));
    iteration_count = iteration_count + (error > .000001); 
end

figure               % Open a new figure window

subplot(2,1,1)       % Divide the plot window into a 2x1 matrix of plots
                     % and make the first plot active                
plot(guess,x)        % Plots x
title('SECANT METHOD ANALYSIS')    % Add a title
ylabel('target of convergence')    % Add a y label
xlabel('initial guess')

subplot(2,1,2)       % Keeps the plot window as a 2x1 matrix of plots
                     % but makes the second plot active
plot(guess,iteration_count)  % Plot the number of iterations
ylabel('number of iterations')     % Add a y label
xlabel('initial guess')

disp(' This plot should look very similar to the one from before....')
disp('When h is small, the secant method should converge to the same root')
disp('as Newton''s method, and it should converge almost as fast.')