tel/TPmatlab/integration/TP3i/demo_trace_compar_gauss.m

   1 % le script trace_compar_gauss permet de comparer les différents\r
   2 % calculs de points et de poids de gauss en comparant l'intégrale exacte \r
   3 % et l'intégrale approchée d'une fonction selon que les poids et points\r
   4 % ont été calculés : \r
   5 % *  par méthode directe;\r
   6 % *  par méthode symbolique;\r
   7 % * ou par méthode par diagonalisation,\r
   8 % et ce pour chacune des quatre méthodes définies par c.\r
   9 % les poids et les points seront préalablements calculés selon chacune des trois méthodes\r
  10 % dans les tableaux : \r
  11 %          pointsleg poidsleg ...\r
  12 %          pointslegsymb poidslegsymb ...\r
  13 %          pointslegdiag poidslegdiag ...\r
  14 %          pointstch poidstch ...\r
  15 %          pointstchsymb poidstchsymb ...\r
  16 %          pointstchdiag poidstchdiag ...\r
  17 %          pointsher poidsher ...\r
  18 %          pointshersymb poidssherymb ...\r
  19 %          pointsherdiag poidsherdiag ...\r
  20 %          pointslag poidslag ...\r
  21 %          pointslagsymb poidslagsymb ...\r
  22 %          pointslagdiag poidslagdiag ...\r
  23 % voir script demo_stocke_points_poids \r
  24 %\r
  25 % \r
  26 % ************ Fonctions auxiliaires utilisées ************\r
  27 %\r
  28 % saisiefonction, int_gauss\r
  29 %\r
  30 % *********************************************************\r
  31 \r
  32 \r
  33 \r
  34 \r
  35 % logarithme de l'erreur nulle prise égale à log10zero\r
  36 log10zero=-35;\r
  37 \r
  38 \r
  39 clear xA yA xB yB xC yC Iexa erA erB erC;\r
  40 disp('choix de la méthode : ')\r
  41 c=input('entrez 1, 2, 3 ou 4 : ');\r
  42 disp(' ');\r
  43 disp('Entrée de la fonction r telle que r=f/w sur (a,b)  où w est le poids ') ;\r
  44 disp(' ');\r
  45 fcn=saisiefonction;\r
  46 Iexa=input('entrez la valeur exacte de l''intégrale de r=f/w sur (a,b) : ') ;\r
  47 erA=zeros(1,nAmax+1);\r
  48 yA=zeros(1,nAmax+1);\r
  49 erB=zeros(1,nBmax+1);\r
  50 yB=zeros(1,nBmax+1);\r
  51 erC=zeros(1,nCmax+1);\r
  52 yC=zeros(1,nCmax+1);\r
  53 switch c\r
  54    case 1\r
  55       pointsA=pointsleg;\r
  56       poidsA=poidsleg;\r
  57       pointsB=pointslegsymb;\r
  58       poidsB=poidslegsymb;\r
  59       pointsC=pointslegdiag;\r
  60       poidsC=poidslegdiag;\r
  61    case 2\r
  62       pointsA=pointstch;\r
  63       poidsA=poidstch;\r
  64       pointsB=pointstchsymb;\r
  65       poidsB=poidstchsymb;\r
  66       pointsC=pointstchdiag;\r
  67       poidsC=poidstchdiag;\r
  68    case 3\r
  69       pointsA=pointsher;\r
  70       poidsA=poidsher;\r
  71       pointsB=pointshersymb;\r
  72       poidsB=poidshersymb;\r
  73       pointsC=pointsherdiag;\r
  74       poidsC=poidsherdiag;\r
  75    case 4\r
  76       pointsA=pointslag;\r
  77       poidsA=poidslag;\r
  78       pointsB=pointslagsymb;\r
  79       poidsB=poidslagsymb;\r
  80       pointsC=pointslagdiag;\r
  81       poidsC=poidslagdiag;\r
  82 end\r
  83 \r
  84 \r
  85 for k=0:nAmax\r
  86    auxi=int_gauss(k,pointsA,poidsA,fcn);\r
  87    erA(k+1)=abs(auxi-Iexa);\r
  88 end\r
  89 for k=0:nBmax\r
  90    auxi=int_gauss(k,pointsB,poidsB,fcn);\r
  91    auxi=vpa(auxi)-Iexa;\r
  92    if ~(isnumeric(auxi))\r
  93       auxi=double(auxi);\r
  94    end\r
  95    erB(k+1)=abs(auxi);\r
  96 end\r
  97 for k=0:nCmax\r
  98    auxi=int_gauss(k,pointsC,poidsC,fcn);\r
  99    erC(k+1)=abs(auxi-Iexa);\r
 100 end\r
 101 indnul=find(erA==0);\r
 102 indnonnul=find(erA~=0);\r
 103 yA(indnul)=log10zero;\r
 104 yA(indnonnul)=log10(erA(indnonnul));\r
 105 indnul=find(erB==0);\r
 106 indnonnul=find(erB~=0);\r
 107 yB(indnul)=log10zero;\r
 108 yB(indnonnul)=log10(erB(indnonnul));\r
 109 indnul=find(erC==0);\r
 110 indnonnul=find(erC~=0);\r
 111 yC(indnul)=log10zero;\r
 112 yC(indnonnul)=log10(erC(indnonnul));\r
 113 xA=log10(1:nAmax+1);\r
 114 xB=log10(1:nBmax+1);\r
 115 xC=log10(1:nCmax+1);\r
 116 clf;\r
 117 hold on;\r
 118 plot(xA,yA,'b',xB,yB,'r',xC,yC,'g');\r
 119 xlabel('logarithme 10 du nombre de points d''intégration');\r
 120 ylabel('logarithme 10 de l''erreur d''intégration');\r
 121 legend('directe','symbolique','diagonalisation',0);\r
 122 hold off;\r
 123 \r
 124 \r
 125    \r
 126 \r
 127 \r
 128 \r
 129 \r
 130 \r
 131 \r
 132 \r
 133 \r
 134 \r