tel/TPmatlab/integration/TP3b/calcul_poly_gauss.m

   1 function P=calcul_poly_gauss(c,n)\r
   2 \r
   3 %       calcul_poly_gauss : calcul des polynômes othogonaux de Gauss\r
   4 %\r
   5 % *********************************************************\r
   6 %\r
   7 % P=calcul_poly_gauss(c,n) renvoie les polynômes de Gauss P_0,P_1,...,P_n\r
   8 %   pour les quatre méthodes génériques (Legendre, Tchebytchev, Hermite \r
   9 %   et Laguerre.\r
  10 %\r
  11 %       variables d'entrées : \r
  12 %    *  n : entier naturel.\r
  13 %     * c définit la méthode : \r
  14 %            c=1 : Gauss-Legendre,\r
  15 %            c=2 : Gauss-Tchebytchev,\r
  16 %            c=3 : Gauss-Hermite,\r
  17 %            c=4 : Gauss-Laguerre.\r
  18 %       variables de sortie\r
  19 %     * P : tableau de taille (n+1)*(n+1) tel que \r
  20 %           pour tout i dans {1,...,n+1}, P(i,1:i) contient\r
  21 %           le polynôme P_{i-1} (coefficient dominant en tête)\r
  22 % \r
  23 % \r
  24 % ************ Fonctions auxiliaires utilisées ************\r
  25 %\r
  26 %       aucune\r
  27 %\r
  28 % *********************************************************\r
  29 %\r
  30 \r
  31 \r
  32 \r
  33 \r
  34 % Contrôles d'entrée\r
  35 % nombre d'arguments\r
  36 if nargin~=2\r
  37    error('nombre d''arguments de la fonction incorrect');\r
  38 end\r
  39 % autres tests éventuels \r
  40 if (fix(n)~=n) | (n<0)\r
  41    error('l''entier n doit être un entier naturel');\r
  42 end\r
  43 if ~((c==1) | (c==2) | (c==3) | (c==4))\r
  44    error('c doit être un entier égal à 1, 2, 3 ou 4');\r
  45 end\r
  46 \r
  47 \r
  48 % Corps d'algorithme.\r
  49 P=zeros(n+1,n+1);\r
  50 if (c~=2)\r
  51    P(1,1)=1;\r
  52    Lb=[];\r
  53    La=[1];\r
  54    for k=0:n-1\r
  55       switch(c)\r
  56          case 1\r
  57             auxi0=1/(k+1);\r
  58             auxi1=(2*k+1)*auxi0;\r
  59             auxi2=-k*auxi0;\r
  60             Lc=[auxi1*La 0]+[0 0 auxi2*Lb];\r
  61          case 3\r
  62             Lc=[2*La 0]-[0 0 2*k*Lb];\r
  63          case 4\r
  64             auxi0=1/(k+1);\r
  65             auxi1=(2*k+1)*auxi0;\r
  66             auxi2=-k*auxi0;\r
  67             Lc=[-auxi0*La 0]+[0 auxi1*La]+[0 0 auxi2*Lb];             \r
  68       end\r
  69       Lb=La;\r
  70       La=Lc;\r
  71       P(k+2,1:k+2)=Lc;    \r
  72    end\r
  73 else\r
  74    P(1,1)=1;\r
  75    if (n>=1)\r
  76       P(2,1:2)=[1 0];\r
  77       Lb=1;\r
  78       La=[1 0];\r
  79       for k=1:n-1      \r
  80          Lc=[2*La 0]-[0 0 Lb];\r
  81          Lb=La;\r
  82          La=Lc;\r
  83          P(k+2,1:k+2)=Lc;    \r
  84       end\r
  85    end\r
  86 end\r