% TP4ex4
% A exécuter avec des << a >> grands ....

clear;
format long;
a=input('La valeur du paramètre a= ');
n=input('La valeur de la taille de la subdivision n= ');
pas=1/n;
xm=[0:pas:(1-pas)];
xp=[pas:pas:1];
c=[pas/2:pas:1-pas/2];

% Méthode des rectangles. (1/n)
A=pas*sum(tp4ex4aux(a,xm));
B=pas*sum(tp4ex4aux(a,xp));

% Méthode du point central. (1/n)
C=pas*sum(tp4ex4aux(a,c));

% Méthode des trapèzes. (1/n^2)
T=(pas/2)*sum(tp4ex4aux(a,xm)+tp4ex4aux(a,xp));

% Méthode de Simpson. (1/n^4)

S=(pas/6)*sum(tp4ex4aux(a,xm)+tp4ex4aux(a,xp)+4*tp4ex4aux(a,c));

% Méthode  de Monte-Carlo.
y=rand(1,n);
MC=mean(tp4ex4aux(a,y));

% Méthode de Monte-Carlo avec rejet sur le rectangle [0;1]*[0;1/4]
z=1/4*rand(1,n);
MCC=mean(z<=tp4ex4aux(a,y))/4;


% Présentation des résultats

disp('Comparaison des évaluations')
disp('Méthode numériques classiques')
[A B C T S]
disp ('Monte-Carlo et Monte-Carlo avec rejet')
[MC MCC]

% Présentation graphique de la fonction tp4ex4aux(a,.)
plot(c,tp4ex4aux(a,c))
hold on
plot(y,z.*(tp4ex4aux(a,y)<z),'r.')
hold on
plot(y,z.*(tp4ex4aux(a,y)>=z),'g.')