tel/TPmatlab/integration/TP3b/points_poids_gauss.m

   1 function [points,poids]=points_poids_gauss(c,n)\r
   2 \r
   3 \r
   4 % points_poids_gauss : calcul des poids et des points de Gauss\r
   5 % \r
   6 % *********************************************************\r
   7 %\r
   8 % [points,poids]=points_poids_gauss(c,n) renvoie les points D_i et les points x_i\r
   9 %                 des formules d'intégration à l+1 points pour 0<=l<=n, pour \r
  10 %         l'une des quatre méthodes : Legendre, Tchebytchev, Hermite ou Laguerre.\r
  11 %\r
  12 %       variables d'entrées : \r
  13 %    *  n : entier naturel.\r
  14 %     * c définit la méthode : \r
  15 %            c=1 : méthode de Gauss-Legendre,\r
  16 %            c=2 : méthode de Gauss-Tchebychev,\r
  17 %            c=3 : méthode de Gauss-Hermite,\r
  18 %            c=4 : méthode de Gauss-Laguerre;\r
  19 %\r
  20 %       variables de sortie\r
  21 %     * points : tableau de taille (n+1)*(n+1) tel que \r
  22 %                pour tout l dans {0,...,n}, points(l+1,1:l+1) contient\r
  23 %                les l+1 points x_i de la formule (à l+1 points)\r
  24 %     * poids  : tableau de taille (n+1)*(n+1) tel que \r
  25 %                pour tout l dans {0,...,n}, poids(l+1,1:l+1) contient\r
  26 %                les l+1 poids D_i de la formule (à l+1 points)\r
  27 %\r
  28 %\r
  29 % ************ Fonctions auxiliaires utilisées ************\r
  30 %\r
  31 %  calcul_poly_gauss\r
  32 %\r
  33 % *********************************************************\r
  34 \r
  35 \r
  36 % Contrôles d'entrée\r
  37 % nombre d'arguments\r
  38 if nargin~=2\r
  39    error('nombre d''arguments de la fonction incorrect');\r
  40 end\r
  41 % autres tests éventuels \r
  42 if (fix(n)~=n) | (n<0)\r
  43    error('l''entier n doit être un entier naturel');\r
  44 end\r
  45 if ~((c==1) | (c==2) | (c==3) | (c==4))\r
  46    error('c doit être un entier égal à 1, 2, 3 ou 4');\r
  47 end\r
  48 \r
  49 % Corps d'algorithme.\r
  50 poids=zeros(n+1,n+1);\r
  51 points=zeros(n+1,n+1);\r
  52 if (c~=2) \r
  53    P=calcul_poly_gauss(c,n+1);\r
  54    for k=0:n\r
  55       C=(sort(roots(P(k+2,1:k+2))))';\r
  56       [points(k+1,1:k+1),ind]=sort(C);   \r
  57       switch c\r
  58          case 1\r
  59             auxipoids=2./((k+1)*polyval(polyder(P(k+2,1:k+2)),C).*polyval(P(k+1,1:k+1),C));\r
  60          case 3\r
  61             auxipoids=(2^(k+2)*sqrt(pi)*factorial(k+1))./((polyval(polyder(P(k+2,1:k+2)),C)).^2);\r
  62          case 4\r
  63             auxipoids=1./(C.*(polyval(polyder(P(k+2,1:k+2)),C)).^2);\r
  64       end\r
  65       poids(k+1,1:k+1)=auxipoids(ind);      \r
  66    end\r
  67 else\r
  68    for k=0:n\r
  69       points(k+1,1:k+1)=cos(((2*(0:k)+1)/(k+1))*pi/2);   \r
  70       poids(k+1,1:k+1)=(pi/(k+1))*ones(1,k+1);\r
  71    end\r
  72 end\r