function [approx , err_abs] = newton_ND_avec_der(F , mat_jac , x0 , nb_it_max , tol_rel , file_name)
% NEWTON_ND_AVEC_DER	Methode de Newton pour la resolution de F(x) = 0, pour F: R^n -> R^n
%
% Syntaxe: [approx , err_abs] = newton_ND_avec_der(F , mat_jac , x0 , nb_it_max , tol_rel , file_name)
%
% Arguments d'entree
%	F			-	String ou fonction handle specifiant la fonction
%					non-lineaire (F: R^n -> R^n)
%	mat_jac		-	String ou function handle specifiant la matrice
%					jacobienne de F (Jac: R^n -> R^{n x n})
%	x0			-	Approximation initiale (x0 dans R^n)
%	nb_it_max	-	Nombre maximum d'iterations 
%	tol_rel		-	Tolerance sur l'approximation de l'erreur relative
%	file_name	-	(Optionnel) Nom du fichier (avec l'extension .txt) dans
%					lequel sera	ecrit les resultats de l'algorithme
%
% Arguments de sortie
%	approx		-	Matrice de taille (n x nb_iter) contenant les 
%					iterations
%	err_abs		-	Vecteur rangee de dimension nb_iter contenant les
%					erreurs absolues
%
% Exemples d'appel
%	[ approx , err_abs ] = newton_ND_avec_der( 'my_sys_nl' , 'my_sys_nl_jac' , [1,1] , 20 , 1e-9 , 'resul_newtonND.txt')
%	[ approx , err_abs ] = newton_ND_avec_der( @(x) [x(1)^2 + x(2)^2 - 1 ; -x(1)^2 + x(2)] , @(x) [2*x(1),2*x(2);-2*x(1),1] , [1,1] , 20 , 1e-9 , 'resul_newtonND.txt')



%%  Verification de la fonction F et de la matrice jacobienne
if isa(F,'char')
	fct			=	str2func(F);
	is_fct_file =	true;
elseif isa(F,'function_handle')
	fct			=	F;
	is_fct_file =	false;
else
	error('L''argument f n''est pas un string ni un function_handle')
end

if isa(mat_jac,'char')
	jac				=	str2func(mat_jac);
	is_jac_file		=	true;
elseif isa(mat_jac,'function_handle')
	jac				=	mat_jac;
	is_jac_file		=	false;
else
	error('L''argument df n''est pas un string ni un function_handle')
end

%% Verification du nb de composantes des approximations initiales et de f
if ~isnumeric(x0) || ~isvector(x0)
	error('L''approximation initiale x0 n''est pas un vecteur')
end

x0_col	=	reshape(x0,[],1);

try 
	fct(x0_col);
catch ME
	rethrow(ME)
end

try 
	jac(x0_col);
catch ME
	rethrow(ME)
end

taille	=	length(x0_col);

if ~isnumeric(fct(x0_col)) || ~isvector(fct(x0_col)) || (length(fct(x0_col))~=taille)
	error(['Le fonction F ne retourne pas un vecteur de meme taille que x0. ',...
		'x0 est de taille %d alors que F(x0) est de taille %d.'],taille,length(fct(x0)))
elseif ~isnumeric(jac(x0_col)) || ~ismatrix(jac(x0_col)) || ~isequal(size(jac(x0_col)),[taille,taille]) 
	error(['Le fonction mat_jac ne retourne pas une matrice de bonne taille; ',...
		   'x0 est de taille %d, alors que mat_jac(x0) est de dimension %d x %d.'],...
										taille,size(jac(x0_col),1),size(jac(x0_col),2))
elseif ~check_derivative(fct,jac,x0_col)
	warning('Il semble y avoir une erreur avec la matrice jacobienne')
end

%% Verification du fichier output
if nargin == 6 && ~isa(file_name,'char')
	error('Le nom du fichier des resultats doit etre de type string')
end

%% Initialisation des matrices app et err
app			=	nan(taille,nb_it_max);
app(:,1)	=	x0_col;
err_rel		=	inf(1,nb_it_max);
arret		=	false;


%% Methode de Newton
for t=1:nb_it_max-1
	
	jacobienne	=	jac(app(:,t));
	delta_x		=	jacobienne\-reshape(fct(app(:,t)),taille,1);
	app(:,t+1)	=	app(:,t) + delta_x;
	
	if any(~isfinite(jacobienne(:)))
		warning(['La matrice jacobienne de f a l''iteration %d est singuliere 0.\n',...
					'Arret de l''algorithme'],t)
		break
	end
	
	if min(abs(eig(jacobienne))) == 0  
		warning(['La matrice jacobienne de f a l''iteration %d est singuliere 0.\n',...
					'Arret de l''algorithme'],t)
		break
	end
	
	
	err_rel(t)	=	norm(app(:,t+1)-app(:,t))/(norm(abs(app(:,t+1))) + eps);
	if (err_rel(t) <= tol_rel) || (norm(fct(app(:,t+1))) == 0)
		arret	=	true;
		break
	end
end

nb_it	=	t+1;
approx	=	app(:,1:nb_it);
err_abs	=	inf(1,nb_it);

if arret
	for t=1:nb_it
		err_abs(t)	=	norm(approx(:,end) - approx(:,t));
	end
else
	warning('La methode de Newton n''a pas convergee')
end

% ecriture des resultats si fichier passe en argument
if nargin == 6
	output_results(file_name , fct , is_fct_file , jac, is_jac_file, ...
				nb_it_max , tol_rel , x0 , approx , err_abs , arret)
end

end


function [test] = check_derivative(f,jac,x0)

	taille	=	length(x0);
	if min(x0) == 0
		h_init	=	1e-6;
	else
		h_init	=	1e-3 * min(x0);
	end	
	h		=	h_init./(2.^(0:4));
	erreur	=	nan(size(h));

	
	for t=1:length(h)
		app			=	zeros(taille,taille);
		for d=1:taille
			delta_h		=	zeros(size(x0));
			delta_h(d)	=	h(t);
			app(:,d)	=	(f(x0+delta_h) - f(x0-delta_h))/(2*h(t));
		end
		erreur(t)	=	norm(jac(x0) - app);
	end
	
	ordre	=	log(erreur(2:end) ./ erreur(1:end-1)) ./ ...
				log(h(2:end)      ./ h(1:end-1));
			
	test	=	(mean( abs(2 - ordre) ) <= 0.25) || (mean(erreur/(norm(jac(x0))+eps)) <=1e-9);	
	
end

function [] = output_results(file_name , fct , is_fct_file , jac, ...
			is_jac_file, it_max , tol_rel , x0 , x , err , status)
						 
	fid		=	fopen(file_name,'w');
	fprintf(fid,'Algorithme de Newton avec la matrice jacobienne\n\n');
	fprintf(fid,'Fonction dont on cherche les racines:\n');
	if is_fct_file
		fprintf(fid,'%s\n\n',fileread([func2str(fct),'.m']));
	else
		fprintf(fid,'%s\n\n',func2str(fct));
	end
	fprintf(fid,'Matrice jacobienne de la fonction dont on cherche les racines:\n');
	if is_jac_file
		fprintf(fid,'%s\n\n',fileread([func2str(jac),'.m']));
	else
		fprintf(fid,'%s\n\n',func2str(jac));
	end
	fprintf(fid,'Arguments d''entree:\n');
	fprintf(fid,'    - Nombre maximum d''iterations: %d\n',it_max);
	fprintf(fid,'    - Tolerance relative: %6.5e\n',tol_rel);
	fprintf(fid,'    - Approximation initiale x0: [');
	[taille, nb_iter]	=	size(x);
	if taille <= 3
		fprintf(fid,'%6.5e ',x0(1:taille));
		fprintf(fid,']\n');
	else
		fprintf(fid,'%6.5e ',x0(1:2));
		fprintf(fid,'... %6.5e]',x0(end));
	end
	
	if status
		fprintf(fid,'\nStatut: L''algorithme de Newton a converge en %d iterations\n\n',nb_iter-1);
	else
		fprintf(fid,'\nStatut: L''algorithme de Newton n''a pas converge\n\n');
	end
	if taille == 1
		fprintf(fid,'#It       x_1           Erreur absolue\n');
		fprintf(fid,'%3d   %16.15e   %6.5e\n',[reshape(0:nb_iter-1,1,[]);reshape(x(:),1,[]);reshape(err,1,[])]);
	elseif taille == 2
		fprintf(fid,'#It       x_1           x_2       Erreur absolue\n');
		fprintf(fid,'%3d   %6.5e   %6.5e   %6.5e\n',[reshape(0:nb_iter-1,1,[]);reshape(x([1,2],:),2,[]);reshape(err,1,[])]);
	elseif taille == 3
		fprintf(fid,'#It       x_1           x_2           x_3       Erreur absolue\n');
		fprintf(fid,'%3d   %6.5e   %6.5e   %6.5e   %6.5e\n',[reshape(0:nb_iter-1,1,[]);reshape(x([1,2,3],:),3,[]);reshape(err,1,[])]);
	else
		fprintf(fid,'#It       x_1           x_2      ...       x_n       Erreur absolue\n');
		fprintf(fid,'%3d   %6.5e   %6.5e  ...  %6.5e   %6.5e\n',[reshape(0:nb_iter-1,1,[]);reshape(x([1,2,taille],:),3,[]);reshape(err,1,[])]);
	end
	
	fclose(fid);

end