O Mathematica se integra muito bem usando o “código” que escrevi

Nov 14 2020

Estou tentando fazer o Mathematica fazer o seguinte: Divida o intervalo [0,1] em $n$intervalos iguais. Em seguida, em cada intervalo, aplique a quadratura gaussiana para 2 pontos.

Tentei usar o método da pergunta mky aqui: Usando regras de integração Newton-Cotes compostas no Mathematica

Mas os resultados que estou obtendo são muito precisos.

Aqui está o código:

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

Então o que isso deveria fazer é apenas fazer a quadratura gaussiana com 2 pontos em $[0,1]$ mas estou obtendo a resposta com pelo menos $10^{-5}$precisão que não deveria estar acontecendo. O que eu fiz errado?

Respostas

9 MichaelE2 Nov 14 2020 at 09:07

Você tem sorte (no que diz respeito a obter poucos erros):

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

O código acima reproduz o NIntegrateresultado:

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

Por que temos sorte neste caso? O erro é igual ao integral da diferença da função e do polinômio de interpolação por meio de abscissas abs, que tem aproximadamente a mesma área acima e abaixo do xeixo:

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