 %stuff
% close all
hold off
clear all
%% 2 Nollstï¿½llen
f=@(x)0.5*(x-2).^2-2*cos(2*x)-1.5;
x=linspace(-3,7);
plot(x, f(x))
grid on

%% 2.1 linjärisering
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);
a=2.5;
x=linspace(-3,7);
plot(x, f(x), 'b', x, L(x, a), 'r')
ylim([-5 10])
grid on

%% 2,2 Newtons Metod

%% Uppgift 1
%Vill hitta 3 lï¿½sningar
%Approx. ett : 0,5
f=@(x) x.^3 - cos(4*x);
Df=@(x) 3*x.^2 + 4*sin(4*x);
L=@(x,a) f(a)+Df(a)*(x-a);
x = linspace(-2, 2);
plot(x, f(x), 'b')
hold on

kmax = 10; tol = 0.5e-8;

x0 = -0.8;
for k = 1:kmax
    d=-f(x0)/Df(x0);
    x0=x0+d;
    disp([x0 d])
    if abs(d)<tol, break, end
end
plot(x0, 0, 'o')

x0 = -0.3;
for k = 1:kmax
    d=-f(x0)/Df(x0);
    x0=x0+d;
    disp([x0 d])
    if abs(d)<tol, break, end
end
plot(x0, 0, 'o')

x0 = 0.5;
for k = 1:kmax
    d=-f(x0)/Df(x0);
    x0=x0+d;
    disp([x0 d])
    if abs(d)<tol, break, end
end
plot(x0, 0, 'o')

grid on
legend('f(x)')

%% 2.4 fzero
f=@(x)0.5*(x-2).^2-2*cos(2*x)-1.5;
x0 = 4;
x=fzero(f, x0)

%% Uppgift 2 - in three different ways!
% Polar function
rho=@(t) (2+sin(3*t))./(sqrt(1+exp(cos(t))));
rhominusone=@(t) (2+sin(3*t))./(sqrt(1+exp(cos(t))))-1; % decrease radius by 1 to bring unit circle to 'zero crossing'.
theta = 0:0.01:2*pi; % Oberoende variabeln theta, i rad.
theta_deg = theta.*180/pi; %omvandla till grader

% Zero crossings (ie. unit circle intercepts) in rad. Based on
% approximations in degrees!
x1=fzero(rhominusone, deg2rad(60))
x2=fzero(rhominusone, deg2rad(100))
x3=fzero(rhominusone, deg2rad(190))
x4=fzero(rhominusone, deg2rad(220))
x5=fzero(rhominusone, deg2rad(300))
x6=fzero(rhominusone, deg2rad(350))
disp(wrapToPi(x5))

subplot(2,2,1) % r vs. theta (linera plot)
plot(theta, rho(theta));
set(gca,'XTick',0:pi/6:2*pi)
set(gca,'XTickLabel',{'0','\pi/6','\pi/3','\pi/2','2\pi/3','5\pi/6','\pi','7\pi/6','4\pi/3','3\pi/2','5\pi/3','11\pi/6','2\pi'})
xlabel('\Theta');ylabel('r');
grid on;hold on

plot(x1, rho(x1), 'ro')
plot(x2, rho(x2), 'ro')
plot(x3, rho(x3), 'ro')
plot(x4, rho(x4), 'ro')
plot(x5, rho(x5), 'ro')
plot(x6, rho(x6), 'ro')
hold off

subplot(2,2,2) %Polar plot 
polarplot(theta, rho(theta));
pax = gca;
pax.ThetaAxisUnits = 'radians';
hold on
polarplot(theta, sin(theta).^2+cos(theta).^2, '--'); %enhetscirkeln
polarplot(x1, rho(x1), 'ro')
polarplot(x2, rho(x2), 'ro')
polarplot(x3, rho(x3), 'ro')
polarplot(x4, rho(x4), 'ro')
polarplot(x5, rho(x5), 'ro')
polarplot(x6, rho(x6), 'ro')
hold off

subplot(2,2,3) % Rec. plot
x=(rho(theta)).*cos(theta);
y=(rho(theta)).*sin(theta);
plot(x, y)
grid on; hold on; xlabel('rcos(\Theta)');ylabel('rsin(\Theta)');
plot(cos(theta),sin(theta), 'r--')
axis equal
r=@(x)rho(x);
plot(r(x1)*cos(x1),r(x1)*sin(x1),'ro')
plot(r(x2)*cos(x2),r(x2)*sin(x2),'ro')
plot(r(x3)*cos(x3),r(x3)*sin(x3),'ro')
plot(r(x4)*cos(x4),r(x4)*sin(x4),'ro')
plot(r(x5)*cos(x5),r(x5)*sin(x5),'ro')
plot(r(x6)*cos(x6),r(x6)*sin(x6),'ro')
hold off

%% 3 extrempunkter
f=@(x)(x+1.5*x.^2+0.1*x.^3).*exp(-0.7*x.^2);
x=linspace(-4,4);
plot(x,f(x)); grid on

x = fminbnd(f, -0.5, 0)
h=@(x) -f(x) % negativa f(x) tillï¿½ter oss att hitta maxpkter med fminbnd().
x=fminbnd(h,-1.6,-1.2)
x=fminbnd(h,1,1.5)

%% Uppgift 3
f=@(x) (1+x.^2-1.5*x.^3+0.5*x.^4)./(1+x.^4);
Df=@(x) (2*x-4.5*x.^2-2*x.^3-2*x.^5+1.5*x.^6)./((1+x.^4).^2);
DDf=@(x) -(3*x.^9 - 6*x.^8 - 10*x.^6 - 36*x.^5 + 24*x.^4 + 6*x.^2 + 9*x - 2)./(1 + x.^4).^3;
x=linspace(-20,20, 1000);
plot(x, f(x), x, Df(x), x, DDf(x))
xticks(-20:1:20);
grid on; hold on; legend('f(x)', 'fï¿½(x)', 'fï¿½ï¿½(x)')
% Hitta extrempunkter (dvs nollstï¿½llen i f') med fminbnd().
h=@(x)-f(x);
x1=fminbnd(h, -1.3, -1) % Glob. max.
x2=fminbnd(f, -.1, .1)  % lok. min.
x3=fminbnd(h, 0, 0.7) % lok. max
x4=fminbnd(f, 1.5, 2.5) % glob. min
% Berï¿½kna max och min vï¿½rde
max([f(x1) f(x2) f(x3) f(x4)])
min([f(x1) f(x2) f(x3) f(x4)])

plot(x1, f(x1), 'ro')
plot(x2, f(x2), 'ro')
plot(x3, f(x3), 'ro')
plot(x4, f(x4), 'ro')

%% 4 Integraler
xi = -3/2; 

a=1; b=2; f=@(x)x.^xi;
n=20;
dx=(b-a)/n;
kontroll = integral(f, a, b) 




q=0;
x=linspace(a,b,n+1);
qv=sum(f(x(1:n)).*dx)

subplot(2,2,1)
plot(x, f(x), 'r-'), hold on, grid on
for i=1:n
    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');
    

q=0;
x=linspace(a,b,n+1);
qh=sum(f(x(2:n+1)).*dx)

subplot(2,2,2)
plot(x, f(x), 'r-'), hold on, grid on
for i=1:n
    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');

q=0;
x=linspace(a,b,n+1);
qm=sum(f((x(1:n)+x(2:n+1))./2).*dx)

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


q=0;
x=linspace(a,b,n+1);
qt=sum(((f(x(1:n))+f(x(2:n+1)))./2).*dx)
%title('Trapetzmetoden')

%% Uppgift 5
L=0.1;g=9.80665; % g=gravity

c=0;d=pi/2; % Integrationsgrï¿½nserna
h=@(t,t0) 1./sqrt(1-sin(t0/2).^2.*sin(t).^2) % Funktionen 

a=deg2rad(10);b=deg2rad(170);n=160; % Intervallgrï¿½nser fï¿½r grafen, antal pkter.
x=linspace(a,b,n); % skapa grafvektor...
f=zeros(size(x))    % och skapa en funktionsvektor i samma lï¿½ngd, fylld med nollor.
for i=1:length(x) % Fï¿½r alla x (dvs. 10-170grader) hitta en lï¿½sning till funktionen.
    f(i)=4*sqrt(L/g)*integral(@(t)h(t,x(i)),c,d)
end
plot(rad2deg(x), f)
xlabel('\theta_0');ylabel('T(\theta_0)');
grid on

%% 5 Derivator
u=@(x)exp(x);
Du=@(x)exp(x);
h=logspace(-15,-1,400);
x=0.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((u(x+h)-u(x))./h))
loglog(h,abs((u(x+h)-u(x-h))./(2*h)))
loglog(h,abs((-u(x+2*h)+8*u(x+h)-8*u(x-h)+u(x-2*h))./(12*h)))


axis square, grid on
ylabel('|D_nu(x)-u''(x)|');xlabel('h')

%% Uppgift 6
tab=load('labtabell');
h=10e-5 % val av h baserad pï¿½ bï¿½sta nogrannhet for O(h^2) visad i uppgiften ovan.
x=1:length(tab.x)

% %kontrollexempel
% u=@(x)x.^3;
u=@(x) tab.y(x)
du=@(x,h)((u(x+h)-u(x-h))./(2*h)) % centraldifferenskvoten D0u(x)
ddu=@(x,h) (u(x+h)-2*u(x)+u(x-h))./h^2 % D+D-u(x)

dx=tab.x(2)-tab.x(1); % equidistant intervals; can be condensed to a single constant.
dy=diff(tab.y);
x=tab.x;y=tab.y;
df = dy./dx;
ddf = diff(df)./dx;

plot(x, tab.y), hold on, grid on
plot(x(2:end)-dx/2, df, 'r--')
plot(x(3:end)-dx, ddf, 'g-.')

plot(tab.x(2:end-1)-dx/2, ddf./(1+df(2:end).^2).^(3/2))

legend('f(x)', 'f''(x)', 'f''''(x)', '\kappa(x)')



%% Annat sï¿½tt at direkt rita derivata utifrï¿½n datapunkter.
% Pilfered from https://se.mathworks.com/matlabcentral/answers/72828-how-do-i-take-the-derivative-of-my-plot
u=@(x)x.^3;
x=linspace(-1,1,100);
dy = diff(u(x))./diff(x)
plot(x, u(x), x(2:end)-diff(x)./2, dy), grid on, hold on
% ^ diff() generates vector of length n-1, therefore shift x by half an
% intervall to compensate.
legend('show')