Re: question about DiracDelta

*To*: mathgroup at smc.vnet.net*Subject*: [mg69624] Re: [mg69610] question about DiracDelta*From*: Andrzej Kozlowski <akoz at mimuw.edu.pl>*Date*: Sun, 17 Sep 2006 22:46:05 -0400 (EDT)

On 17 Sep 2006, at 19:58, dimmechan at yahoo.com wrote: > Dear everyone, > > Mathematica is able to integrate numerically functions with peaks. > E.g. > > > Clear["Global*"] > > Plot[Exp[-x^2], {x, -10, 10}, PlotRange -> All] > > Integrate[Exp[-x^2], {x, -1000, 1000}] > N[%] > NIntegrate[Exp[-x^2], {x, -1000, 1000}] > Sqrt[Pi]*Erf[1000] > 1.7724538509055159 > 1.7724538509054621 > > With the default settings for all options, NIntegrate misses the peak > at x=0 here. > > NIntegrate[Exp[-x^2], {x, -10000, 10000}] > NIntegrate::ploss: Numerical integration stopping due to loss of > precision. \ > Achieved neither the requested PrecisionGoal nor AccuracyGoal; suspect > one of \ > the following: highly oscillatory integrand or the true value of the > integral \ > is 0. If your integrand is oscillatory on a (semi-)infinite interval > try \ > using the option Method->Oscillatory in NIntegrate. > 0. > > However, > > NIntegrate[Exp[-x^2], {x, -10000, 10000}, MinRecursion -> 6, > MaxRecursion -> 12] > NIntegrate[Exp[-x^2], {x, -10000, -10, 0, 10, 10000}] > 1.7724538509054268 > 1.7724538509055054 > > My question is if, somehow, can NIntegrate treat functions with > "spikes" like DiracDelta? > > Integrate[DiracDelta[x - 1], {x, -1, 3}] > Block[{Message},NIntegrate[DiracDelta[x - 1], {x, -1,1, 3}]] > Block[{Message},NIntegrate[DiracDelta[x - 1], {x, -1, 1, 3}, > WorkingPrecision -> 50, MinRecursion -> 10, MaxRecursion -> 20]] > 1 > 0. > 0. > > I know that I am asking essentially Mathematica to do a numerical > integral of a function that will be zero everywhere except x=1. So, > unless NIntegrate causes the function to be sampled at x=1, it will > evaluate to zero. > > Anyway, here are the sampled points used by NIntegrate. > > Block[{Message}, ListPlot[Reap[a1=NIntegrate[DiracDelta[x - 1], {x, > -1, > 1, 3}, EvaluationMonitor :> Sow[x]]][[2,1]]]] > Block[{Message}, ListPlot[a2=Reap[NIntegrate[DiracDelta[x - 1], {x, > -1, > 1, 3}, WorkingPrecision->50,MinRecursion -> 10, MaxRecursion -> 20, > EvaluationMonitor :> Sow[x]]][[2,1]]]] > > Length[Select[a1, 0.95 < #1 < 1.05 & ]] > 2 > Length[Select[a2, 0.95 < #1 < 1.05 & ]] > 1688 > > Thanks in advance for any help. > > Dimitris Anagnostou > Numerical functions are not suitable for dealing with distributions for rather obvious reasons. So you have to use Integrate and not NIntegrate when you have an expression involving DiracDelta. Integrate[DiracDelta[x - 1], {x, -1, 3}] 1 However, you will get much better results if you download from MathSource Maxim Rytin's package Piecewise and use his function PieciwiseIntegrate. It pretty much leaves Integrate in the dust in this kind of problems. Andrzej Kozlowski