%% TMA671 Lösning av uppgift N.4.16
% Under storgruppsövningen som hölls onsdagen den 24e April 2018
% så demonstrerades uppgift 16 i Kapitel 4 i Numerisk Analys
% boken. Uppgiften handlar om att approximera en integral med hjälp
% av två olika kvadraturregler
clc,clear;
% Definiera det exakta värdet
exactValue = log(2) + (pi/2) -2;
% Utskrift till användaren
fprintf('========================================================================================\n\n');
fprintf('\nUPPGIFT N.4.16\n\n');
fprintf('\tApproximera integralen $I=\\int_{0}^{1}\\mathrm{ln}(1+x^2)\\mathrm{d}x$\n');
fprintf('\tvars värde är\n\t\t''''log(2)+ pi/2 - 2''''\t=\t%0.5f\n\n',exactValue);
fprintf('========================================================================================\n\n');
%% Delupggift (a)
% Beräkna integralen \int_{0}^{1}ln(1+x^2)dx med MATLABs
% ''quad'' and ''quadl''. Jämför nogrannheet och tidsåtgång
%-----------------------------------------------------------------------
% Skriv ut till användaren
fprintf('\nDELUPPGIFT (a)\n');
% Okajjj... så vi definierar ett funktionshandtag
f = @(x) log(1+x.^2);
% Tidsåtgång quad 
tStart = tic;
% Beräkna integral med ''quad''
IQuad = quad(f,0,1);
% Sluta ta tid
tQuad = toc(tStart);
% Beräkna felet för quad
felQuad = abs(IQuad - exactValue);
% Skriv ut svaret
fprintf('\tFör quad har vi:\n\t\tTid\t=\t%0.6f millisekunder\t\n\t\tFel\t=\t%0.3f*10^(-7)\n',tQuad*10^3,felQuad*10^7);
% Samma grej för quadl
% Tidsåtgång quadl 
tStart = tic;
% Beräkna integral med ''quad''
IQuadl = quadl(f,0,1);
% Sluta ta tid
tQuadl = toc(tStart);
% Beräkna felet för quad
felQuadl = abs(IQuadl - exactValue);
% Skriv ut svaret
fprintf('\tFör quadl har vi:\n\t\tTid\t=\t%0.6f millisekunder\t\n\t\tFel\t=\t%0.3f*10^(-7)\n',tQuadl*10^3,felQuadl*10^7);
%-----------------------------------------------------------------------
%% Delupggift (b)
% Beräkna integralen approximativt med trapetsmetoden, dvs. med trapz
% i MATLAB. Välj olika steglängder h och jämför resultatet med teorin
fprintf('\nDELUPPGIFT (b)\n');
fprintf('\n\tVi kollar nu hur steglängden h påverkar felet för trapetsmetoden...\n\n');
% Steglängder
h = [0.1; 0.01; 0.001; 0.0001];
% Integraler
I = zeros(size(h));
% Fel
Fel = zeros(size(h));
%  Vi skriver ut lite grejor
fprintf('\t\t\t--------------------------------------------------------------\n\n');
fprintf('\t\t\t\th\t\tFel*10^(-5)\t\tFel/h^2\n\n');
fprintf('\t\t\t--------------------------------------------------------------\n');
% Vi loopar och beräknar integralen
for i = 1:length(h)
   % Vi definierar en variabelvektor
   xv = (0:h(i):1)';
   % Vi definierar funktionen vi vill integrera
   y = f(xv);
   % Vi beräknar integralen
   I(i) = trapz(xv,y);    
   % Fel för integraler
   Fel(i) = abs(I(i)-exactValue).*10^5;
   % Vi skriver ut
   fprintf('\t\t\t\t%0.4f\t\t%0.10f\t\t%0.5f\n\n',h(i),Fel(i),Fel(i)/(h(i)^2));
end

fprintf('\t\t\t-------------------------------------------------------------\n');

fprintf('========================================================================================\n\n');