% EX 2.10

% Ekvationen f(x) = ax - sin(x) - ln(1 + x) = 0, a = 1.96 har rot x = 0, 
% och också en andra liten, positiv rot. 

% Plotta
a = 1.96;
x = 0:0.001:0.1;
plot(x, a*x - sin(x) - log(1 + x))
grid on

% Använd serieutveckling (Taylorutveckling av varje term runt 0) för att 
% hitta roten > 0 med 1 siffras noggrannhet.
% sin(x) = 0 + 1x/1! - 0x^2/2! - 1x^3/3! + ... 
% ln(1 + x) = 0 + 1x/1! - 1x^2/2! + 2x^3/3! - ...
% f_3(x) = ax - (x - x^3/6) - (x - x^2/2 + x^3/3) = x((a - 2) + x/2 - x^2/6)
% x = 0 trivial rot
% Två termer: 0.04 - x/2 = 0                 =>  x = 0.08
% Tre termer: 0.04 - x/2 + x^2/6 = 0 
%             x^2 - 3x + 0.24 = 0 
%             x = 1.5 +- sqrt(1.5^2 - 0.24)  =>  x = 0.0823
% En siffras noggrannhet => x = 0.08


% Beräkna nu roten med minst 5 siffors noggrannhet
disp(' ')
disp('Exakt a:')

% Startpunkt för positiv rot
x_start = 0.08;

% Maximalt relativfel
h_max = 5e-6;

% Skapa funktioner f och f'
f = @(x) a*x - sin(x) - log(1 + x);
f_prim = @(x) a - cos(x) - 1/(1 + x);
[x, h] = newton_raphson(x_start, h_max, f, f_prim);
	

% Hur noggrannt kan roten anges om a inte är exakt, men korrekt avrundat?
% Störningsräkning
% a med två siffrors noggrannhet => a = 1.96 +- 0.005
% Instängning: lös ekvationen med a = 1.955 och a = 1.965
disp(' ')
disp('Störningsräkning på a:')

% Min
a = 1.955;
f = @(x) a*x - sin(x) - log(1 + x);
f_prim = @(x) a - cos(x) - 1/(1 + x);
x_min = newton_raphson(x_start, h_max, f, f_prim);

% Max
a = 1.965;
f = @(x) a*x - sin(x) - log(1 + x);
f_prim = @(x) a - cos(x) - 1/(1 + x);
x_max = newton_raphson(x_start, h_max, f, f_prim);

disp(sprintf('Om a = 1.96 +- 0.05, så är x = %0.4f+-%0.4f', ...
	(x_min + x_max)/2, abs(x_min - x_max)/2));


% Hur noggrant måste man ange a om roten ska ha fem gällande siffror?
% Felanalys med felfortplantingsformeln
% absolutfel e_a < |da/dx| e_x
% 5 gällande siffror ger absolutfel e_x < 5e-7, eftersom x = 0.082
% Implicit derivering: om f(x,a) = 0, så gäller da/dx = -(df/dx)/(df/da)
% da/dx = -(a - cos(x) - 1/(1 + x))/x
% da/dx(a = 1.96, x = 0.082) = -0.48
% e_a < 0.48*5e-7 = 2.4e-07    
%     =>  a bör ha 7 decimaler <=> 8 siffors precision
