TP MP* n°4 (SCILAB) Intégration numérique

// on affiche suffisamment de chiffres significatifs
format(20) 
funcprot(0)
 
// taux de variation pour différents h
function L = diff_avant(f,x,d)
h = 1
for n = 1:d do
  h = h/2
  L(n,1) = h  
  L(n,2) = (f(x+h)-f(x))/h
end
endfunction
 
// x -> sqrt(1+x)
function y = sqa(x)
  y = sqrt(1+x)
endfunction
 
// on observe
//diff_avant(sqa,0,30)
 
 
 
 
// différence avant d'ordre 3  pour différents h
function L = diff_avant3(f,x,d)
L = []
h=1
for n = 1:d do
  h = h/2
  L(n,1) = h 
  L(n,2) = (-f(x+2*h)/6 +f(x+h)-0.5*f(x) -f(x-h)/3)/h
end
endfunction
 
// on observe
//diff_avant3(sqa,0,30)
 
 
/////////////////////////////////////////////////////////////
//   NEWTON COTES
/////////////////////////////////////////////////////////////
 
 
deff('[y] = p(x)', 'y = %pi  - 4/(1+x^2)')
 
 
// rectangles
function R = Rectangles(f,a,b,N)
  h = (b-a)/N
  pivots = a + h*[1:N]
  R = h*sum(feval(pivots,f))
endfunction
 
 
 
 
// trapèzes
function T = Trapezes(f,a,b,N)
  h = (b-a)/N
  if N == 1 then pivots = []
  else pivots = a + h*[1:N-1]
  end
  T = h*(sum(feval(pivots,f))+0.5*f(a)+0.5*f(b))
endfunction
 
 
 
// Simpson
function S = Simpson(f,a,b,N)
  h = (b-a)/N
  pivots_ext = a + h*[1:N-1]
  pivots_mil = a + 0.5*h*(2*[1:N]-1)
  S = (h/6)*(f(a)+f(b)+2*sum(feval(pivots_ext,f))+4*sum(feval(pivots_mil,f)))
endfunction
 
T=[];  
for k = 1:10 do
T = T; [Rectangles(p,0,1,2^k) Trapezes(p,0,1,2^k) Simpson(p,0,1,2^k)]
end;
T
 
 
//////////////////////////////////////////////////////////////
//             ROMBERG
///////////////////////////////////////////////////////////
 
 
function R = Romberg(f,a,b,m)
  R = zeros(m,m)
  N = 1
  for k = 1:m do
    //Ap = h*(sum(feval(a+h*[1:2:2^k],f)))
    A(k,1)=Trapezes(f,a,b,N)//[0.5*A(k-1)+Ap]
    N = 2*N
  end
  for t=2:m do
    for k=t:m do
      A(k,t) = (4^{t-1}*A(k,t-1)-A(k-1,t-1))/(4^{t-1}-1)
    end
  end
  R = A 
endfunction
  
  
  Romberg(p,0,1,5)
  
  deff('[y] = a(x)', 'y = 2/(1+4*x^2)')
  
  deff('[y] = e(x)', 'y = exp(x)+1-exp(1)')
  deff('[y] = c(x)', 'y = cos(x^2)')
  
  
//   Romberg(a,-1,2,6)
//   Romberg(e,0,1,6)
//   Romberg(c,0,1.25,6)
  
//////////////////////////////////////////////////////////////
//          MONTE CARLO
///////////////////////////////////////////////////////////