Numerical integration with Dirac delta

Clash Royale CLAN TAG#URR8PPP

up vote
1
down vote

favorite

I have some complicated function depending on many arguments $x,y,z$ and parameter $a$ multiplied by Dirac delta of another function,

$$
tag 1 f(a,x,y,z) = g(a,x,y,z)delta(t(a,x,y,z))
$$

I want to perform numerical integration over all the variables. Then, independently on the parameter a which has been chosen, the result is 0. However, evaluating the same integral analytically, I obtain non-zero result.

What is a reason that Mathematica doesn't evaluate the numerical integral correctly, and how to force it to do this? It seems that I can't use limit representations of the Dirac delta because of the numerical integration.

Edit. My example is

f[a_,x_,y_,z_] =(a^4 + 6400 a^2 - 81920000) Exp[-0.01 x^2] Sqrt[y^2 - a^2]Sin[z] UnitStep[11 Pi/45 - z] UnitStep[z - 11 Pi/90] DiracDelta[a^2 + 2 Sqrt[y^2 - a^2] x Cos[z] - 2 y Sqrt[x^2 + 6400] + 6400]

I want to integrate it over the region $x in (0,3000), yin (a,3000), z in (0,pi)$. I can perform the first step - evaluate one of the integrals by rewriting the Dirac delta as, say,
$$
delta (t(a,x,y,z)) = fracdelta(x-x_1)t'(x_1)+fracdelta(x-x_2)t'(x_2),
$$
where $x_1,2 equiv x_1,2(a,y,z)$ are solutions of $t(a,x,y,z) = 0$, and then to introduce reduced function depending only on $a,y,z$,
$$
tag 2 f_1(a,y,z) = fracg(a,x_1,y,z)t'(x_1)+fracg(a,x_2,y,z)t'(x_2),
$$
where $g$ is defined in $(1)$. However, this step introduces extra work which I would like to avoid.

Simple numerical integration

NIntegrate[f[2, x, y, z], x, 0, 3000, y, 2, 3000, z, 0, Pi]

gives zero, but the integration performed with $(2)$ is non-zero.

numerical-integration special-functions

share|improve this question

edited 2 hours ago

asked 6 hours ago

John Taylor

557211

• 
  
  
  

  
  3
  

  
  
  
  
  
  

  
  The problem with DiracDelta is that it will evaluate to 0 when fed by a nonzero numerical argument. This way, $delta(t(a,x,y,z))$ will quite likely be interpreted as function that is zero almost everywhere (but in reality, $delta$ is not a function). It is better not to use it for numerical code. If I understand you correctly, you want to intergrate over the hypersurface defined by the equation $t(a,x,y,z) = 0$. Please, provide the concrete equations.
  – Henrik Schumacher
  5 hours ago
  
  

  
  
  
  


• 
  
  
  

  
  1
  

  
  
  
  
  
  

  
  In some applications it shoudl also be noted, that KroneckerDelta and DiscreteDelta and DiracDelta are not the same thing.
  – Johu
  3 hours ago
  

  
  
  
  


• 
  
  
  

  
  1
  

  
  
  
  
  
  

  
  Maybe something like NIntegrate[ g[a, x, y, z], a, x, y, z [Element] ImplicitRegion[t[a, x, y, z] == 0, a, x, y, z]]
  – Szabolcs
  3 hours ago
  
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  @HenrikSchumacher : I've added my example. Please take a look on it.
  – John Taylor
  2 hours ago
  
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  @Szabolcs : applied to my example (see the update of my question please) it doesn't give the result (only displays the code in the output).
  – John Taylor
  2 hours ago
  

  
  
  
  

add a comment | 

up vote
1
down vote

favorite

I have some complicated function depending on many arguments $x,y,z$ and parameter $a$ multiplied by Dirac delta of another function,

$$
tag 1 f(a,x,y,z) = g(a,x,y,z)delta(t(a,x,y,z))
$$

I want to perform numerical integration over all the variables. Then, independently on the parameter a which has been chosen, the result is 0. However, evaluating the same integral analytically, I obtain non-zero result.

What is a reason that Mathematica doesn't evaluate the numerical integral correctly, and how to force it to do this? It seems that I can't use limit representations of the Dirac delta because of the numerical integration.

Edit. My example is

f[a_,x_,y_,z_] =(a^4 + 6400 a^2 - 81920000) Exp[-0.01 x^2] Sqrt[y^2 - a^2]Sin[z] UnitStep[11 Pi/45 - z] UnitStep[z - 11 Pi/90] DiracDelta[a^2 + 2 Sqrt[y^2 - a^2] x Cos[z] - 2 y Sqrt[x^2 + 6400] + 6400]

I want to integrate it over the region $x in (0,3000), yin (a,3000), z in (0,pi)$. I can perform the first step - evaluate one of the integrals by rewriting the Dirac delta as, say,
$$
delta (t(a,x,y,z)) = fracdelta(x-x_1)t'(x_1)+fracdelta(x-x_2)t'(x_2),
$$
where $x_1,2 equiv x_1,2(a,y,z)$ are solutions of $t(a,x,y,z) = 0$, and then to introduce reduced function depending only on $a,y,z$,
$$
tag 2 f_1(a,y,z) = fracg(a,x_1,y,z)t'(x_1)+fracg(a,x_2,y,z)t'(x_2),
$$
where $g$ is defined in $(1)$. However, this step introduces extra work which I would like to avoid.

Simple numerical integration

NIntegrate[f[2, x, y, z], x, 0, 3000, y, 2, 3000, z, 0, Pi]

gives zero, but the integration performed with $(2)$ is non-zero.

numerical-integration special-functions

share|improve this question

edited 2 hours ago

asked 6 hours ago

John Taylor

557211

• 
  
  
  

  
  3
  

  
  
  
  
  
  

  
  The problem with DiracDelta is that it will evaluate to 0 when fed by a nonzero numerical argument. This way, $delta(t(a,x,y,z))$ will quite likely be interpreted as function that is zero almost everywhere (but in reality, $delta$ is not a function). It is better not to use it for numerical code. If I understand you correctly, you want to intergrate over the hypersurface defined by the equation $t(a,x,y,z) = 0$. Please, provide the concrete equations.
  – Henrik Schumacher
  5 hours ago
  
  

  
  
  
  


• 
  
  
  

  
  1
  

  
  
  
  
  
  

  
  In some applications it shoudl also be noted, that KroneckerDelta and DiscreteDelta and DiracDelta are not the same thing.
  – Johu
  3 hours ago
  

  
  
  
  


• 
  
  
  

  
  1
  

  
  
  
  
  
  

  
  Maybe something like NIntegrate[ g[a, x, y, z], a, x, y, z [Element] ImplicitRegion[t[a, x, y, z] == 0, a, x, y, z]]
  – Szabolcs
  3 hours ago
  
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  @HenrikSchumacher : I've added my example. Please take a look on it.
  – John Taylor
  2 hours ago
  
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  @Szabolcs : applied to my example (see the update of my question please) it doesn't give the result (only displays the code in the output).
  – John Taylor
  2 hours ago
  

  
  
  
  

add a comment | 

up vote
1
down vote

favorite

up vote
1
down vote

favorite

I have some complicated function depending on many arguments $x,y,z$ and parameter $a$ multiplied by Dirac delta of another function,

$$
tag 1 f(a,x,y,z) = g(a,x,y,z)delta(t(a,x,y,z))
$$

I want to perform numerical integration over all the variables. Then, independently on the parameter a which has been chosen, the result is 0. However, evaluating the same integral analytically, I obtain non-zero result.

What is a reason that Mathematica doesn't evaluate the numerical integral correctly, and how to force it to do this? It seems that I can't use limit representations of the Dirac delta because of the numerical integration.

Edit. My example is

f[a_,x_,y_,z_] =(a^4 + 6400 a^2 - 81920000) Exp[-0.01 x^2] Sqrt[y^2 - a^2]Sin[z] UnitStep[11 Pi/45 - z] UnitStep[z - 11 Pi/90] DiracDelta[a^2 + 2 Sqrt[y^2 - a^2] x Cos[z] - 2 y Sqrt[x^2 + 6400] + 6400]

I want to integrate it over the region $x in (0,3000), yin (a,3000), z in (0,pi)$. I can perform the first step - evaluate one of the integrals by rewriting the Dirac delta as, say,
$$
delta (t(a,x,y,z)) = fracdelta(x-x_1)t'(x_1)+fracdelta(x-x_2)t'(x_2),
$$
where $x_1,2 equiv x_1,2(a,y,z)$ are solutions of $t(a,x,y,z) = 0$, and then to introduce reduced function depending only on $a,y,z$,
$$
tag 2 f_1(a,y,z) = fracg(a,x_1,y,z)t'(x_1)+fracg(a,x_2,y,z)t'(x_2),
$$
where $g$ is defined in $(1)$. However, this step introduces extra work which I would like to avoid.

Simple numerical integration

NIntegrate[f[2, x, y, z], x, 0, 3000, y, 2, 3000, z, 0, Pi]

gives zero, but the integration performed with $(2)$ is non-zero.

numerical-integration special-functions

share|improve this question

edited 2 hours ago

asked 6 hours ago

John Taylor

557211

I have some complicated function depending on many arguments $x,y,z$ and parameter $a$ multiplied by Dirac delta of another function,

$$
tag 1 f(a,x,y,z) = g(a,x,y,z)delta(t(a,x,y,z))
$$

I want to perform numerical integration over all the variables. Then, independently on the parameter a which has been chosen, the result is 0. However, evaluating the same integral analytically, I obtain non-zero result.

What is a reason that Mathematica doesn't evaluate the numerical integral correctly, and how to force it to do this? It seems that I can't use limit representations of the Dirac delta because of the numerical integration.

Edit. My example is

f[a_,x_,y_,z_] =(a^4 + 6400 a^2 - 81920000) Exp[-0.01 x^2] Sqrt[y^2 - a^2]Sin[z] UnitStep[11 Pi/45 - z] UnitStep[z - 11 Pi/90] DiracDelta[a^2 + 2 Sqrt[y^2 - a^2] x Cos[z] - 2 y Sqrt[x^2 + 6400] + 6400]

I want to integrate it over the region $x in (0,3000), yin (a,3000), z in (0,pi)$. I can perform the first step - evaluate one of the integrals by rewriting the Dirac delta as, say,
$$
delta (t(a,x,y,z)) = fracdelta(x-x_1)t'(x_1)+fracdelta(x-x_2)t'(x_2),
$$
where $x_1,2 equiv x_1,2(a,y,z)$ are solutions of $t(a,x,y,z) = 0$, and then to introduce reduced function depending only on $a,y,z$,
$$
tag 2 f_1(a,y,z) = fracg(a,x_1,y,z)t'(x_1)+fracg(a,x_2,y,z)t'(x_2),
$$
where $g$ is defined in $(1)$. However, this step introduces extra work which I would like to avoid.

Simple numerical integration

NIntegrate[f[2, x, y, z], x, 0, 3000, y, 2, 3000, z, 0, Pi]

gives zero, but the integration performed with $(2)$ is non-zero.

numerical-integration special-functions

numerical-integration special-functions

share|improve this question

edited 2 hours ago

asked 6 hours ago

John Taylor

557211

share|improve this question

edited 2 hours ago

asked 6 hours ago

John Taylor

557211

share|improve this question

share|improve this question

edited 2 hours ago

edited 2 hours ago

edited 2 hours ago

asked 6 hours ago

John Taylor

557211

asked 6 hours ago

John Taylor

557211

asked 6 hours ago

John Taylor

557211

557211

• 
  
  
  

  
  3
  

  
  
  
  
  
  

  
  The problem with DiracDelta is that it will evaluate to 0 when fed by a nonzero numerical argument. This way, $delta(t(a,x,y,z))$ will quite likely be interpreted as function that is zero almost everywhere (but in reality, $delta$ is not a function). It is better not to use it for numerical code. If I understand you correctly, you want to intergrate over the hypersurface defined by the equation $t(a,x,y,z) = 0$. Please, provide the concrete equations.
  – Henrik Schumacher
  5 hours ago
  
  

  
  
  
  


• 
  
  
  

  
  1
  

  
  
  
  
  
  

  
  In some applications it shoudl also be noted, that KroneckerDelta and DiscreteDelta and DiracDelta are not the same thing.
  – Johu
  3 hours ago
  

  
  
  
  


• 
  
  
  

  
  1
  

  
  
  
  
  
  

  
  Maybe something like NIntegrate[ g[a, x, y, z], a, x, y, z [Element] ImplicitRegion[t[a, x, y, z] == 0, a, x, y, z]]
  – Szabolcs
  3 hours ago
  
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  @HenrikSchumacher : I've added my example. Please take a look on it.
  – John Taylor
  2 hours ago
  
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  @Szabolcs : applied to my example (see the update of my question please) it doesn't give the result (only displays the code in the output).
  – John Taylor
  2 hours ago
  

  
  
  
  

add a comment | 

• 
  
  
  

  
  3
  

  
  
  
  
  
  

  
  The problem with DiracDelta is that it will evaluate to 0 when fed by a nonzero numerical argument. This way, $delta(t(a,x,y,z))$ will quite likely be interpreted as function that is zero almost everywhere (but in reality, $delta$ is not a function). It is better not to use it for numerical code. If I understand you correctly, you want to intergrate over the hypersurface defined by the equation $t(a,x,y,z) = 0$. Please, provide the concrete equations.
  – Henrik Schumacher
  5 hours ago
  
  

  
  
  
  


• 
  
  
  

  
  1
  

  
  
  
  
  
  

  
  In some applications it shoudl also be noted, that KroneckerDelta and DiscreteDelta and DiracDelta are not the same thing.
  – Johu
  3 hours ago
  

  
  
  
  


• 
  
  
  

  
  1
  

  
  
  
  
  
  

  
  Maybe something like NIntegrate[ g[a, x, y, z], a, x, y, z [Element] ImplicitRegion[t[a, x, y, z] == 0, a, x, y, z]]
  – Szabolcs
  3 hours ago
  
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  @HenrikSchumacher : I've added my example. Please take a look on it.
  – John Taylor
  2 hours ago
  
  

  
  
  
  


• 
  
  
  

  
  

  
  
  
  
  
  

  
  @Szabolcs : applied to my example (see the update of my question please) it doesn't give the result (only displays the code in the output).
  – John Taylor
  2 hours ago
  

  
  
  
  

3

3

The problem with DiracDelta is that it will evaluate to 0 when fed by a nonzero numerical argument. This way, $delta(t(a,x,y,z))$ will quite likely be interpreted as function that is zero almost everywhere (but in reality, $delta$ is not a function). It is better not to use it for numerical code. If I understand you correctly, you want to intergrate over the hypersurface defined by the equation $t(a,x,y,z) = 0$. Please, provide the concrete equations.
– Henrik Schumacher
5 hours ago

The problem with DiracDelta is that it will evaluate to 0 when fed by a nonzero numerical argument. This way, $delta(t(a,x,y,z))$ will quite likely be interpreted as function that is zero almost everywhere (but in reality, $delta$ is not a function). It is better not to use it for numerical code. If I understand you correctly, you want to intergrate over the hypersurface defined by the equation $t(a,x,y,z) = 0$. Please, provide the concrete equations.
– Henrik Schumacher
5 hours ago

1

1

In some applications it shoudl also be noted, that KroneckerDelta and DiscreteDelta and DiracDelta are not the same thing.
– Johu
3 hours ago

In some applications it shoudl also be noted, that KroneckerDelta and DiscreteDelta and DiracDelta are not the same thing.
– Johu
3 hours ago

1

1

Maybe something like NIntegrate[ g[a, x, y, z], a, x, y, z [Element] ImplicitRegion[t[a, x, y, z] == 0, a, x, y, z]]
– Szabolcs
3 hours ago

Maybe something like NIntegrate[ g[a, x, y, z], a, x, y, z [Element] ImplicitRegion[t[a, x, y, z] == 0, a, x, y, z]]
– Szabolcs
3 hours ago

@HenrikSchumacher : I've added my example. Please take a look on it.
– John Taylor
2 hours ago

@HenrikSchumacher : I've added my example. Please take a look on it.
– John Taylor
2 hours ago

@Szabolcs : applied to my example (see the update of my question please) it doesn't give the result (only displays the code in the output).
– John Taylor
2 hours ago

@Szabolcs : applied to my example (see the update of my question please) it doesn't give the result (only displays the code in the output).
– John Taylor
2 hours ago

add a comment | 

2 Answers
2

active

oldest

votes

up vote
3
down vote

You can integrate symbolically of z and then use the result in NIntegrate.

newInt = Integrate[f[2, x, y, z], z, 0, [Pi]]

(* -(1/x)
 40947192 E^(-(x^2/
 100)) (HeavisideTheta[
 3202 - Sqrt[6400 + x^2] y + 
 x Sqrt[-4 + y^2] Cos[(11 [Pi])/90]] - 
 HeavisideTheta[
 3202 - Sqrt[6400 + x^2] y + x Sqrt[-4 + y^2] Cos[(11 [Pi])/45]]) *)

NIntegrate[newInt, x, 0, 3000, y, 2, 3000, 
 MinRecursion -> 4, Method -> "GaussKronrodRule", PrecisionGoal -> 4]

(* -4.85146*10^6 *)

NIntegrate[newInt, x, 0, 3000, y, 2, 3000, 
 MinRecursion -> 4, Method -> "GaussKronrodRule", PrecisionGoal -> 4, 
 WorkingPrecision -> 30]

(* -4.19942269494560346742215774341*10^7 *)

Take note of NIntegrate's messages.
Is this producing results you expect?

share|improve this answer

answered 1 hour ago

Anton Antonov

21.9k164108

add a comment | 

up vote
2
down vote

Here's an extension of Anton's approach:

newInt = Integrate[f[2, x, y, z], z, 0, π] /. 
 HeavisideTheta -> UnitStep; (* replacement by UnitStep is optional *)

The integrand has the form k * Exp[_] * (A-B) / x, where A and B are UnitStep (or HeavisideTheta) functions.

Block[A, B, (* there are two cases of UnitStep/HeavisideTheta *)
 A, B = Cases[newInt, UnitStep[u_] :> u, Infinity];
 reg = Reduce[ (* find reg where UnitStep[A]-UnitStep[B] is nonzero *)
 (A > 0 && B < 0 || A < 0 && B > 0) && 
 0 < x < 3000 && 2 < y < 3000, x, y]
 ];

Cases[DeleteCases[reg, _Equal && _], (* delete zero-measure subsets *)
 _[a_, ___, x, ___, b_] && _[c_, ___, y, ___, d_] :> 
 NIntegrate[newInt, x, a, b, y, c, d]]
Total[%]

NIntegrate::izero: Integral and error estimates are 0 on all integration subregions....

(*
 -4.23643*10^7, 0.
 -4.23643*10^7
*)

On the second part of reg, the exponential factor of newInt underflows.

share|improve this answer

edited 28 mins ago

answered 35 mins ago

Michael E2

141k11191457

add a comment | 

Your Answer

StackExchange.ifUsing("editor", function ()
return StackExchange.using("mathjaxEditing", function ()
StackExchange.MarkdownEditor.creationCallbacks.add(function (editor, postfix)
StackExchange.mathjaxEditing.prepareWmdForMathJax(editor, postfix, [["$", "$"], ["\\(","\\)"]]);
);
);
, "mathjax-editing");

StackExchange.ready(function()
var channelOptions =
tags: "".split(" "),
id: "387"
;
initTagRenderer("".split(" "), "".split(" "), channelOptions);

StackExchange.using("externalEditor", function()
// Have to fire editor after snippets, if snippets enabled
if (StackExchange.settings.snippets.snippetsEnabled)
StackExchange.using("snippets", function()
createEditor();
);

else
createEditor();

);

function createEditor()
StackExchange.prepareEditor(
heartbeatType: 'answer',
convertImagesToLinks: false,
noModals: false,
showLowRepImageUploadWarning: true,
reputationToPostImages: null,
bindNavPrevention: true,
postfix: "",
onDemand: true,
discardSelector: ".discard-answer"
,immediatelyShowMarkdownHelp:true
);



);

 

draft saved

draft discarded

Sign up or log in

StackExchange.ready(function ()
StackExchange.helpers.onClickDraftSave('#login-link');
);

Sign up using Google

Sign up using Facebook

Sign up using Email and Password

Post as a guest

Name Email

StackExchange.ready(
function ()
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmathematica.stackexchange.com%2fquestions%2f182855%2fnumerical-integration-with-dirac-delta%23new-answer', 'question_page');

);

Post as a guest

Name Email

2 Answers
2

active

oldest

votes

2 Answers
2

active

oldest

votes

active

oldest

votes

active

oldest

votes

up vote
3
down vote

You can integrate symbolically of z and then use the result in NIntegrate.

newInt = Integrate[f[2, x, y, z], z, 0, [Pi]]

(* -(1/x)
 40947192 E^(-(x^2/
 100)) (HeavisideTheta[
 3202 - Sqrt[6400 + x^2] y + 
 x Sqrt[-4 + y^2] Cos[(11 [Pi])/90]] - 
 HeavisideTheta[
 3202 - Sqrt[6400 + x^2] y + x Sqrt[-4 + y^2] Cos[(11 [Pi])/45]]) *)

NIntegrate[newInt, x, 0, 3000, y, 2, 3000, 
 MinRecursion -> 4, Method -> "GaussKronrodRule", PrecisionGoal -> 4]

(* -4.85146*10^6 *)

NIntegrate[newInt, x, 0, 3000, y, 2, 3000, 
 MinRecursion -> 4, Method -> "GaussKronrodRule", PrecisionGoal -> 4, 
 WorkingPrecision -> 30]

(* -4.19942269494560346742215774341*10^7 *)

Take note of NIntegrate's messages.
Is this producing results you expect?

share|improve this answer

answered 1 hour ago

Anton Antonov

21.9k164108

add a comment | 

up vote
3
down vote

You can integrate symbolically of z and then use the result in NIntegrate.

newInt = Integrate[f[2, x, y, z], z, 0, [Pi]]

(* -(1/x)
 40947192 E^(-(x^2/
 100)) (HeavisideTheta[
 3202 - Sqrt[6400 + x^2] y + 
 x Sqrt[-4 + y^2] Cos[(11 [Pi])/90]] - 
 HeavisideTheta[
 3202 - Sqrt[6400 + x^2] y + x Sqrt[-4 + y^2] Cos[(11 [Pi])/45]]) *)

NIntegrate[newInt, x, 0, 3000, y, 2, 3000, 
 MinRecursion -> 4, Method -> "GaussKronrodRule", PrecisionGoal -> 4]

(* -4.85146*10^6 *)

NIntegrate[newInt, x, 0, 3000, y, 2, 3000, 
 MinRecursion -> 4, Method -> "GaussKronrodRule", PrecisionGoal -> 4, 
 WorkingPrecision -> 30]

(* -4.19942269494560346742215774341*10^7 *)

Take note of NIntegrate's messages.
Is this producing results you expect?

share|improve this answer

answered 1 hour ago

Anton Antonov

21.9k164108

add a comment | 

up vote
3
down vote

up vote
3
down vote

You can integrate symbolically of z and then use the result in NIntegrate.

newInt = Integrate[f[2, x, y, z], z, 0, [Pi]]

(* -(1/x)
 40947192 E^(-(x^2/
 100)) (HeavisideTheta[
 3202 - Sqrt[6400 + x^2] y + 
 x Sqrt[-4 + y^2] Cos[(11 [Pi])/90]] - 
 HeavisideTheta[
 3202 - Sqrt[6400 + x^2] y + x Sqrt[-4 + y^2] Cos[(11 [Pi])/45]]) *)

NIntegrate[newInt, x, 0, 3000, y, 2, 3000, 
 MinRecursion -> 4, Method -> "GaussKronrodRule", PrecisionGoal -> 4]

(* -4.85146*10^6 *)

NIntegrate[newInt, x, 0, 3000, y, 2, 3000, 
 MinRecursion -> 4, Method -> "GaussKronrodRule", PrecisionGoal -> 4, 
 WorkingPrecision -> 30]

(* -4.19942269494560346742215774341*10^7 *)

Take note of NIntegrate's messages.
Is this producing results you expect?

share|improve this answer

answered 1 hour ago

Anton Antonov

21.9k164108

You can integrate symbolically of z and then use the result in NIntegrate.

newInt = Integrate[f[2, x, y, z], z, 0, [Pi]]

(* -(1/x)
 40947192 E^(-(x^2/
 100)) (HeavisideTheta[
 3202 - Sqrt[6400 + x^2] y + 
 x Sqrt[-4 + y^2] Cos[(11 [Pi])/90]] - 
 HeavisideTheta[
 3202 - Sqrt[6400 + x^2] y + x Sqrt[-4 + y^2] Cos[(11 [Pi])/45]]) *)

NIntegrate[newInt, x, 0, 3000, y, 2, 3000, 
 MinRecursion -> 4, Method -> "GaussKronrodRule", PrecisionGoal -> 4]

(* -4.85146*10^6 *)

NIntegrate[newInt, x, 0, 3000, y, 2, 3000, 
 MinRecursion -> 4, Method -> "GaussKronrodRule", PrecisionGoal -> 4, 
 WorkingPrecision -> 30]

(* -4.19942269494560346742215774341*10^7 *)

Take note of NIntegrate's messages.
Is this producing results you expect?

share|improve this answer

answered 1 hour ago

Anton Antonov

21.9k164108

share|improve this answer

share|improve this answer

answered 1 hour ago

Anton Antonov

21.9k164108

answered 1 hour ago

Anton Antonov

21.9k164108

answered 1 hour ago

Anton Antonov

21.9k164108

21.9k164108

add a comment | 

add a comment | 

up vote
2
down vote

Here's an extension of Anton's approach:

newInt = Integrate[f[2, x, y, z], z, 0, π] /. 
 HeavisideTheta -> UnitStep; (* replacement by UnitStep is optional *)

The integrand has the form k * Exp[_] * (A-B) / x, where A and B are UnitStep (or HeavisideTheta) functions.

Block[A, B, (* there are two cases of UnitStep/HeavisideTheta *)
 A, B = Cases[newInt, UnitStep[u_] :> u, Infinity];
 reg = Reduce[ (* find reg where UnitStep[A]-UnitStep[B] is nonzero *)
 (A > 0 && B < 0 || A < 0 && B > 0) && 
 0 < x < 3000 && 2 < y < 3000, x, y]
 ];

Cases[DeleteCases[reg, _Equal && _], (* delete zero-measure subsets *)
 _[a_, ___, x, ___, b_] && _[c_, ___, y, ___, d_] :> 
 NIntegrate[newInt, x, a, b, y, c, d]]
Total[%]

NIntegrate::izero: Integral and error estimates are 0 on all integration subregions....

(*
 -4.23643*10^7, 0.
 -4.23643*10^7
*)

On the second part of reg, the exponential factor of newInt underflows.

share|improve this answer

edited 28 mins ago

answered 35 mins ago

Michael E2

141k11191457

add a comment | 

up vote
2
down vote

Here's an extension of Anton's approach:

newInt = Integrate[f[2, x, y, z], z, 0, π] /. 
 HeavisideTheta -> UnitStep; (* replacement by UnitStep is optional *)

The integrand has the form k * Exp[_] * (A-B) / x, where A and B are UnitStep (or HeavisideTheta) functions.

Block[A, B, (* there are two cases of UnitStep/HeavisideTheta *)
 A, B = Cases[newInt, UnitStep[u_] :> u, Infinity];
 reg = Reduce[ (* find reg where UnitStep[A]-UnitStep[B] is nonzero *)
 (A > 0 && B < 0 || A < 0 && B > 0) && 
 0 < x < 3000 && 2 < y < 3000, x, y]
 ];

Cases[DeleteCases[reg, _Equal && _], (* delete zero-measure subsets *)
 _[a_, ___, x, ___, b_] && _[c_, ___, y, ___, d_] :> 
 NIntegrate[newInt, x, a, b, y, c, d]]
Total[%]

NIntegrate::izero: Integral and error estimates are 0 on all integration subregions....

(*
 -4.23643*10^7, 0.
 -4.23643*10^7
*)

On the second part of reg, the exponential factor of newInt underflows.

share|improve this answer

edited 28 mins ago

answered 35 mins ago

Michael E2

141k11191457

add a comment | 

up vote
2
down vote

up vote
2
down vote

Here's an extension of Anton's approach:

newInt = Integrate[f[2, x, y, z], z, 0, π] /. 
 HeavisideTheta -> UnitStep; (* replacement by UnitStep is optional *)

The integrand has the form k * Exp[_] * (A-B) / x, where A and B are UnitStep (or HeavisideTheta) functions.

Block[A, B, (* there are two cases of UnitStep/HeavisideTheta *)
 A, B = Cases[newInt, UnitStep[u_] :> u, Infinity];
 reg = Reduce[ (* find reg where UnitStep[A]-UnitStep[B] is nonzero *)
 (A > 0 && B < 0 || A < 0 && B > 0) && 
 0 < x < 3000 && 2 < y < 3000, x, y]
 ];

Cases[DeleteCases[reg, _Equal && _], (* delete zero-measure subsets *)
 _[a_, ___, x, ___, b_] && _[c_, ___, y, ___, d_] :> 
 NIntegrate[newInt, x, a, b, y, c, d]]
Total[%]

NIntegrate::izero: Integral and error estimates are 0 on all integration subregions....

(*
 -4.23643*10^7, 0.
 -4.23643*10^7
*)

On the second part of reg, the exponential factor of newInt underflows.

share|improve this answer

edited 28 mins ago

answered 35 mins ago

Michael E2

141k11191457

Here's an extension of Anton's approach:

newInt = Integrate[f[2, x, y, z], z, 0, π] /. 
 HeavisideTheta -> UnitStep; (* replacement by UnitStep is optional *)

The integrand has the form k * Exp[_] * (A-B) / x, where A and B are UnitStep (or HeavisideTheta) functions.

Block[A, B, (* there are two cases of UnitStep/HeavisideTheta *)
 A, B = Cases[newInt, UnitStep[u_] :> u, Infinity];
 reg = Reduce[ (* find reg where UnitStep[A]-UnitStep[B] is nonzero *)
 (A > 0 && B < 0 || A < 0 && B > 0) && 
 0 < x < 3000 && 2 < y < 3000, x, y]
 ];

Cases[DeleteCases[reg, _Equal && _], (* delete zero-measure subsets *)
 _[a_, ___, x, ___, b_] && _[c_, ___, y, ___, d_] :> 
 NIntegrate[newInt, x, a, b, y, c, d]]
Total[%]

NIntegrate::izero: Integral and error estimates are 0 on all integration subregions....

(*
 -4.23643*10^7, 0.
 -4.23643*10^7
*)

On the second part of reg, the exponential factor of newInt underflows.

share|improve this answer

edited 28 mins ago

answered 35 mins ago

Michael E2

141k11191457

share|improve this answer

share|improve this answer

edited 28 mins ago

edited 28 mins ago

edited 28 mins ago

answered 35 mins ago

Michael E2

141k11191457

answered 35 mins ago

Michael E2

141k11191457

answered 35 mins ago

Michael E2

141k11191457

141k11191457

add a comment | 

add a comment | 

 

draft saved

draft discarded

 

draft saved

draft discarded

Sign up or log in

StackExchange.ready(function ()
StackExchange.helpers.onClickDraftSave('#login-link');
);

Sign up using Google

Sign up using Facebook

Sign up using Email and Password

Post as a guest

Name Email

StackExchange.ready(
function ()
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmathematica.stackexchange.com%2fquestions%2f182855%2fnumerical-integration-with-dirac-delta%23new-answer', 'question_page');

);

Post as a guest

Name Email

Sign up or log in

StackExchange.ready(function ()
StackExchange.helpers.onClickDraftSave('#login-link');
);

Sign up using Google

Sign up using Facebook

Sign up using Email and Password

Post as a guest

Name Email

Sign up or log in

StackExchange.ready(function ()
StackExchange.helpers.onClickDraftSave('#login-link');
);

Sign up using Google

Sign up using Facebook

Sign up using Email and Password

Post as a guest

Name Email

Sign up or log in

StackExchange.ready(function ()
StackExchange.helpers.onClickDraftSave('#login-link');
);

Sign up using Google

Sign up using Facebook

Sign up using Email and Password

Sign up using Google

Sign up using Facebook

Sign up using Email and Password

Post as a guest

Name Email

Name Email

Name

Email

Name Email

Name

Email