# DiracDelta in numerical integration

I need to integrate a complicated expression numerically. However, I am not sure how to implement the DiracDelta and its corresponding derivative into code that makes sense.

Here is an example

func =  2 (3-2 x) x^2 HeavisideTheta[1-x]-4 (11/2 DiracDelta[-1+x]-x^2 (6+5 x) HeavisideTheta[1-x])-2 (1/6 DiracDelta[-1+x]-5/3 x^3 HeavisideTheta[1-x]+1/6 (DiracDelta^\[Prime])[1-x])

NIntegrate[func,{x,0,1}]/.HeavisideTheta[0]->1//Simplify


The result of the integral is -15/2 with symbolic integration. I know that DiracDelta is tricky in numerical integration since it will evaluated as zero, however, I really need a trick how to deal with this numerically, especially the derivate

• I am not sure it will suit your purpose but you could replace the Dirac by a Gaussian of finite small thickness and make that thickness smaller to check that your integral has converged? Make sure to ask the integration to sample the position of the Dirac though... – chris Aug 27 '20 at 20:18
• The integral under consideration makes no sense (e. g. see encyclopediaofmath.org/wiki/Generalized_function ). That was noticed several times in this forum. – user64494 Aug 27 '20 at 23:29
• @user64494 this was noted several times in this forum by you, specifically. For everybody else, integrals of this sort make perfect sense. – Roman Aug 28 '20 at 14:23

You can use it as you tried, I only corrected the derivative and modified the integration range:

Integrate[2 (3 - 2 x) x^2 HeavisideTheta[1 - x] -4 (11/2 DiracDelta[-1 + x] - x^2 (6 + 5 x) HeavisideTheta[1 - x]) -2 (1/6 DiracDelta[-1 + x] - 5/3 x^3 HeavisideTheta[1 - x] +
1/6 DiracDelta'[1 - x]), {x, 0, 2}]
(*-(15/2)*)


With the integration range {x,0,1} the integration evaluates to

(*89/6 - (67 HeavisideTheta[0])/3*)


For numerically applications you could use limit definition 1/Sqrt[2 Pi \[CurlyEpsilon]] Exp[-(x^2/(2 \[CurlyEpsilon]))] , \[CurlyEpsilon] sufficiently small

• I know how to do this symbolically. That's not the issue. I want to compute it numerically to understand how to deal with the diracdelta function – NeAr Aug 27 '20 at 20:16
• @NeAr The delta function cannot be numerically integrated: it is zero everywhere except at one point where it isn't a number. In numerical work, you must use a function that is nonzero over some suitable interval. – John Doty Aug 27 '20 at 22:04

In general this works. Use a function like:

d[x_] = a E^(-100000 x^2)


and we know the integral is 1.

Integrate[d[x], {x, -Infinity, Infinity}] == 1

Solve[%, a] // Flatten

a = a /. %


also

dd[x_] = d'[x]


Change your func using d.

func = 2 (3 - 2 x) x^2 HeavisideTheta[1 - x] -
4 (11/2 d[-1 + x] - x^2 (6 + 5 x) HeavisideTheta[1 - x]) -
2 (1/6 d[-1 + x] - 5/3 x^3 HeavisideTheta[1 - x] + 1/6 dd'[1 - x])


Unfortunately, because the numerical function d has some spread, integrating it when the spike is at an end point is not too accurate:

NIntegrate[func, {x, 0, 1}, WorkingPrecision -> 30]
(*3.66666666666666666666666666622*)


We can expand the integral a little, however.

NIntegrate[func, {x, 0, 1.1}, WorkingPrecision -> 30]
(*-7.49914231815272796510999142299*)


or expanding further:

Integrate[func, {x, 0, 2}] // Simplify // N[#, 20] &
(*-7.5000000000000000000*)


It is fairly accurate when the integral nearly covers the entire spike.

• This is the answer that I was lookin for. Thank you sir – NeAr Aug 28 '20 at 7:20