//calcul de l'integrale d'une fonction f 
// 
//choix de la fonction 
fonction=x_dialog('fonction a integrer ','y=x.*sin(x)'); 
a=evstr(x_dialog('borne min ','0'));
b=evstr(x_dialog('borne max ','1'));
pmax=evstr(x_dialog('nombre de sous intervalles (2 puissance p)','13'));
//
deff('y=f(x)',fonction);
//
// trapezes et Romberg
//
write(%io(2),'  n      I(trapezes)      I(Romberg): ') 
val3a=0;val3=(f(a)+f(b))/2*(b-a);A3=[];A4=[];R=[val3];
for p=1:pmax
n=2^p;
h=(b-a)/n; 
for i=1:2:n-1
  val3a=val3a+f(a+i*h); 
end 
val3=(val3a+(f(a)+f(b))/2)*h;
R=[R;val3];
//
R1=R;
for i=1:p
  for j=p:-1:i
    R1(j+1)=(4^i*R1(j+1)-R1(j))/(4^i-1);
  end
end
//
disp([n,val3,R1(p+1)]);
end 
//
// Monte Carlo (probabilite '1.96sigma'=95%) 
n=1;
xb=[a+(b-a)*rand()];
xint=f(xb);
write(%io(2),' n     I(Monte Carlo)     intervalle de confiance à 95% ') 
for p=1:pmax
n0=n;
n=n+(2^p-2^(p-1)); 
xb=[xb;a+(b-a)*rand(n-n0+1,1)];
val4=(b-a)/n*sum(f(xb)); 
xvar4=(b-a)/(n-1)*sum((f(xb)-val4).^2);
xerr4=1.96*xvar4/sqrt(n);
disp([n,val4,xerr4]);
end
//
//valeur approchee Scilab
val1=intg(a,b,f); 
disp('valeur calculee par Scilab:');disp(val1);
