clear 
close all

% dati del problema
a = input('inserisci numero ');

% iterato iniziale
b = input('inserisci approssimazione iniziale ');

% criteri di arresto
tol = input('inserisci la tolleranza per il criterio di arresto ');
kmax = 100;

[x,k,bvec,hvec] = erone(a,b,kmax,tol);
fprintf('\nRadice calcolata: %4.16f\n',x);
fprintf('Numero di iterazioni=%d \n', k);
fprintf('Errore=%d \n', abs(x-sqrt(a)));

figure; plot(1:k,bvec,'*r',1:k,hvec,'*b',1:k,sqrt(a)*ones(size(bvec)),'b-'); 
xlabel('Iterazioni');
legend('Approssimazioni b','Approssimazioni h','valore esatto')
title('Approssimazioni');

err_rel = abs(bvec-sqrt(a))/abs(b-sqrt(a));
kk=0:length(err_rel)-1;
figure; plot(kk,err_rel,'*-'); hold on; 
plot(kk,1./(2.^kk)); 
hold off
legend('|x_{k+1}-sqrt(a)|/|x_k-sqrt(a)|','1/2^k');
title('andamento dell''errore');

figure; plot(abs(diff(bvec)),'-*'); title('|x_{k+1}-x_k|');



%--------------------------------------------------------------------------
function [x,k,bvec,hvec] = erone(a,b,kmax,tol)
% erone: calcola la radice quadrata di un numero
% INPUT: a = numero del quale si vuole determinare la radice quadrata 
%        b = approssimazione iniziale
%        kmax = massimo numero di iterazioni
%        tol = tolleranza per il criterio di arresto sulla differenza tra
%        due approssimazioni successive
% OUTPUT: x = approssimazione calcolata
%         k = numero di iterazioni
%         bvec = vettore delle approssimazioni successive


if a < 0
    disp('Errore: il numero deve essere positivo');
    x=[]; k=[]; bvec=[];
else
    % Inizializza 
    bvec = zeros(kmax+1,1);
    hvec = zeros(kmax+1,1);
    bvec(1) = b;
    hvec(1) = a/b;
    fprintf('k=1, b=%4.8f, h=%4.8f \n', bvec(1),hvec(1));
    
    continua = true; k=1;
    while continua 

        % aggiorna la stima
        bvec(k+1) = 0.5*(bvec(k)+hvec(k));
        hvec(k+1) = a/bvec(k+1);
        
        % visualizza le stime calcolate
        fprintf('k=%d, b=%4.16f, h=%4.16f, err=%4.16f \n', k+1, bvec(k+1),hvec(k+1),abs(bvec(k+1)-sqrt(a)));

        % controlla condizioni di arresto 
        diff = abs(bvec(k+1)-bvec(k));
        continua = k<kmax & diff >tol;
        
        k = k+1;

    end
    bvec = bvec(1:k);
    hvec = hvec(1:k);
    x = bvec(end);
end
end