Mixing
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
asked 6 hours ago
John Taylor
557211
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
asked 6 hours ago
John Taylor
557211
3
The problem withDiracDeltais that it will evaluate to0when 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, thatKroneckerDeltaandDiscreteDeltaandDiracDeltaare not the same thing.
– Johu
3 hours ago
1
Maybe something likeNIntegrate[ 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
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
asked 6 hours ago
John Taylor
557211
asked 6 hours ago
John Taylor
557211
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 withDiracDeltais that it will evaluate to0when 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, thatKroneckerDeltaandDiscreteDeltaandDiracDeltaare not the same thing.
– Johu
3 hours ago
1
Maybe something likeNIntegrate[ 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 withDiracDeltais that it will evaluate to0when 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, thatKroneckerDeltaandDiscreteDeltaandDiracDeltaare not the same thing.
– Johu
3 hours ago
1
Maybe something likeNIntegrate[ 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?
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.
answered 35 mins ago
Michael E2
141k11191457
add a comment |Â
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?
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?
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?
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?
answered 1 hour ago
Anton Antonov
21.9k164108
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.
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.
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.
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.
answered 35 mins ago
Michael E2
141k11191457
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 |Â
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
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
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
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
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
3
The problem with
DiracDeltais that it will evaluate to0when 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
KroneckerDeltaandDiscreteDeltaandDiracDeltaare 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