tel/TPmatlab/integration/TP3g/points_poids_gauss_vand.m

   1 function [points,poids]=points_poids_gauss_vand(c,n)\r
   2 \r
   3 % points_poids_gauss_vand : calcul des poids et des points de Gauss (avec Vandermonde)\r
   4 % \r
   5 % *********************************************************\r
   6 %\r
   7 % [points,poids]=points_poids_gauss_vand(c,n) renvoie les points D_i et les points x_i\r
   8 %                    des formules d'intégration à l+1 points pour 0<=l<=n, pour \r
   9 %                l'une des quatre méthodes : Legendre, Tchebytchev, Hermite ou Laguerre,\r
  10 %                éventuellement en formel  \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 %   Pour ces deux tableau, évaluation numérique (augmenter le nombre\r
  29 %   de chiffre avec la fonction digits) ou formelle si possible.\r
  30 %   méthode :  points : recherche de zéro,\r
  31 %              poids : calcul à partir des points en inversant la matrice de Vandermonde\r
  32 %\r
  33 %\r
  34 % ************ Fonctions connexes utilisées ***************\r
  35 %\r
  36 %  calcul_poly_gauss_symb, polyval_symb, expression_part_symb\r
  37 %\r
  38 % *********************************************************\r
  39 \r
  40 % Contrôles d'entrée\r
  41 % nombre d'arguments\r
  42 if nargin~=2\r
  43    error('nombre d''arguments de la fonction incorrect');\r
  44 end\r
  45 % autres tests éventuels \r
  46 if (fix(n)~=n) | (n<0)\r
  47    error('l''entier n doit être un entier naturel');\r
  48 end\r
  49 \r
  50 % Corps d'algorithme.\r
  51 poids=zeros(n+1,n+1);\r
  52 points=zeros(n+1,n+1);\r
  53 poids=sym(poids);\r
  54 points=sym(points);\r
  55 spi=sym(pi);\r
  56 if (c~=2) \r
  57    L=calcul_poly_gauss_symb(c,n+1);\r
  58    if ((1<=c) & (c<=3))\r
  59       points(1,1:1)=0;\r
  60    else\r
  61       points(1,1)=1;\r
  62    end\r
  63    switch c\r
  64      case 1\r
  65          poids(1,1)=2;\r
  66       case 3\r
  67          poids(1,1)=sqrt(spi);\r
  68       case 4\r
  69          poids(1,1)=1;\r
  70    end\r
  71    U=zeros(n+1,1);\r
  72    U=sym(U);\r
  73    syms x;\r
  74    for k=1:n\r
  75       C=polyval_symb(L(k+2,1:k+2),x);\r
  76       E=solve(C);\r
  77       points(k+1,1:k+1)=E;   \r
  78       for i=0:k\r
  79          U(i+1)=expression_part_symb(c,i);\r
  80       end\r
  81       M(1,1:k+1)=sym(ones(1,k+1));\r
  82       for i=2:1+k\r
  83          for j=1:k+1\r
  84             M(i,j)=(points(k+1,j))^(i-1);\r
  85          end\r
  86       end\r
  87       poids(k+1,1:k+1)=M\U(1:k+1);\r
  88    end\r
  89 else\r
  90    for k=0:n\r
  91       sk=sym(k);\r
  92       points(1+k,1:1+k)=cos((2*sym(0:k)+1)/(sk+1)*spi/2);\r
  93       poids(1+k,1:1+k)=spi/(sk+1)*ones(1,k+1);\r
  94    end\r
  95 end\r