// Analyse numérique : TP n°1
// - Calcul des zéros de fonctions
// © EISTI 2013-2014 - Semestre 2

// Auteurs : INGOUFF Christian, TYNDAL Pierre-Alexandre
// Groupe : C06
// Date de création : 10/03/2014
// Dernière modification : 23/03/2014

//
// ################################################# Liste de fonctions à analyser
//

// 5a. f(x) = x^2 - a, a > 0 (Calcul de la racine carrée de a)
// Pré-conditions :
//		x : réel
//		a : réel positif
function y = f1(a, x)
	if (a > 0)
		y = x*x - a;
	else
		disp("Erreur : a est négatif");
		y = %nan;
	end
endfunction

// 5b. f(x) = x^2 + 17 (Point initial : x0 = 1)
// Pré-conditions :
//		x : réel
function y = f2(x)
	y = x*x + 17;
endfunction

// 5c. f(x) = x - 2 sin(x) (Point initial : x0 = 1.1)
// Pré-conditions :
// 		x : réel
function y = f3(x)
	y = x - 2*sin(x);
endfunction

// 8a. f(x) = sin(3x) + x^2 - 2 (Point initial : x0 = 1)
// Pré-conditions :
//		x : réel dans [-3, 3]
function y = f4(x)
	if (x < -3)
		y = -%inf;
	elseif (x > 3)
		y = %inf;
	else
		y = sin(3*x) + x*x - 2;
	end
endfunction

// 8b. f(x) = x cos(x) (Point initial : x0 = 0.1 et 2.5)
// Pré-conditions :
//		x : réel dans [-2pi, 2pi]

function y = f5(x)
	if (x < -2*%pi)
		y = -%inf;
	elseif (x > 2*%pi)
		y = %inf;
	else
		y = x*cos(x);
	end
endfunction

//
// ################################################# Dérivées analytiques des fonctions
//

function y = df1(x)
	y = 2*x;
endfunction

function y = df2(x)
	y = 2*x;
endfunction

function y = df3(x)
	y = 1 - 2*cos(x);
endfunction

function y = df4(x)
	if (x < -3)
		y = -%inf;
	elseif (x > 3)
		y = %inf;
	else
		y = 3*cos(3*x) + 2*x;
	end
endfunction

function y = df5(x)
	if (x < -2*%pi)
		y = -%inf;
	elseif (x > 2*%pi)
		y = %inf;
	else
		y = cos(x) - x*sin(x);
	end
endfunction

// Dérivées secondes

function y = ddf4(x)
	y = -9*cos(3*x) + 2;
endfunction

function y = ddf5(x)
	y = - 2*sin(x) - x*cos(x);
endfunction

//
// ################################################# Dérivée approchée
//

// Pré-conditions :
//		f : chaîne de caractère dénotant la fonction choisie (de la forme "f1(a," ou "fi(")
//		x : réel
//		h : réel strictement positif, proche de 0
//		isDouble : boolean, vrai si en mode dérivée seconde
function y = deriv(f, x, h, isDouble)
	if (isDouble)
		fxh = deriv(f, x+h, h, %f);
		fx = deriv(f, x, h, %f);
		y = (fxh - fx) / h;
	else
		y = (eval(f + string(x) + "+" + string(h) + ")") - eval(f + string(x) + ")")) / h;
	end
endfunction

//
// ################################################# Utilitaires
//

// Affichage de la fonction avec la racine représentée
// Pré-conditions :
//		f : chaîne de caractère dénotant la fonction choisie (de la forme "f1(a," ou "fi(")
//		x0 : réel, point initial
//		racine : réel, racine trouvée
function res = plotRoot(f, x0, racine)
	clf();
	X = (racine - 10):0.01:(racine + 10);// axe des x
	l = length(X);
	for i = 1:l
		Y(i) = eval(f + string(X(i)) + ")");// axe des y
	end
	xset("window", 0);
	xtitle("Représentation d''une racine pour " + part(f,1:strindex(f,"(")-1), "x", f + "x)");
	plot([(racine - 10) (racine + 10)]', [0 0]', 'g'); // axe des abscisses
	plot(X,Y,'r'); // Représentation de la fonction
	plot([racine racine]',[-1 1]','b'); // Marquage de la racine
endfunction

// Affichage de la fonction avec le minimum représenté
// Pré-conditions :
//		f : chaîne de caractère dénotant la fonction choisie (de la forme "f1(a," ou "fi(")
//		x0 : réel, point initial
//		mini : réel, x minimum trouvé
function res = plotMini(f, x0, mini)
	clf();
	if (grep(f,"f4") == 1)
		bounds = 3;
	else
		bounds = 2*%pi;
	end
	minValue = eval(f + string(mini) + ")");// Valeur du minimum
	X = (-bounds):0.01:(bounds);// axe des x
	l = length(X);
	for i = 1:l
		Y(i) = eval(f + string(X(i)) + ")");// axe des y
	end
	xset("window", 0);
	xtitle ("Représentation du minimum pour " + part(f,1:strindex(f,"(")-1), "x", f+"x)");
	plot([-bounds bounds]', [minValue minValue]', 'b'); // Axe horizontal du minimum
	plot(X,Y,'r'); // Représentation de la fonction
	plot([mini mini]',[minValue-1 minValue+1]','b'); // Marquage du minimum
endfunction

// Demande oui ou non
function res = yesNoInput()
	res = -1;
	while (res < 0)
		choice = input("(o/n) : ","string");
		choice = part(choice,1);
		if (choice == "o")
			res = 1;
		elseif (choice == "n")
			res = 0;
		end
	end
endfunction

//
// ################################################# Méthode Newton-Raphson
//

// Application du domaine de définition pour f4 et f5
// Pré-conditions :
//		f : chaîne de caractère dénotant la fonction choisie (de la forme "f1(a," ou "fi(")
//		x : réel
function x1 = enforceBounds(f, x)
	if (grep(f,"f4") == 1)
		if (x < -3)
			x1 = -3;
		elseif (x > 3)
			x1 = 3;
		else
			x1 = x;
		end
	elseif (grep(f,"f5") == 1)
		if (x < -2*%pi)
			x1 = -2*%pi;
		elseif (x > 2*%pi)
			x1 = 2*%pi;
		else
			x1 = x;
		end
	end
endfunction

// Intersection de la droite tangente de f au point f(x0)
// Pré-conditions :
//		f : chaîne de caractère dénotant la fonction choisie (de la forme "f1(a," ou "fi(")
//		x0 : réel (point initial)
function x1 = interAnalytic(f, x0)
	if (grep(f,"f1") == 1)
		df = "d"+part(f,1:strindex(f,"("));
	else
		df = "d"+f;
	end

	if (eval(df + string(x0) + ")") <> 0)
		x1 = x0 - eval(f + string(x0) + ")") / eval(df + string(x0) + ")");
	else
		disp("Erreur : la dérivée est égale à 0");
		x1 = x0;
	end
	
	x1 = enforceBounds(f, x1);
endfunction

// Intersection : Méthode approximation Taylor
// Pré-conditions :
//		f : chaîne de caractère dénotant la fonction choisie (de la forme "f1(a," ou "fi(")
//		x0 : réel (point initial)
function x1 = interApprox(f, x0)
	h = 0.0001;
	df = deriv(f, x0, h, %f);

	if (grep(f,"/f[45]/",'r') == 1)
		ddf = deriv(f, df, h, %t);
		if (ddf <> 0)
			x1 = x0 - df / ddf;
		else
			disp("Erreur : la dérivée est égale à 0");
			x1 = x0;
		end
	else
		if (df <> 0)
			x1 = x0 - eval(f + string(x0) + ")") / df;
		else
			disp("Erreur : la dérivée est égale à 0");
			x1 = x0;
		end
	end
	
	x1 = enforceBounds(f, x1);
endfunction

// Pré-conditions :
//		f : chaîne de caractère dénotant la fonction choisie (de la forme "f1(a," ou "fi(")
//		x0 : réel (point initial)
//		precision : réel strictement positif (précision), proche de 0
//		isAnalytic : booléen, vrai si mode analytique, faux sinon
// Post-conditions : Le but est de calculer la racine de fonctions non linéaires réelles
function res = NewtonRaphson(f, x0, precision, isAnalytic)
	i = 1; // Compteur d'itérations
	max_count = 10000; // Limite supérieure du compteur
	
	X(1) = x0; // Tableau des x au fur et à mesure du calcul des intersections
	ecart = %inf; // Ecart entre X(n) et X(n-1)
	
	while ((ecart > precision) & (i < max_count))
		if (isAnalytic)
			X(i+1) = interAnalytic(f, X(i));
		else
			X(i+1) = interApprox(f, X(i));
		end	
		ecart = abs(X(i+1) - X(i));
		i = i+1;
	end
	
	fprintfMat("log_iteration.txt", X, "%10.10f"); // Enregistre les valeurs du tableau X dans un fichier texte
	disp("Nombre d''itérations : " + string(i));
	disp("La liste des valeurs itérées est disponible dans log_iteration.txt");
	
	if (grep(f,"/f[123]/","r") == 1)
		if (ecart < precision)
			disp("Racine : x = " + string(X(i)));
			disp("Représenter la solution graphiquement ?");
			choice = yesNoInput();
			if (choice == 1)
				plotRoot(f, x0, X(i));
			end
		else
			disp("Cette fonction n''a pas de racine");
		end
	else
		if (ecart < precision)
			disp("Minimum : x = " + string(X(i)));
			disp("Représenter la solution graphiquement ?");
			choice = yesNoInput();
			if (choice == 1)
				plotMini(part(f,strindex(f,"f"):length(f)), x0, X(i));
			end
		else
			disp("Le minimum n''a pas pu être trouvé");
		end			
	end
	
	res = 0;
endfunction

//
// ################################################# Fonctions : entrées utilisateur
//

function res = pickFunction(choice,precision)	
	if (choice >= 1) & (choice <= 6)
		disp("Utiliser la formule analytique de la dérivée pour les racines ?");
		disp("(l''approximation de Taylor au premier ordre sera utilisé le cas échéant)");
		analytic = yesNoInput();
	end
	
	select choice,
		case 0 then res = 0;,
		case 1 then
			a = -1;
			while (a < 0)
				disp("Spécifier l''entier a");
				a = input("a = ");
				if (a < 0)
					disp("Erreur : le a entré est négatif");
				end
			end
			disp("Spécifier le point initial x0");
			x0 = input("x0 = ");
			disp("f(x) = x^2 - " + string(a) + " (Point initial : x0 = " + string(x0) + ")");
			NewtonRaphson("f1("+string(a)+",", x0, precision, analytic == 1);
		,
		case 2 then
			disp("f(x) = x^2 + 17 (Point initial : x0 = 1)");
			NewtonRaphson("f2(", 1, precision, analytic == 1);
		,
		case 3 then 
			disp("f(x) = x - 2 sin(x) (Point initial : x0 = 1.1)");
			NewtonRaphson("f3(", 1.1, precision, analytic == 1);
		,
		case 4 then
			disp("f(x) = sin(3x) + x^2 - 2 (Point initial : x0 = 1)");
			if (analytic == 1)
				NewtonRaphson("df4(", 1, precision, analytic == 1);
			else
				NewtonRaphson("f4(", 1, precision, analytic == 1);
			end
		,
		case 5 then
			disp("f(x) = x cos(x) (Point initial : x0 = 0.1)");
			if (analytic == 1)
				NewtonRaphson("df5(", 0.1, precision, analytic == 1);
			else
				NewtonRaphson("f5(", 0.1, precision, analytic == 1);
			end
		,
		case 6 then
			disp("f(x) = x cos(x) (Point initial : x0 = 2.5)");
			if (analytic == 1)
				NewtonRaphson("df5(", 2.5, precision, analytic == 1);
			else
				NewtonRaphson("f5(", 2.5, precision, analytic == 1);
			end
		,
		case 7 then res = 2;,
		else disp("Veuillez entrer un entier entre 0 et 7"); res = 1;,
	end
	
	if (choice >= 1) & (choice <= 6)
		disp("Relancer le programme ?");
		res = yesNoInput();
	end
endfunction

function res = main()
	precision = 0.000001; // Précision de la méthode Newton-Raphson
	res = 1;
	while (res > 0)
		disp("Veuillez entrer l''entier correspondant à la fonction :");
		disp("Calcul de racines :");
		disp("1. f(x) = x^2 - a, a > 0");
		disp("2. f(x) = x^2 + 17 (Point initial : x0 = 1)");
		disp("3. f(x) = x - 2 sin(x) (Point initial : x0 = 1.1)");
		disp("Calcul de minimums :");
		disp("4. f(x) = sin(3x) + x^2 - 2 (Point initial : x0 = 1)");
		disp("5. f(x) = x cos(x) (Point initial : x0 = 0.1)");
		disp("6. f(x) = x cos(x) (Point initial : x0 = 2.5)");
		disp("Programme :");
		disp("7. Régler la précision de la méthode Newton-Raphson (actuellement "+string(precision)+")");
		disp("0. Quitter le programme");
		choice = input("Réponse : ");
		res = pickFunction(choice,precision);
		if (res == 2)
			precision = input("Nouvelle précision :");
		end
	end
endfunction
