0
$\begingroup$

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

Improve this question
$\endgroup$
  • 1
    $\begingroup$ 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... $\endgroup$ – chris Aug 27 '20 at 20:18
  • $\begingroup$ The integral under consideration makes no sense (e. g. see encyclopediaofmath.org/wiki/Generalized_function ). That was noticed several times in this forum. $\endgroup$ – user64494 Aug 27 '20 at 23:29
  • 2
    $\begingroup$ @user64494 this was noted several times in this forum by you, specifically. For everybody else, integrals of this sort make perfect sense. $\endgroup$ – Roman Aug 28 '20 at 14:23
3
$\begingroup$

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

Improve this answer
$\endgroup$
  • $\begingroup$ 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 $\endgroup$ – NeAr Aug 27 '20 at 20:16
  • $\begingroup$ @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. $\endgroup$ – John Doty Aug 27 '20 at 22:04
1
$\begingroup$

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.

Improve this answer
$\endgroup$
  • $\begingroup$ This is the answer that I was lookin for. Thank you sir $\endgroup$ – NeAr Aug 28 '20 at 7:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.