Fonctions utiles

load(orthopoly)$

support_et_poids_de_gauss_legendre(n):=
/* renvoie le tableau [[x1,w1],...,[xn,wn]] composés des tableaux [points de support,poids]. 
adaptation de gauleg--Rybicki du numerical recipes */
block([polynome,derivee,racines,poids,x,precision:10^(-15)],
polynome:legendre_p(n,x),derivee:diff(polynome,x),
racines:realroots(polynome,precision),
map(lambda([u],ev(float([x,2/((1-x^2)*derivee^2)]),u)),racines)
)$
SPGL:support_et_poids_de_gauss_legendre;

Le format de sortie de support_et_poids_de_gauss_legendre :
[[w1,x1],[w2,x1],...]

SPGL(3);

Exemple d'utilisation

SP:SPGL(3);
f(x):=x^2;
block([P:length(SP),p:1,res:0],
      while (p <= P) do {res:res+SP[p][2]*f(SP[p][1]),p:p+1},
      res);

## Corrigé

P

P(x,sp):=block([P:length(sp),res:0,p:1],
 while (p<=P) do block([],
  res:res+sp[p][2]*exp(-x^2*(1+sp[p][1])^2/8),
  p:p+1),
 res:1/2+x/2*res/float(sqrt(2*%pi)))$

sp:SPGL(5)$

P(-1.96,sp);

P^(-1)

InvP(p,sp,eps):=block([x0:0,x1:2*eps],
 if (p < 0) then print(concat("p = ",p," est négatif"))
 else if (p >1) then print(concat("p = ",p," est > 1"))
 else if p=0 then x0:-inf
 else if p=1 then x0:inf
 else block([],
  while abs(x0-x1) > eps do block([],
   x0:x1,
   x1:float(x0-(P(x0,sp)-p)/(exp(-x0^2/2)/sqrt(2*%pi))))),
 x0)$

InvP(0.975,sp,1.0e-06);

InvP(0.999,sp,1.0e-06);