tel/TPmatlab/integration/TP3i/verifie_ordre_gauss.m

   1 function eta=verifie_ordre_gauss(c,n,points,poids)\r
   2 \r
   3 \r
   4 %       verifie_ordre_gauss : vérification de l'ordre des formules Gaussiennes\r
   5 %\r
   6 % *********************************************************\r
   7 %\r
   8 %       eta=verifie_ordre_gauss(c,n,points,poids) calcule numériquement le max\r
   9 %       de |I_{n,k}-int_a ^b x^k w(x)dx| pour 0<=k<=2n+1, où I_{n,k}\r
  10 %       est l'intégrale approchée de x^kw(x) sur l'intervalle conventionnel (a,b) \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 : Gauss-Legendre,\r
  16 %            c=2 : Gauss-Tchebytchev,\r
  17 %            c=3 : Gauss-Hermite,\r
  18 %            c=4 : Gauss-Laguerre.\r
  19 %     * points : tableau de taille (q+1)*(q+1) (avec n<=q) tel que \r
  20 %                pour tout l dans {0,...,q}, points(l+1,1:l+1) contient\r
  21 %                les l+1 points x_i de la formule (à l+1 points)\r
  22 %     * poids  : tableau de taille (n+1)*(n+1) tel que \r
  23 %                pour tout l dans {0,...,n}, poids(l+1,1:l+1) contient\r
  24 %                les l+1 poids D_i de la formule (à l+1 points)\r
  25 %       variable de sortie\r
  26 %   *   eta= max_{0<=k<=2n+1} |I_{n,k}-int_a ^b x^k w(x)dx|\r
  27 %\r
  28 %   eta doit être nul.\r
  29 %\r
  30 % \r
  31 % ************ Fonctions auxiliaires utilisées ************\r
  32 %\r
  33 %  expression_part_symb\r
  34 %\r
  35 % *********************************************************\r
  36 \r
  37 % Contrôles d'entrée\r
  38 % nombre d'arguments\r
  39 if nargin~=4\r
  40    error('nombre d''arguments de la fonction incorrect');\r
  41 end\r
  42 % autres tests éventuels \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 p=length(poids);\r
  47 if (2*n+1>=p)\r
  48    error('n est trop grand');\r
  49 end\r
  50 \r
  51 \r
  52 % Corps d'algorithme.\r
  53 er=0;\r
  54 for k=0:2*n+1\r
  55    I=sum(poids(k+1,1:k+1).*((points(k+1,1:k+1)).^k));\r
  56    erloc=expression_part_symb(c,k)-I;\r
  57    if ~(isnumeric(erloc))\r
  58       erloc=double(vpa(erloc));\r
  59    end\r
  60    erloc=abs(erloc);\r
  61    er=max(er,erloc);\r
  62 end\r
  63 eta=er;