Mathematica si integra troppo bene usando il "codice" che ho scritto

Nov 14 2020

Sto cercando di fare in modo che Mathematica esegua quanto segue: Dividi l'intervallo [0,1] in $n$intervalli uguali. Quindi su ogni intervallo applica la quadratura gaussiana per 2 punti.

Ho provato a utilizzare il metodo dalla domanda di mky qui: Utilizzo delle regole di integrazione Composite Newton-Cotes in Mathematica

Ma i risultati che sto ottenendo sono troppo accurati.

Ecco il codice:

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

Quindi quello che dovrebbe fare è semplicemente fare la quadratura gaussiana con 2 punti $[0,1]$ ma almeno sto ottenendo la risposta $10^{-5}$accuratezza che non dovrebbe accadere. Cos'ho fatto di sbagliato?

Risposte

9 MichaelE2 Nov 14 2020 at 09:07

Sei solo fortunato (per quanto riguarda l'errore basso):

{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)
*)

Il codice sopra riproduce il NIntegraterisultato:

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

Perché siamo fortunati in questo caso? L'errore è pari all'integrale della differenza della funzione e del polinomio interpolante per ascisse abs, che ha all'incirca la stessa area sopra e sotto l' xasse:

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