//
//  programme 3
//
write(%io(2),'calcul de l''integrale de f(x)=sinx -x') 
write(%io(2),'entre 0 et 1 par la methode des trapezes ') 
write(%io(2),'et par la methode de Monte Carlo')   
write(%io(2),'//')   
//
nmax=evstr(x_dialog('nombre de points maximum ','2000'));
// 
//choix de la fonction 
//
a=0; 
b=1; 
deff('y=f(x)','y=sin(x)-x') 
//
//valeur approchee calculee par Scilab
//
val1=integrate('sin(x)-x','x',a,b); 
//
//valeur exacte 
//
val2=-cos(b)-(b^2)/2+cos(a)+(a^2)/2; 
//
//trapezes sur 2^p points
//
p=1;n=1;
val3a=0;A3=[];A4=[];
//
write(%io(2),'  trapezes sur 2^p points: ') 
while (2*n<nmax) 
  p=p+1; 
  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;
  
  write (%io(2),[n,abs(val3-val2)],'(f6.0,1x,e14.6)')  
  A3=[A3;n];A4=[A4;abs(val3-val2)];
end 
//
//
// Monte Carlo (probabilite '1.96sigma'=95%) 
//
n=0;
xint=0; xint2=0;xerr4=1E5;
A1=[];A2a=[];A2b=[];
write(%io(2),'  Monte Carlo: ') 
while (n<nmax)
  n0=n;
  n=ceil(nmax/10)+n; 
  for i=n0:n
    xb=a+(b-a)*rand();
    xint=xint+f(xb); 
    xint2=xint2+f(xb)^2;
  end 
  val4=xint/n; 
  xvar4=max(0,sqrt(xint2/n-val4^2));
  xerr4=1.96*xvar4/sqrt(n);
  write (%io(2),[n,xerr4,val4-val2],'(f6.0,1x,2(e14.6))')  
  A1=[A1;n];A2a=[A2a;xerr4];A2b=[A2b;abs(val4-val2)];
end
//
// trace des resultats
//
xbasc()
plot2d1('gll',[A1,A1],[A2a,A2b],[1,2],'151','erreur MC a 96%@erreur MC calculee',[2,min([A2a;A2b;A4]),max([A1;A3]),max([A2a;A2b;A4])]);
plot2d1('gll',A3,A4,3,'100','erreur trapezes');
xtitle('calcul approche de l integrale de sin(x)-x entre 0 et 1','n','erreur')
 
