% EX 6.4

% a
% Beräkna I_1 med trapetsregeln + extrapolation
f = @(x) 30./(1 + x.^4 + sqrt(1 + x.^3));

% Bara två siffror: Kan göra för hand! Men vi gör det med dator ändå för
% att öva.

% Grövsta indelningen av intervallet [0 1]
h = 1;
x = 0:h:1;

% Plotta f
figure(1)
clf
hold on
patch([x x(end) x(1)], [f(x) 0 0], [0.8 0.8 0.8])
plot(x, f(x), 'ko')
xlabel('x');
ylabel('f(x)');

% Beräkna med olika stora steg h för att uppskatta felet
disp('    h         T         S      trunkfel');
for steg = 1:5
	T(steg) = h*(0.5*f(x(1)) + sum(f(x(2:(end-1)))) + 0.5*f(x(end)));
	if steg == 1
		S(steg) = T(steg);
		trunkfel = 1;
	else
		S(steg) = T(steg) + (T(steg) - T(steg - 1))/3;
		trunkfel = (T(steg) - T(steg - 1))/3;
	end
	disp(sprintf('%0.6f  %0.6f  %0.6f   %0.6f   %0.8f', ...
		h, T(steg), S(steg), trunkfel))
	h = h/2;
	x = 0:h:1;
end

  
% Körning:
%    h         T         S      trunkfel
%1.000000  11.893398  11.893398   1.000000   
%0.500000  13.011640  13.384387   0.372747   
%0.250000  13.223472  13.294083   0.070611   
%0.125000  13.274844  13.291968   0.017124   
%0.062500  13.287594  13.291844   0.004250   
% Dvs integralen = 13 +- 0.4 vilket fås med h = 0.5


% b
% Beräkna integral ?_0^inf f med 5 korrekta decimaler
% Oändlig integral = kapa svansen
%  ?_0^inf f =  ?_0^A f +  ?_A^inf f
% Fem korrekta decimaler ger int_A^inf < 5e-6
% int_A^inf < ?_A^inf 30/x^4 = 10/A^3 < 5e-6  =>  A^3 > 2e6 =>  A = 200 tex


% Gräns A
A = 200;

% Beräkna integral med quadl, ta till lite i osäkerhetskraven
Q = quadl(f, 0, A, 1e-9)

% Körning:
%Q = 18.590863466463780
% dvs Q = 18.59086 med 5 decimalers noggrannhet
	

	