function [res] = tp2()

%% question 3
% question 3-2
figNum = 1;
[res]  = question32(figNum);

% question 3-3
figNum = 2;
[res]  = question33(figNum);

% question 3-4 
  theta  = 0;   figNum = 3;
  [res]  = question34CasUnTheta(theta, figNum);
  theta  = 0.5; figNum = 4;
  [res]  = question34CasUnTheta(theta, figNum);
  theta  = 1;   figNum = 5;
  [res]  = question34CasUnTheta(theta, figNum);
 
% question 3-5  
  theta  = 0;   figNum = 6;
  [res]  = question35CasUnTheta(theta, figNum);
  theta  = 0.5; figNum = 7;
  [res]  = question35CasUnTheta(theta, figNum);
  theta  = 1;   figNum = 8;
  [res]  = question35CasUnTheta(theta, figNum);

%% question 4
  figNum = 9;
  [res]  = question4(figNum);
res =1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%  question 3  %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%==============================================================
function [res]  = question32(figNum)
theta=[0,0.25,0.5,0.75,1];
col     = ['r';'*b';'+g';'m';'--k'];
tol     = 1e-3;
IterMax = 10;
t0      = 0;
tf      = pi;
x0      = [1;0];
N       = 50 ;

figure(figNum); clf; hold off;

% solutions approchees
for i=1:length(theta)
 monSchema = @(f, t0, x0, h0) monThetaSChemaNewton(f, t0, x0, h0, tol, IterMax,theta(i));
 [t,x,h]   = solveurEDO(@monEDO,t0,tf,x0,N,monSchema);
 plot(x(1,:), x(2,:),col(i,:)); hold on;
end
% solution exacte
x = cos(t); y = -sin(t);
plot(x,y,'k'); 

legend('theta = 0','theta = 0.25','theta =0.5','theta =0.75', 'theta = 1', 'sol exacte');
title('Solution avec methode de Newton');
hold off;

% impression ndu graphique
print('Q32SolutionNewton.eps','-depsc');
res =1;
end
%%==================================================================
function [res]  = question33(figNum)
theta=[0,0.25,0.5,0.75,1];
col     = ['r';'*b';'+g';'m';'--k'];
tol     = 1e-3;
IterMax = 10;
t0      = 0;
tf      = pi;
x0      = [1;0];
N       = 50 ;

figure(figNum); clf; hold off;

% solutions approchees
for i=1:length(theta)
 monSchema = @(f, t0, x0, h0) monThetaSChemaPicard(f, t0, x0, h0, tol, IterMax,theta(i));
 [t,x,h]   = solveurEDO(@monEDO,t0,tf,x0,N,monSchema);
 plot(x(1,:), x(2,:),col(i,:)); hold on;
end

% solution exacte
x = cos(t); y = -sin(t);
plot(x,y,'k'); 

legend('theta = 0','theta = 0.25','theta =0.5','theta =0.75', 'theta = 1', 'sol exacte');
title('Solution avec methode de Point fixe');
hold off;

% impression ndu graphique
print('Q33SolutionPointFixe.eps','-depsc');
res =1;
end
%%==================================================================
function [res]  = question34CasUnTheta(theta, figNum)
% theta  : valeur de theta
% figNum : numero de la figure
%----------------------------------------------------

tol       = 1e-3;
IterMax   = 10;
t0        = 0;
tf        = pi;
x0        = [1;0];
N         = 50 ;
monSchema = @(f, t0, x0, h0) monThetaSChemaNewton(f, t0, x0, h0, tol, IterMax,theta);

H = []; % pour stocker les pas 
E = []; % pour stocker les erreurs

for k=1:12
  N = 2^k ;
  [t,x,h] = solveurEDO(@monEDO,t0,tf,x0,N,monSchema);
  H = [H,h];
  E = [E, norm(x(:,end) - [-1;0])];
end

% calcul des log
logH = log(H);
logE = log(E);

% droites de pente 1 et 2 passant par le dernier point
x   = linspace(min(log(H)), max(log(H)), 30);
y1  = (x - log(H(end))) + log(E(end));
y2  = 2 * (x - log(H(end))) + log(E(end));

%
figure(figNum); clf; hold off;
plot(logH,logE,'-*'); hold on;     % courbe log log
plot(x,y1,'g',x,y2,'--r');         % droites de pente 1 et 2
xlabel('Log(h)'); ylabel('Log(E)');
legend('log(E)', strcat('droite de pente = ','1'), strcat('droite de pente =', '2'));
title(strcat('Courbe pour theta = ', num2str(theta)));

% Impression du graphique au format .eps couleur
nomFichier = strcat('Q34CourbeTheta', num2str(theta),'.eps');
print(nomFichier,'-depsc');
res=1;
end
%%==================================================================
function [res]  = question35CasUnTheta(theta, figNum)
% theta  : valeur de theta
% figNum : numero de la figure
%----------------------------------------------------

tol       = 1e-3;
IterMax   = 1;          % C'EST ICI  QU'IL Y A EU MODIFICATION
t0        = 0;
tf        = pi;
x0        = [1;0];

monSchema = @(f, t0, x0, h0) monThetaSChemaPicard(f, t0, x0, h0, tol, IterMax,theta);

H = []; % pour stocker les pas 
E = []; % pour stocker les erreurs

for k=1:12
  N = 2^k ;
  [t,x,h] = solveurEDO(@monEDO,t0,tf,x0,N,monSchema);
  H = [H,h];
  E = [E, norm(x(:,end) - [-1;0])];
end

% calcul des log
logH = log(H);
logE = log(E);

% droites de pente 1 et 2 passant par le dernier point
x   = linspace(min(log(H)), max(log(H)), 30);
y1  = (x - log(H(end))) + log(E(end));
y2  = 2 * (x - log(H(end))) + log(E(end));

%
figure(figNum); clf; hold off;
plot(logH,logE,'-*'); hold on;     % courbe log log
plot(x,y1,'g',x,y2,'--r');         % droites de pente 1 et 2
xlabel('Log(h)'); ylabel('Log(E)');
legend('log(E)', strcat('droite de pente = ','1'), strcat('droite de pente =', '2'));
title(strcat('Courbe pour theta = ', num2str(theta)));

% Impression du graphique au format .eps couleur
nomFichier = strcat('Q35CourbeTheta', num2str(theta),'.eps'); % MODIFICATION FAITE ICI AUSSI
print(nomFichier,'-depsc');
res=1;
end
%%%%%%%%%%%%%%%%%% monThetaSChemaNewton %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [t1, x1, h1] = monThetaSChemaNewton(f, t0, x0, h0, tol, IterMax,theta)
t1 = t0 + h0;
gPourNewton = @(Z) Ncursif(Z,f,t0,x0,h0,theta);
[x1,N] = Newton(gPourNewton,x0,tol,IterMax);
h1 = h0 ;
end
%------- fonction Ncursif ------
function [Nc, DzNc] = Ncursif(Z,f,t0,X0,h,theta)
[F0, DF0] = f(t0,X0);
[F1,   DF1] = f(t0+h,Z);
dim = length(Z);
Nc  = Z - X0  - h * ( (1-theta)* F0  + theta * F1 );
DzNc = eye(dim)  - h * theta * DF1;
end
%%%%%%%%%%%%%%%%%% monThetaSChemaPicard %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [t1, x1, h1] = monThetaSChemaPicard(f, t0, x0, h0, tol, IterMax, theta)
t1 = t0 + h0;
[F0, DF0] = f(t0, x0);
gPourPicard = @(Z) Pcursif(Z,f,t0,x0,h0,theta);  
[x1,N] = Picard(gPourPicard,x0,tol,IterMax);
h1 = h0 ;
end
%------- fonction Pcursif ------
function [Pc] = Pcursif(Z,f,t0,X0,h,theta)
F0  = f(t0,X0);

F1  = f(t0+h,Z);

Pc  = X0  + h * ( (1-theta)* F0  + theta * F1 );
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%% question 4 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [res] = question4(figNum)
t0 = 0;
tf = pi;
x0 =[1;0];
H  = [];
E1 = []; %Pour schema 1
E2 = []; %Pour schema 2
E3 = []; %Pour schema 3
E4 = []; %Pour schema 4

for k=1:12
  N = 2^k ;
  [t,x1,h] = solveurEDO(@monEDO,t0,tf,x0,N,@Schema1);
  H  = [H,h];
  E1 = [E1, norm(x1(:,end) - [-1;0])];
  [t,x2,h] = solveurEDO(@monEDO,t0,tf,x0,N,@Schema2);
  E2 = [E2, norm(x2(:,end) - [-1;0])];
  [t,x3,h] = solveurEDO(@monEDO,t0,tf,x0,N,@Schema3);
  E3 = [E3, norm(x3(:,end) - [-1;0])];
  [t,x4,h] = solveurEDO(@monEDO,t0,tf,x0,N,@Schema4);
  E4 = [E4, norm(x4(:,end) - [-1;0])];
end

logH  = log(H);
logE1 = log(E1);
logE2 = log(E2);
logE3 = log(E3);
logE4 = log(E4);

x  = linspace(min(log(H)), max(log(H)), 30);
y2 = (x - log(H(end))) + log(E1(end));
y4 = 2 * (x - log(H(end))) + log(E1(end));

figure(figNum);
hold off;

plot(logH,logE1,'*r', logH,logE1,'+b', logH,logE1,'og', logH(2:end),logE4(2:end),'sm');

hold on
plot(x,y2,'g',x,y4,'--r');
xlabel('Log(h)');
ylabel('Log(E)');
legend('Schema1','Schema2','Schema3','Schema4', strcat('droite de pente = ','1'), strcat('droite de pente =', '2'));
title('Graphique de comparaison des 4 schemas');

hold off;
print('Q4Comparaison.eps','-depsc');
res=1;
end
%===========================================================================
function [t1, x1, h1] = Schema1(f, t0, x0, h0)
tol     = 1e-3;
IterMax = 10;
theta   = 0.5;
[t1, x1, h1] = monThetaSChemaNewton(f, t0, x0, h0, tol, IterMax,theta);
end
%-----------------------------------------------------------------------------
function [t1, x1, h1] = Schema2(f, t0, x0, h0)
t1     = t0 + h0;
x1star = x0 + h0 * f(t0,x0);
x1     = x0 + (h0/2) * (f(t0,x0) + f(t1,x1star) );
h1     = h0 ;
end
%----------------------------------------------------------------------------
function [t1, x1, h1] = Schema3(f, t0, x0, h0)
t1       = t0 + h0;
x1star   = x0 + (h0/2) * (f(t0,x0) + f(t1,x0) );
x1       = x0 + (h0/2) * (f(t0,x0) + f(t1,x1star) );
h1       = h0 ;
end
%----------------------------------------------------------------------------
function [t1, x1, h1] = Schema4(f, t0, x0, h0)
t1       = t0 + h0;
x1star   = x0 + h0 * (f(t0 + h0/2.0, x0 + (h0/2.0) * f(t0,x0)) );
x1       = x0 + (h0/2) * (f(t0,x0) + f(t1,x1star) );
h1       = h0 ;
end
%----------------------------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%% UTILITAIRES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%% Solveur d'EDO %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [t,x,h]= solveurEDO(f,t0,tf,x0,N,schema)
t_temp=linspace(t0,tf,N+1);
h=t_temp(2) - t_temp(1);
t_temp =[];
x= [x0];
t= [t0];
for k=1:N
  [t1,x1,hh] = schema(f,t(k),x(:,k),h);
  x=[x,x1];
  t=[t,t1];
end
end
%%%%%%%%%%%%%%%%% Methode de Newton %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x,N] = Newton(g,x0,tol,iterMax)
x = x0;
iter = 0;
conv = false;
while(false == conv)
 [gn,Dgn] = g(x);
 delta    = Dgn\gn ;
 x = x - delta;
 iter = iter + 1;
 conv = (norm(delta)< tol * norm(x)) ||(iter > iterMax); 
end
N = iter;
end
%%%%%%%%%%%%%%%% Methode de Point Fixe %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x,N] = Picard(g,x0,tol,iterMax)
x = x0;
iter = 0;
conv = false;
while(false == conv)
 %calcul
 xx = g(x);
 delta = xx - x;
 x = xx;
 iter = iter + 1;
 conv = (norm(delta)< tol * norm(x)) ||(iter > iterMax); 
end
N = iter;
end
%%%%%%%%%%%%%%%%%%%% monSchemaEulerExplicite %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [t1, x1, h1] = monSchemaEulerExplicite(f, t0, x0, h0)
t1 = t0 + h0;
x1 = x0 + h0 * f(t0,x0);
h1 = h0 ;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [f,df] = monEDO(t,x)
x1 = x(1);
x2 = x(2);
u       = x1.^2 + x2.^2;
f  = [u * x2; - u * x1 ];
df(1,:) = [ 2 * x1 * x2, (u + 2*x2.^2)];
df(2,:) = [(-u - 2*x1.^2) , -2*x2*x1];
end
%%%%%%%%%%%%%%%%%%%%%%%% FIN %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%