Mathematica s'intègre trop bien en utilisant le «code» que j'ai écrit

Nov 14 2020

J'essaye de faire en sorte que Mathematica fasse ce qui suit: Divisez l'intervalle [0,1] en $n$intervalles égaux. Ensuite, sur chaque intervalle, appliquez la quadrature gaussienne pour 2 points.

J'ai essayé d'utiliser la méthode de mky question ici: Utilisation des règles d'intégration composites Newton-Cotes dans Mathematica

Mais les résultats que j'obtiens sont bien trop précis.

Voici le code:

n = 1;                                                                                                            
NIntegrate[Sin[x]/x, Evaluate@Flatten@{x, Subdivide[0., 1., n]}, 
 Method -> {"GaussBerntsenEspelidRule", "Points" -> 2}, 
 MaxRecursion -> 0]

Donc, ce que cela devrait faire, c'est simplement faire la quadrature gaussienne avec 2 points sur $[0,1]$ mais j'obtiens la réponse avec au moins $10^{-5}$précision qui ne devrait pas se produire. Qu'ai-je fait de mal?

Réponses

9 MichaelE2 Nov 14 2020 at 09:07

Vous avez juste de la chance (en ce qui concerne une erreur faible):

{abs, wts, err} = 
 NIntegrate`GaussBerntsenEspelidRuleData[2, MachinePrecision]
(*
  {{0.0469101, 0.230765, 0.5, 0.769235, 0.95309},
   {0.118463, 0.239314, 0.284444, 0.239314, 0.118463},
   {0.155257, -0.439701, 0.568889, -0.439701, 0.155257}}
*)
(Sin[x]/x /. x -> abs).wts
(Sin[x]/x /. x -> abs).err
(*
  0.946083       <-- integral estimate
  0.0000639286   <-- estimated error bound
*)
(Sin[x]/x /. x -> abs).wts - Integrate[Sin[x]/x, {x, 0, 1}]
(*
  3.31957*10^-14  <-- actual error (less than the bound)
*)

Le code ci-dessus reproduit le NIntegraterésultat:

(Sin[x]/x /. x -> abs).wts -
 NIntegrate[Sin[x]/x, Evaluate@Flatten@{x, Subdivide[0., 1., nn]}, 
  Method -> {"GaussBerntsenEspelidRule", "Points" -> 2}, 
  MaxRecursion -> 0]
(*
  0.
*)

Pourquoi avons-nous de la chance dans ce cas? L'erreur est égale à l'intégrale de la différence de la fonction et du polynôme d'interpolation par les abscisses abs, qui a à peu près la même aire au-dessus et au-dessous de l' xaxe:

Plot[
 InterpolatingPolynomial[Transpose@{abs, (Sin[x]/x /. x -> abs)}, x] -
  Sin[x]/x // Evaluate,
 {x, 0, 1}]