tel/TPmatlab/integration/TP3g/points_poids_gauss_symb.m

   1 function [points,poids]=points_poids_gauss_symb(c,n)\r
   2 \r
   3 % points_poids_gauss_symb : calcul des poids et des points de Gauss (symbolique)\r
   4 % \r
   5 % *********************************************************\r
   6 %\r
   7 % [points,poids]=points_poids_gauss_symb(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 %       en symbolique \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 : utilisation des formules du cours (dérivée) \r
  32 %\r
  33 %\r
  34 % ************ Fonctions auxiliaires utilisées ************\r
  35 %\r
  36 %  calcul_poly_gauss_symb, polyval_symb, polyder_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    if (c==3)\r
  72       auxis=sym(4)*sqrt(spi);\r
  73    end\r
  74    syms x;\r
  75    for k=1:n\r
  76       sk=sym(k);\r
  77       C=polyval_symb(L(k+2,1:k+2),x);\r
  78       E=solve(C);\r
  79       points(k+1,1:k+1)=E;   \r
  80       switch c\r
  81          case 1\r
  82             Pp=polyder_symb(L(k+2,1:k+2));\r
  83             for j=1:k+1\r
  84                poids(k+1,j)=2./((sk+1)*polyval_symb(Pp,E(j)).*polyval_symb(L(k+1,1:k+1),E(j)));\r
  85             end\r
  86          case 3\r
  87             Pp=polyder_symb(L(k+2,1:k+2));\r
  88             auxis=2*auxis*(sk+1);\r
  89             for j=1:k+1\r
  90                poids(k+1,j)=auxis/(polyval_symb(Pp,E(j)))^2;\r
  91             end\r
  92          case 4\r
  93             Pp=polyder_symb(L(k+2,1:k+2));\r
  94             for j=1:k+1\r
  95                poids(k+1,j)=1/(E(j)*(polyval_symb(Pp,E(j)))^2);\r
  96             end\r
  97       end\r
  98    end\r
  99 else\r
 100    for k=0:n\r
 101       sk=sym(k);\r
 102       points(1+k,1:1+k)=cos(((2*sym(0:k))+1)/(sk+1)*spi/2);\r
 103       poids(1+k,1:1+k)=spi/(sk+1)*ones(1,k+1);\r
 104    end\r
 105 end\r
 106 \r