%% Lab 1
% Table of Contents
% Summary of Lab 1
% Linearization
% Newton's Metod
%     Uppgift 1 - Find roots using Newton's method
% fzero()
%     Uppgift 2 - find roots using fzero()
% Extrema (Stationary Points/Maxima and Minima)
%     Uppgift 3 - find stationary points using fminbnd()
% Numerical Integration
%     Uppgift 4 - approximate integral using Riemann sums
%     Illustration of different Riemann Sums
%     Convergence of Riemann sums
%     Error Bounds
%     Uppgift 5 - Undamped Pendulum (solve parameterized function)
% Numerical Derivation
%     Uppgift 6
%
%% Summary of Lab 1
% Roots using Newton's method and fzero().
% Extrema using fminbnd().
% Integrals using Riemann sums and integral(). We consider the rate of convergence O\left(h^n \right) and the error bounds of each method.
% Derivatives using difference quotients. We consider the convergence rates depending on step-size h.
%
%% Linearization
% We wish to linearize a differentiable function around a point a. That is, we take the tangent at a. More generally, the first two terms of the Taylor series.
% L\left(x\right)=f\left(a\right)+f^{\prime } \left(a\right)\left(x-a\right)
% With this, we get a line that approximates f\left(x\right) near a.

close all
f  = @(x) 0.5*(x-2).^2-2*cos(2*x)-1.5;
Df = @(x) x-2+4*sin(2*x);
L  = @(x,a) f(a)+Df(a)*(x-a); % Tangent line of curve at a.
a=4;

x = linspace(-3,7);
plot(x, f(x), 'b', x, L(x, a), 'r')
ylim([-5 10]); grid on;
title(['Linearization of f(x) around a=' num2str(a)])

%
%% Newton's Metod
% Newton's method follows directly from the linearization of a function, were solve for L\left(x\right)=0.
% It is an interative algorithm were we initially solve for a = x0, an estimate, then use the result as the new point to linearize about.
% L\left(x\right)=f\left(x_k \right)+f^{\prime } \left(x_k \right)\left(x-x_k \right)=0
% \Rightarrow x_{k+1} =x_k -\frac{f\left(x_k \right)}{f^{\prime } \left(x_k \right)}
% We stop when we have converged (to within a predefined tolerance).
% {|x}_{k+1} -x_k |\le \textrm{tol}
% There are some obvious limitations: the initial estimate must be close to the actual root or the algorithm may diverge, stationary points always result in divergence.
% We implement this in a function (at the bottom of this file) as follows:
% function x = Newton(f, fprim, x0, tol, kmax)
% % Newton's method
% % f: f(x)
% % f: f'(x)
% % x0: initial value
% % tol: convergence tolerance (break condition)
% % kmax: maximum number of iterations
%     xk = x0;
%     for k = 1:kmax
%         x = xk-f(xk)/fprim(xk);
%         if abs(x-xk) <= tol
%             break
%         else
%             xk = x;
%         end
%     end
% end
%% Uppgift 1 - Find roots using Newton's method
% Find the roots of x^3 -\cos \left(4x\right) using newton's method.

close all
f  = @(x) x.^3 - cos(4*x);      % Set up the functions
Df = @(x) 3*x.^2 + 4*sin(4*x);

x = linspace(-2, 2, 100);       % Plot f(x) and estimate roots.
plot(x, f(x), 'b')
grid on, hold on

kmax = 10;      % Set termination parameters for Newton's method. 
tol = 0.5e-8;

x0_vec = [-1.2, -0.5, 0.5];     % Vector containing estimated roots. 
roots = [zeros(size(x0_vec))];  % Empty vector for results.

for k = 1:length(x0_vec)
    roots(k) = Newton(f, Df, x0_vec(k), tol, kmax); % Newton() is described above, and implemented at the bottom of the file.
end
plot(roots, zeros(size(roots)), 'ko')   % Plot roots.

% Just for fun, here are the iterations visualized (use the slider to adjust the inital estimate x0):

x0=-1.8;
for k = 1:kmax
    plot(x0, 0, 'r*')
    plot([x0 x0], [0 f(x0)], 'r--')
    d=-f(x0)/Df(x0);
    xold = x0;
    x0=x0+d;
    
    plot([xold x0], [f(xold) 0], 'r-')
    
    if abs(d)<tol, break, end
end
title('Newton''s Method (with visualization)')

%
%% fzero()
% fzero(fun, x0) is an inbuilt function that is similar to Newton's method. x0 can be a scalar, or a vector containing two points where the sign of f(x) changes.

close all
f = @(x)0.5*(x-2).^2-2*cos(2*x)-1.5;

x = linspace(-5,5,100);
plot(x, f(x))
grid on

root = fzero(f, [-1 0]); % Use interval where endpoints change sign (ie. the interval passes zero).
hold on
plot(root, 0, 'ro')
title('fzero demo')

%
%% Uppgift 2 - find roots using fzero()
% Find where polar function r\left(\phi \right), 0\le \phi \le 2\piintersects unit circle (ie. r\left(\phi \right)=1)

close all

rho = @(t) (2+sin(3*t))./(sqrt(1+exp(cos(t))));
f   = @(t) rho(t) - 1;  % Modify function to give r(phi) = 1 intersection.

phi = linspace(0, 2*pi, 1000);
polarplot(phi, rho(phi));
pax = gca;                      % Get Current Axis. This lets us modify graph properties.
pax.ThetaAxisUnits = 'radians'; % Radians are way nicer.
hold on

x0    = deg2rad([0 69 100 200 225 309]);    % Visual inspection gives the following approximate roots.
roots = zeros(size(x0));                    % Allocate vector for solved roots.

for k = 1:length(roots)
    roots(k) = fzero(f, x0(k)); % Find the roots using fzero.
end

polarplot(roots, ones(size(roots)), 'ro')
polarplot(phi, ones(size(phi)), 'r--')

%
%% Extrema (Stationary Points/Maxima and Minima)
% Previously, we linearized a function in order to describe the slope. To discriminate between minima and maxima, we need a higher order approximation in order to describe the curvature. This information is, again, contained in the Taylor series approximation.
% f\left(x\right)\approx f\left(a\right)+f^{\prime } \left(a\right)\left(x-a\right)+\frac{{f^{\prime } }^{\prime } \left(a\right)}{2}\left(x-a\right)²
% where \textrm{sign}\left\lbrack {f^{\prime } }^{\prime } \left(a\right)\right\rbrack is the relevant discriminant.
% {f^{\prime } }^{\prime } \left(a\right)>0: minimum
% {f^{\prime } }^{\prime } \left(a\right)<0: maximum
%% Uppgift 3 - find stationary points using fminbnd()
% Now we'll disregard all that and use the builtin function fminbnd(fun, xa, xb), where xa and xb are interval bounds.

close all 

f = @(x) (1 + x.^2 - 1.5.*x.^3 + 0.5.*x.^4)./(1 + x.^4);

x = linspace(-4,4);
plot(x, f(x)); 
grid on; 
hold on;

[x, fval] = fminbnd(f, -0.1, 0.1); % Find minimum.
plot(x, fval, 'r*')

[x, fval] = fminbnd(f, 1, 3);
plot(x, fval, 'r*')

h = @(x) -f(x);     % We invert function f in order to use fminbnd to find maxima instead.
                    % We can not modify a function handle directly, we must create a new function handle. 
[x, fval] = fminbnd(h,-1.2,-0.8); % Return coordinate of extrema.
plot(x, -fval, 'm*') % We found local minima of INVERTED function, thus flip it again (-fval).

[x, fval] = fminbnd(h, 0.3, 0.4);
plot(x, -fval, 'm*')

%
%% Numerical Integration
% Here we examine how integrals can be approximated numerically using Riemann sums.
% All the following methods divide the problem-domain into n fixed-size intervals \Delta x=\frac{b-a}{n}. The basic idea is:
% \int_a^b f\left(x\right)\textrm{dx}=\sum_{i=1}^n \int_{x_{i-1} }^{x_i } f\left(x\right)\textrm{dx}
% We are now interested in how best to approximate the integral for each segment. The basic methods are derived geometrically:
% Left rectangle L\left(x\right)=\Sigma_i f\left(x_{i-1} \right)\Delta x
% Right rectangle R\left(x\right)=\Sigma_i f\left(x_i \right)\Delta x
% Midpoint M\left(x\right)=\Sigma_i f\left(\frac{x_{i-1} +x_i }{2}\right)\Delta x
% Trapezeum T\left(x\right)=\Sigma_i \left(\frac{f\left(x_{i-1} \right)+f\left(x_i \right)}{2}\right)\Delta x
% We implement these as functions directly:

L = @(f,x,dx) dx.*sum(f(x(1:end-1)));
R = @(f,x,dx) sum(f(x(2:end)).*dx);
M = @(f,x,dx) sum(f((x(1:end-1) + x(2:end))./2).*dx);
T = @(f,x,dx) sum(((f(x(1:end-1)) + f(x(2:end)))./2).*dx);

%% Uppgift 4 - approximate integral using Riemann sums
% Approximate \int_0^1 \textrm{xsin}\left(x\right)\textrm{dx} using the above Riemann sums.

f = @(x) x.*sin(x);

a=0; 
b=1; 
n=10;
dx = (b-a)/n;
x = linspace(a, b, n+1);

format long
["L", L(f,x,dx); "R",R(f,x,dx); "M",M(f,x,dx); "T",T(f,x,dx); "Control", integral(f,a,b) ]

%% Illustration of different Riemann Sums

a=0; 
b=10; 
n=10;

x = linspace(a,b,n+1);
xfine = linspace(a,b,1000);

subplot(4,2,1:2)
plot(xfine, f(xfine), 'r-'), hold on, grid on
for i=1:n
    plot(x(i), f(x(i)), 'g*')
    plot([x(i) x(i+1)], [f(x(i)) f(x(i))], 'g-')
    plot([x(i+1) x(i+1)], [f(x(i)) f(x(i+1))], 'g-')
end
title('Vänster rektangelregel');

subplot(4,2,3:4)
plot(xfine, f(xfine), 'r-'), hold on, grid on
for i=1:n
    plot(x(i+1), f(x(i+1)), 'g*')
    plot([x(i) x(i+1)], [f(x(i+1)) f(x(i+1))], 'g-')
    plot([x(i) x(i)], [f(x(i)) f(x(i+1))], 'g-')
end
title('Höger rektangelregel');

subplot(4,2,5:6)
hold on
for i=1:n
    plot((x(i)+x(i+1))/2, f((x(i)+x(i+1))/2), 'g*')
    plot([x(i) x(i)], [f(x(i)), f((x(i)+x(i+1))/2)], 'g-')
    plot([x(i) x(i+1)], [f((x(i)+x(i+1))/2),  f((x(i)+x(i+1))/2)], 'g-')
    plot([x(i+1) x(i+1)], [f((x(i)+x(i+1))/2), f(x(i+1))], 'g-')
end
plot(xfine, f(xfine), 'r-'), grid on
title('Mittpunktsregeln rektangelregel');

subplot(4,2,7:8)
plot(xfine, f(xfine), 'r-'), hold on, grid on
for i=1:n
    plot([x(i) x(i+1)], [f(x(i)) f(x(i+1))], 'g-')
end
title('Trapetzregel');

%
%% Convergence of Riemann sums

close all
f = @(x) x.*sin(x);

n = 10:10:100;
qv = zeros(size(n));
qh = zeros(size(n));
qm = zeros(size(n));
qt = zeros(size(n));

for k = 1:length(n)
    x = linspace(a,b,n(k)+1);
    dx = (b-a)/n(k);

    qv(k) = L(f,x,dx);
    qh(k) = R(f,x,dx);
    qm(k) = M(f,x,dx);
    qt(k) = T(f,x,dx);
end

kontroll = integral(f,a,b);
plot(n, qh-kontroll, n,qv-kontroll, n,qm-kontroll, n,qt-kontroll)
legend('L(n)', 'R(n)', 'M(n)', 'T(n)')
title('Error vs step-size')
xlabel('n'); ylabel('Err(n)')
grid on

% Conclusion: H and L converge at same rate, though H overestimates and L underestimates.
% Likewise, M and T converge with n^2 , though M overestimates and T underestimates (see Stewart pg. 509, the rectangle formed by the Midpoint rule is equal in area to a trapezoid that is more accurate than that of the Trapezoidal rule, ironically. Twice as accurate.)
%% Error Bounds
% L and H use rectangles to approximate f, thus the error will be proportionate to the slope of the function, f^{\prime }.
% M and T can both be interpreted as trapezoidal, and thus compensate for slope f^{\prime }, but not for curvature {f^{\prime } }^{\prime }. The error is thus proportionate to {f^{\prime } }^{\prime }.
% Maximum error bounds can be given in terms of maximum slope and maximum curvature.
%
%% Uppgift 5 - Undamped Pendulum (solve parameterized function)
% Undamped pendulum. Calculate the period time resulting from different initial angles T\left(\phi_0 \right).
%
% Here, we are asked to calclate a function of the form:
% f\left(x\right)=\int_c^d h\left(t,x\right)\textrm{dt}
% h(t, x) is a parameterized function ie. it has a controllable parameter x that is different from the variable of integration.
% In Matlab,  integral() expects a function with only one parameter: the variable of integration. Thus, we must create a new (anonymous) function that restricts it to only a single parameter (that is, we assign fixed values to the other parameters).
% h = @(t,x) ...
% integral(@(t) h(t,0), c, d)
% The above will now integrate over h(t), with x=0.
% So, in our example:
% T\left(\phi_0 \right)=\int_c^d h\left(\phi ,\phi_0 \right)\;\;d\phi

close all 

h = @(phi, phi0) 1./sqrt(1 - sin(phi0/2).^2.*sin(phi).^2); % The integrand with two parameters, the former of which is the variable of integration. 

L = 0.1;
g = 9.80665; 
c = 0; 
d = pi/2;

phi0 = 10:1:170;
phi0 = deg2rad(phi0);

T = zeros(size(phi0)); 
for i = 1:length(phi0) 
    T(i) = 4*sqrt(L/g).*integral(@(phi) h(phi, phi0(i)), c, d); % phi is symbolic, whereas phi0 is a fixed constant that we have defined. We must tell integral() which variable to integrate over. In other words, we create a new anonymous function that restricts it to only one parameter, and give the other parameters fixed values.
end

plot(rad2deg(phi0), T)
xlabel('\phi_0');ylabel('T(\phi_0)');
grid on

%
%% Numerical Derivation
% We compute the derivatives numerically using difference quotients. Remember, the definition of the derivative is usually:
% f^{\prime } \left(x\right)=\lim_{h\to 0} \frac{u\left(x+h\right)-u\left(x\right)}{h}
% We are forced to choose a finite value for h. Many different difference quotients with differing convergence rates (as h becomes smaller) can be directly derived from Taylor Series.
% Recap: the Taylor Series is the expansion of a function f\left(x\right) into an infinite polynomial, around a specific point a. It is a weighted sum of the sum of a function and its derivatives at a single point that gives the function itself around that point. The Taylor Series of a function around point a is as follows:
% f\left(x\right)=\sum_{n=0}^{\infty } \frac{f^{\left(n\right)} \left(a\right)}{n!}{\left(x-a\right)}^n
% The most trivial are the backward- and forwards difference quotients:
% D_+ \left(x\right)=\frac{u\left(x+h\right)-u\left(x\right)}{h} (1)
% D_- \left(x\right)=\frac{u\left(x\right)-u\left(x-h\right)}{h} (2)
% The Taylor Series of u\left(x+h\right)at point x is
% u\left(x+h\right)=u\left(x\right)+u^{\prime } \left(x\right)\left(x+h-x\right)+\frac{{u^{\prime } }^{\prime } \left(x\right)}{2}{\left(x+h-x\right)}^2 +\ldotp \ldotp \ldotp
% We substitute this result into (1) to obtain:
% D_+ \left(x\right)=u^{\prime } \left(x\right)+{u^{\prime } }^{\prime } \left(x\right)\frac{h}{2}+\ldotp \ldotp \ldotp =u^{\prime } \left(x\right)+O\left(h\right)
% The same result is obtained for D_- \left(x\right).
%
% Here we approximate the derivative \frac{d}{\textrm{dx}}e^x at x=0\ldotp 4, for different step-sizes h.
% In the final plot we see the effect that rounding error has on small h.

close all

u  = @(x) exp(x);
Du = @(x) exp(x);
DDu = @(x) exp(x);

h = logspace(-15,-1,400);
x=0.4;
subplot(1,4,1:2)
loglog(h, abs((u(x+h)-u(x))./h - Du(x)), 'r'), hold on
loglog(h,abs((u(x+h)-u(x-h))./(2*h) - Du(x)),'g')
loglog(h,abs((-u(x+2*h)+8*u(x+h)-8*u(x-h)+u(x-2*h))./(12*h) - Du(x)),'b')

DD = @(x, h) (u(x+h)-2.*u(x)+u(x-h))./(h.^2);
loglog(h, abs(DD(x,h)-DDu(x)))

legend('Du', 'D_0', 'D_*', 'DD')
xlabel('h'); ylabel('Err(h)')
title('Error vs h')
ylim([10^-20, 10^-1])

digits(32)
h = vpa(logspace(-15,-1,400)); % vpa() is Variable-Precision Arithmetics, which allows us to specify how many digits of precision we want (standard is 16).
subplot(1,4,3:4)
loglog(h, abs((u(x+h)-u(x))./h - Du(x)), 'r'), hold on
loglog(h,abs((u(x+h)-u(x-h))./(2*h) - Du(x)),'g')
loglog(h,abs((-u(x+2*h)+8*u(x+h)-8*u(x-h)+u(x-2*h))./(12*h) - Du(x)),'b')
loglog(h, abs(DD(x,h)-DDu(x)))
legend('Du', 'D_0', 'D_*', 'DD')
xlabel('h'); ylabel('Err(h)')
title('Error vs h (32 bit precision)')
ylim([10^-20, 10^-1])

%
%% Uppgift 6
% Load the table of (x,y) values, find the first and second derivatives, and plot the curvature:
%
% Here we approximate the second derivative with
% D_+ D_- u\left(x\right)=\frac{u\left(x+h\right)-2u\left(x\right)+u\left(x-h\right)}{h^2 }
% This works since, if we substitute the Taylor expansions of u\left(x\pm h\right), we are left with {u^{\prime } }^{\prime } \left(x\right)+O\left(h^2 \right).

close all

tab = load('labtabell.mat');
x = tab.x;
u = tab.y;
h = x(2) - x(1);

% Forwards difference
Dforward = (u(2:end) - u(1:end-1))./h;
% Backward difference quotient is no different in this situation.

% Central difference quotient
D0 = (u(3:end) - u(1:end-2))./(2*h);


% The derivative is essentially evaluated at the midpoint of each pair of
% points. Thus we create a new x-axis with only the midpoints (and one less
% point overall).
xm = (x(1:end-1) + x(2:end))./2;
plot(x, u), hold on, grid on
%plot(xm, Dforward, 'b--')
plot(x(2:end-1), D0, '--')

DD = (u(3:end) - 2.*u(2:end-1) + u(1:end-2))./h.^2;
plot(x(2:end-1), DD, '-.')

k = DD./( (1+D0.^2).^(3/2));
plot(x(2:end-1), k)

%

function x = Newton(f, fprim, x0, tol, kmax)
% Newton's method
% f: f(x)
% f: f'(x)
% x0: initial value
% tol: convergence tolerance (break condition)
% kmax: maximum number of iterations
    xk = x0;
    for k = 1:kmax
        x = xk-f(xk)/fprim(xk);
        if abs(x-xk) <= tol
            break
        else
            xk = x;
        end
    end
end