tel/TPmatlab/integration/TP3h/points_poids_gauss_diago.m

   1 function [points,poids]=points_poids_gauss_diago(c,n)\r
   2 \r
   3 % points_poids_gauss_diago : calcul des poids et des points de Gauss par diagonalisation\r
   4 % \r
   5 % *********************************************************\r
   6 %\r
   7 % [points,poids]=points_poids_gauss_diago(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 %                calculés par diagonalisation d'une matrice tridiagonale.   \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 % \r
  30 % ************ Fonctions auxiliaires utilisées ************\r
  31 %\r
  32 %       aucune\r
  33 %\r
  34 % *********************************************************\r
  35 \r
  36 \r
  37 % Contrôles d'entrée\r
  38 % nombre d'arguments\r
  39 if nargin~=2\r
  40    error('nombre d''arguments de la fonction incorrect');\r
  41 end\r
  42 % autres tests éventuels \r
  43 if (fix(n)~=n) | (n<0)\r
  44    error('l''entier n doit être un entier naturel');\r
  45 end\r
  46 \r
  47 % Corps d'algorithme.\r
  48 points=zeros(n+1,n+1);\r
  49 poids=zeros(n+1,n+1);\r
  50 switch c\r
  51    case 1\r
  52       points(1,1)=0;\r
  53       poids(1,1)=2;\r
  54    case 2\r
  55       points(1,1)=0;\r
  56       poids(1,1)=pi;\r
  57    case 3\r
  58       points(1,1)=0;\r
  59       poids(1,1)=sqrt(pi);\r
  60    case 4\r
  61       points(1,1)=1;\r
  62       poids(1,1)=1;\r
  63 end\r
  64 if (n>=1)\r
  65    for k=1:n\r
  66       if (c<=3)\r
  67          B=zeros(1,k+1);\r
  68       else\r
  69          B=2*(0:k)+1;\r
  70       end\r
  71       switch c\r
  72          case 1\r
  73             gamma=1./(2*sqrt(1-1./(4*(1:k).^2)));\r
  74          case 2\r
  75             gamma=[1/sqrt(2) ones(1,k-1)/2];\r
  76          case 3\r
  77             gamma=sqrt((1:k)/2);\r
  78          case 4\r
  79             gamma=(1:k);\r
  80       end\r
  81       auxi=gallery('tridiag',gamma,B,gamma);\r
  82       if (k<=4) \r
  83          [Q,R]=eigs(auxi,k+1);\r
  84       else\r
  85          [Q,R]=eig(full(auxi));\r
  86       end\r
  87       R=diag(R);\r
  88       U=(sqrt(sum(Q.^2)));\r
  89       switch c\r
  90          case 1\r
  91             auxik=2;\r
  92          case 2\r
  93             auxik=pi;\r
  94          case 3\r
  95             auxik=sqrt(pi);\r
  96          case 4\r
  97             auxik=1;\r
  98       end\r
  99       Q=auxik*(Q(1,:)./U).^2;\r
 100       [auxipoints,ind]=sort((flipud(R))');\r
 101       auxi=fliplr(Q);\r
 102       auxipoids=auxi(ind);\r
 103       if (c<=3) \r
 104          p=fix(k/2);\r
 105          if ~(rem(k,2))\r
 106             auxi=(auxipoints(p+2:k+1)-auxipoints(p:-1:1))/2;      \r
 107             auxipoints(p+2:k+1)=auxi;\r
 108             auxipoints(1:p)=-fliplr(auxi);\r
 109             auxipoints(p+1)=0;\r
 110             auxi=(auxipoids(p+2:k+1)+auxipoids(p:-1:1))/2;      \r
 111             auxipoids(p+2:k+1)=auxi;\r
 112             auxipoids(1:p)=fliplr(auxi);\r
 113          else\r
 114             auxi=(auxipoints(p+2:k+1)-auxipoints(p+1:-1:1))/2;      \r
 115             auxipoints(p+2:k+1)=auxi;\r
 116             auxipoints(1:p+1)=-fliplr(auxi);\r
 117             auxi=(auxipoids(p+2:k+1)+auxipoids(p+1:-1:1))/2;      \r
 118             auxipoids(p+2:k+1)=auxi;\r
 119             auxipoids(1:p+1)=fliplr(auxi);\r
 120          end\r
 121       end\r
 122       points(k+1,1:k+1)=auxipoints;\r
 123       poids(k+1,1:k+1)=auxipoids;\r
 124    end\r
 125 end\r
 126 \r