]> AND Private Git Repository - cours-mesi.git/blob - tel/TPmatlab/integration/TP3i/verifie_ordre_gauss.m
Logo AND Algorithmique Numérique Distribuée

Private GIT Repository
j
[cours-mesi.git] / 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;