## assumptions – Why is this integrand not integrating to a Bessel function?

I know from the identities of Bessel functions that the following is true:
$$J_{m}left( x right) = frac{ 1 }{ 2 pi i^{m} } int_{0}^{2 pi} dphi e^{i left( x cos{phi} – m phi right)} tag{0}$$

However, when I try to integrate the following, Mathematica doesn’t actually integrate the result. It just spits back the input.

``````Assuming(Element(k, Reals) && Element(m, Integers), Integrate(Exp(I (m x + k Cos(x))),{x, 0, 2 Pi}))
``````

Should I allow `m` to be more general, i.e., any numeric value (not necessarily fond of fractional Bessel function indices) or am I missing some other assumption that is preventing Mathematica from giving the proper result?

(I doubt this matters but for completeness, I have v12.1.1.0 and am using a Macbook Pro laptop.)

## Problem statement

I have some experimental data of measurement values vs. temperature.
Now, I want to find the spectrum which causes the measured temperature dependence.

The measurement values are the integrals of

## The forward calculation

I will now illustrate the problem with an example. There is an yet unknown sensitivity spectrum.

``````spectrum(lam_) := Exp(-((lam - 5)/2)^4);
Plot(spectrum(lam), {lam, 0, 10}, AxesLabel -> {"wavelength", "sensitivity"})
``````

And there is a known spectral intensity

``````intensity(lam_, T_) := 1/(T^4 lam^5) 1/(Exp(1/(lam *T)) - 1)
Plot(Table(intensity(lam, T), {T, 0.03, 0.1, 0.01}) // Evaluate, {lam, 0, 10},
AxesLabel -> {"wavelength", "intensity"}) // Quiet
``````

The measured signals are the areas under these curves

``````Plot(Table(spectrum(lam)*intensity(lam, T), {T, 0.03, 0.1, 0.01}) // Evaluate, {lam, 0, 10},
AxesLabel -> {"wavelength", "intensity"}) // Quiet

data = Table({T,NIntegrate(spectrum(lam)*intensity(lam, T), {lam, 1, 10})}, {T, 0.03, 0.1, 0.01});
ListPlot(data, AxesLabel -> {"temperature", "signal"})
``````

## What’s the solution?

I thought of nonlinear curve fitting for the control points of an interpolation function as was proposed
in this question by Hugh and Alexey Popkov.
Using the sample data, the result should look like the `spectrum(lam)` from above.

But how to make this code work?

``````nOfControlPoints = 5;
controlPoints = Subdivide(1, 10, nOfControlPoints - 1)

model(y : {__Real}) := Interpolation(Transpose({controlPoints, y}))
fctn = NIntegrate(model(Array(y, nOfControlPoints))(lam)*intensity(lam, T), {lam, 1, 10});

nlm = NonlinearModelFit(data, fctn, Array(y, nOfControlPoints), T)
``````

## numerical integration – Why does the Kernel CRASH without any Message when I use “NIntegrate” for the following complicated compiled Integrand?

When I NIntegrate the Integrand with the particular parameter showing below the Kernel crash without any Message.

But if I change the “range” parameter a little bit, the Kernel won’t crash. The crash problem seems to occur only under certain parameters.

Because I want to make a density plot of the integral as a function of the parameters, this crash problem under certain parameters always breaks the density plot process.

PS: This crash problem happens both on Mac and Win10

``````(*Integrand*)

expression1 = {{-0.03783737745791316` +
52.07594879537986` (k1^2 + k2^2) -
0.0006566020385220937` (GothicCapitalB)3,
-0.00020215783672119025` ((GothicCapitalB)1 -
I (GothicCapitalB)2), (0.` + 0.` I) +
1/2 (273.35324964239567` (k1 - I k2)^3 -
22.668570452264646` (k1 + I k2)^3),
0.` + I (k1 - I k2) (-1.0301133331497063` -
370.16816003256605` (k1^2 +
k2^2))}, {-0.00020215783672119025` ((GothicCapitalB)1 +
I (GothicCapitalB)2), -0.03783737745791316` +
52.07594879537986` (k1^2 + k2^2) +
0.0006566020385220937` (GothicCapitalB)3,
0.` - I (k1 + I k2) (-1.0301133331497063` -
370.16816003256605` (k1^2 + k2^2)), (0.` + 0.` I) +
1/2 (22.668570452264646` (k1 - I k2)^3 -
273.35324964239567` (k1 + I k2)^3)}, {(0.` + 0.` I) +
1/2 (-22.668570452264646` (k1 - I k2)^3 +
273.35324964239567` (k1 + I k2)^3),
0.` + I (k1 - I k2) (-1.0301133331497063` -
370.16816003256605` (k1^2 + k2^2)), -0.03304414880767263` +
95.54303239292248` (k1^2 + k2^2) -
0.0015031216532114644` (GothicCapitalB)3,
-0.00017807130082869512` ((GothicCapitalB)1 -
I (GothicCapitalB)2)}, {0.` -
I (k1 + I k2) (-1.0301133331497063` -
370.16816003256605` (k1^2 + k2^2)), (0.` + 0.` I) +
1/2 (-273.35324964239567` (k1 - I k2)^3 +
22.668570452264646` (k1 +
I k2)^3), -0.00017807130082869512` ((GothicCapitalB)1 +
I (GothicCapitalB)2), -0.03304414880767263` +
95.54303239292248` (k1^2 + k2^2) +
0.0015031216532114644` (GothicCapitalB)3}};
expression2 = {{104.15189759075972` k1s, 0,
1/2 (820.059748927187` (k1s - I k2)^2 -
68.00571135679394` (k1s + I k2)^2), (0.` -
740.3363200651321` I) k1s (k1s - I k2) +
I (-1.0301133331497063` -
370.16816003256605` (k1s^2 + k2^2))}, {0,
104.15189759075972` k1s, (0.` + 740.3363200651321` I) k1s (k1s +
I k2) - I (-1.0301133331497063` -
370.16816003256605` (k1s^2 + k2^2)),
1/2 (68.00571135679394` (k1s - I k2)^2 -
820.059748927187` (k1s + I k2)^2)}, {1/
2 (-68.00571135679394` (k1s - I k2)^2 +
820.059748927187` (k1s + I k2)^2), (0.` -
740.3363200651321` I) k1s (k1s - I k2) +
I (-1.0301133331497063` - 370.16816003256605` (k1s^2 + k2^2)),
191.08606478584497` k1s,
0}, {(0.` + 740.3363200651321` I) k1s (k1s + I k2) -
I (-1.0301133331497063` - 370.16816003256605` (k1s^2 + k2^2)),
1/2 (-820.059748927187` (k1s - I k2)^2 +
68.00571135679394` (k1s + I k2)^2), 0,
191.08606478584497` k1s}};
expression3 = {{104.15189759075972` k2s, 0,
1/2 ((0.` - 820.059748927187` I) (k1 - I k2s)^2 - (0.` +
68.00571135679394` I) (k1 +
I k2s)^2), -1.0301133331497063` - (0.` +
740.3363200651321` I) (k1 - I k2s) k2s -
370.16816003256605` (k1^2 + k2s^2)}, {0,
104.15189759075972` k2s, -1.0301133331497063` + (0.` +
740.3363200651321` I) (k1 + I k2s) k2s -
370.16816003256605` (k1^2 + k2s^2),
1/2 ((0.` - 68.00571135679394` I) (k1 - I k2s)^2 - (0.` +
820.059748927187` I) (k1 + I k2s)^2)}, {1/
2 ((0.` + 68.00571135679394` I) (k1 - I k2s)^2 + (0.` +
820.059748927187` I) (k1 +
I k2s)^2), -1.0301133331497063` - (0.` +
740.3363200651321` I) (k1 - I k2s) k2s -
370.16816003256605` (k1^2 + k2s^2), 191.08606478584497` k2s,
0}, {-1.0301133331497063` + (0.` + 740.3363200651321` I) (k1 +
I k2s) k2s - 370.16816003256605` (k1^2 + k2s^2),
1/2 ((0.` + 820.059748927187` I) (k1 - I k2s)^2 + (0.` +
68.00571135679394` I) (k1 + I k2s)^2), 0,
191.08606478584497` k2s}};
HamiltCompiled =
Compile({{k1, _Real}, {k2, _Real}, {(GothicCapitalB)1, _Real}, {
(GothicCapitalB)2, _Real}, {(GothicCapitalB)3, _Real}},
Evaluate(expression1),
CompilationOptions -> {"ExpressionOptimization" -> True});
Dk1HamiltCompiled =
Compile({{k1s, _Real}, {k2, _Real}, {(GothicCapitalB)1, _Real}, {
(GothicCapitalB)2, _Real}, {(GothicCapitalB)3, _Real}},
Evaluate(expression2),
CompilationOptions -> {"ExpressionOptimization" -> True});
Dk2HamiltCompiled =
Compile({{k1, _Real}, {k2s, _Real}, {(GothicCapitalB)1, _Real}, {
(GothicCapitalB)2, _Real}, {(GothicCapitalB)3, _Real}},
Evaluate(expression3),
CompilationOptions -> {"ExpressionOptimization" -> True});
OccupiedCurvatureCompiled((Mu)_?NumericQ, k1_?NumericQ,
k2_?NumericQ, (GothicCapitalB)1_, (GothicCapitalB)2_,
(GothicCapitalB)3_) :=
Block({Mx1, Mx2, values, vectors, curvaturesList},
{values, vectors} =
Eigensystem@
HamiltCompiled(k1,
k2, (GothicCapitalB)1, (GothicCapitalB)2,
(GothicCapitalB)3);
{values, vectors} = {values((#)), vectors((#))} &(
Ordering@values);
Mx1 = vectors(Conjugate).(Dk1HamiltCompiled(k1,
k2, (GothicCapitalB)1, (GothicCapitalB)2,
(GothicCapitalB)3)).vectors(Transpose);
Mx2 = vectors(Conjugate).(Dk2HamiltCompiled(k1,
k2, (GothicCapitalB)1, (GothicCapitalB)2,
(GothicCapitalB)3)).vectors(Transpose);
curvaturesList =
Re(I Total((# - #(Transpose)))) &(
Mx1*Mx2(Transpose)*
Table(If(n1 == n2, 0, 1/(values((n1)) - values((n2)))^2), {n1,
1, 4}, {n2, 1, 4}));
);
``````
``````(*No Crash*)
With({range =
0.0006` +
0.00025*19.`, (Mu) = -0.034637500000000015`, (GothicCapitalB)1 =
0., (GothicCapitalB)2 = 19.`, (GothicCapitalB)3 = 0.},
NIntegrate(
OccupiedCurvatureCompiled((Mu), k1,
k2, (GothicCapitalB)1, (GothicCapitalB)2, (GothicCapitalB)3),
{k1, -range, range}, {k2, -range, range}, MinRecursion -> 2,
PrecisionGoal -> 3, AccuracyGoal -> 3))
``````
``````(*Crash*)
With({range =
0.0005` +
0.00025*19.`, (Mu) = -0.034637500000000015`, (GothicCapitalB)1 =
0., (GothicCapitalB)2 = 19.`, (GothicCapitalB)3 = 0.},
NIntegrate(
OccupiedCurvatureCompiled((Mu), k1,
k2, (GothicCapitalB)1, (GothicCapitalB)2, (GothicCapitalB)3),
{k1, -range, range}, {k2, -range, range}, MinRecursion -> 2,
PrecisionGoal -> 3, AccuracyGoal -> 3))
``````

## calculus and analysis – Why are the results of integrations of a DiracDelta related integrand and its approximations strikingly different?

Let us consider a double integral, treating it as a functional which corresponds a number to a function,

``````Integrate( Exp(2 x + 3 y)*DiracDelta(x - UnitStep(x + y)), {x, -3/2, 3/2}, {y, -3/2, 3/2})
``````

`1/E^(9/2) - 2 E^2 + E^(13/2)`

In fact, Mathematica calculates the iterated integral, so we also consider another iterated integral

``````Integrate(Exp(2 x + 3 y)*DiracDelta(x - UnitStep(x + y)), {y, -3/2,3/2}, {x, -3/2, 3/2})
``````

`-((2 DiracDelta(0))/E^(9/2))`

As we see, the results are strikingly different. Let us consider approximations of `DiracDelta(x - UnitStep(x + y))`
in the weak topology (Exactly saying, the functions which are associated with that approximations. I’d like to recall that
the $$delta$$-distribution is not associated with any usual function.) ant the limits when those approximations tend to
`DiracDelta(x - UnitStep(x + y))` in the weak topology:

``````Integrate( Exp(2 x + 3 y)*UnitBox((x - UnitStep(x + y))/eps)/eps, {x, -3/2, 3/2},
{y, -3/2, 3/2}, Assumptions -> eps > 0 && eps < 1)
``````

`(-2 E^(7/2) Sinh(eps/2) + 2 E^(9/2) Sinh(eps/2) - Sinh(eps) + E^11 Sinh(eps))/(3 E^(9/2) eps)`

``````Limit(%, eps -> 0, Direction -> "FromAbove")
``````

`(-1 - E^(7/2) + E^(9/2) + E^11)/(3 E^(9/2))`

and

``````Integrate(Exp(2 x + 3 y)*UnitBox((x - UnitStep(x + y))/eps)/eps, {y, -3/2, 3/2},
{x, -3/2, 3/2}, Assumptions -> eps > 0 && eps < 1)
``````

`1/(3 eps) E^(-(9/2) - (3 eps)/2) (-2 E^(7/2 + (3 eps)/2) Sinh(eps/2) + 3 E^(9/2 + (3 eps)/2) Sinh(eps/2) + E^(7/2 + (5 eps)/2) Sinh(eps/2) + E^(7/2 + (7 eps)/2) Sinh(eps/2) + E^(3/2) Sinh(eps) + E^(9/2) Sinh(eps) + E^(13/2) Sinh(eps) - E^(3 eps/2) Sinh(eps) + E^(11 + (3 eps)/2) Sinh(eps) - E^(3/2 + 3 eps) Sinh(eps) - E^(7/2 + 3 eps) Sinh(eps) - E^(13/2 + 3 eps) Sinh(eps) - E^(9/2 + eps/2) Sinh((3 eps)/2) + 2 E^(3/2 + (3 eps)/2) Sinh(eps) Sinh((3 eps)/2) + 2 E^(13/2 + (3 eps)/2) Sinh(eps) Sinh((3 eps)/2))`

``````Limit(%, eps -> 0, Direction -> "FromAbove")
``````

`(-1 - E^(7/2) + E^(9/2) + E^11)/(3 E^(9/2))`.

To be sure, let us also calculate those with non-symmetric approximations:

``````Integrate(Exp(2 x + 3 y)* UnitBox((x - UnitStep(x + y) - Sqrt(eps))/eps)/eps, {x, -3/2,3/2},
{y, -3/2, 3/2}, Assumptions -> eps > 0 && eps < 1);
Limit(%, eps -> 0, Direction -> "FromAbove")
``````

`(-1 - E^(7/2) + E^(9/2) + E^11)/(3 E^(9/2)) `

and

``````eps = 0.005; NIntegrate(Exp(2*x + 3*y)*eps/((x - UnitStep(x + y) - eps)^2 + eps^2)/Pi,
{y, -3/2, 3/2}, {x, -3/2, 3/2})
``````

`223.53`

(compare the latter with `N((-1 - E^(7/2) + E^(9/2) + E^11)/(3 E^(9/2)))` which results in
`221.921` and with `N(1/E^(9/2) - 2 E^2 + E^(13/2))` which results in `650.375`).

It should be noticed the case under consideration is substantially different from
that case:
two iterated integrals of ` Exp(2 x + 3 y)*DiracDelta(x - UnitStep(x + y))` and limits of the integrals of its approximation
in the weak topology `Exp(2 x + 3 y)*UnitBox((x - UnitStep(x + y))/eps)/eps` as `eps` tends to zero from above produce three
different results; if the limits of the integration are slightly perturbed, then the intersection of the support of `DiracDelta(x - UnitStep(x + y))`
with the boundary of the set of the integration does not qualitativily change as

``````ContourPlot(x - UnitStep(x + y) == 0, {x, -2, 2}, {y, -2, 2})
``````

shows.

The question is open: what does the result of the integration of `Exp(2 x + 3 y)*DiracDelta(x - UnitStep(x + y))`
over the square `{x, -3/2, 3/2}, {y, -3/2, 3/2}` mean? The only reference I know where the integration of distributions over
bounded sets is defined is Antosik, P., Mikusinski, J., Sikorski, R. Theory of distributions. The sequential approach. Reprints. (English) Zbl 0267.46028
Amsterdam: Elsevier Scientific Publishing Company; Warszawa: PWN-Polish Scientific Publishers. XIV, 273 p.(1973)
,
but a smooth integral introduced there does not include the integral under consideration.

## numerical integration – Is there a way to speed up Integrate when the integrand contains a product of polynomials each of which having a large degree?

I have integrals of the form
$$int_0^inftymathrm d x; e^{-x^2};;(textrm{Polynomial of degree n})times(textrm{Polynomial of degree m}),$$
where $$n,m$$ can be as large as 300. (If this is helpful, these polynomials are Laguerre polynomials, provided by the function `LaguerreL(...)`).

When I use `Integrate(f(x), {x, 0, Infinity})`, the calculation is too slow for large polynomial degrees.

Is there a way to speed up this calculation? (I tried to use NIntegrate, but unfortunately, since the integrand is highly oscillatory at large $$n,m$$, I’m unable to get reliable results from numerical methods).

## calculus and analysis – Numerical integration strategy for oscillatory integrand

What is the best procedure for calculating the following kind of oscillatory integrals?

``````NIntegrate(f(t) BesselJ(n, t) Sin(t x - n Pi/2), {t, 0, Infinity},
Assumptions -> {-1 < x < 1, n > 0})
``````

where f(t) is bounded at infinity

For instance, for f(t) = 1 there exists a closed form solution for validation of the numerical results, namely

``````W(n_, x_) := -(Sin(n ArcCos(x))/Sqrt(1 - x^2))
``````

I tried three different approaches:

``````W1(n_, x_?NumericQ) :=
NIntegrate(
BesselJ(n, ArcTanh(t)) /(1 - t^2)
Sin(ArcTanh(t) x - n Pi/2), {t, 0, 1},
Method -> {"GlobalAdaptive", "SymbolicProcessing" -> 0});

W2(n_, x_?NumericQ) :=
NIntegrate( BesselJ(n, a) Sin(a x - n Pi/2), {a, 0, Infinity},
Method -> {"ExtrapolatingOscillatory", "SymbolicProcessing" -> 0});

W3(n_, x_?NumericQ) :=
NIntegrate(BesselJ(n, a) Sin(a x - n Pi/2), {a, 0, Infinity},
Method -> {"DoubleExponentialOscillatory",
"SymbolicProcessing" -> 0} );
``````

Here the numerical results:

``````In(1):= Quiet(Timing(W1(3, .3)))

Out(1)= {0.109375, 0.376386}

In(2):= Quiet(Timing(W2(3, .3)))

Out(2)= {2.07813, 0.662773}

In(3):= Quiet(Timing(W3(3, .3)))

Out(3)= {0.015625, 1.14447*10^123}
``````

The exact solution is

``````In(4):= W(3, .3)

Out(4)= 0.64
``````

How can I improve the numerical solution?
Thanks

## calculus and analysis – How does mathematica integrate around points where the integrand evaluates to ”indeterminate?”

I am interested in studying a numerical integral where the integrand, when evaluated naively at some points in the domain of integration, is indeterminate. However, when one carefully takes limits to these points, the integrand turns out to be finite. Since my actual integrand is multi-dimensional and complicated, let me pose my question in the context of a simple example instead. If I ask mathematica to evaluate $$frac{sin x}{x}$$ at $$x=0$$, it says it’s indeterminate (although we know that the actual value as $$x to 0$$ is finite). However, if I evaluate the numerical integral of this function from $$-1 leq x leq 1$$, it gives me the right result, i.e, a value that matches the result of analytic integration.

So my question is: in numerical integration, how does mathematica deal with points where if you ask mathematica to evaluate the integrand, it gives ”indeterminate.” Are these points omitted? If not, how are they accounted for? In the example above, mathematica was somehow able to deal with $$x=0$$ correctly when integrating numerically even though it said that the function evaluated at $$x=0$$ was indeterminate.

## real analysis – Deducing finiteness of the integrand from finiteness of integral

Suppose we know that:
$$int_{(-pi, pi)^d} , g(k) , h^2(k) , dk leqint_{ (-pi, pi)^d} , h(k) , dk < infty.$$
where $$h(k) = sin(frac{k_1}{2})^2 + sum_{i=2}^{d} sin(k_i)^2,$$
where $$k = (k_1, ldots, k_d)$$ and the finiteness of the integral in the RHS holds only when $$d > 2$$ (since $$h(k) = O(1/|k|^2)$$ in the limit as $$|k| rightarrow 0$$.
Is it possible from this to deduce that,
$$int_{(-pi, pi)^d} , g(k) < infty?$$

## real analysis – changing argument of integrand of a periodic function inside an integral

Let $$f: mathbb{R} rightarrow mathbb{R}$$ be a $$2pi$$-periodic function that is integrable on $$(0,2pi)$$. Then for all $$x in mathbb{R}$$:
begin{align}int^{x+2pi}_x f(y) dy =int^{2pi}_0 f(y) dy = int^{2pi}_0 f(y+x) dy. end{align}
I get the first equation, because $$F(x) = int^{x+2pi}_x f(y) dy$$ is a constant function. However, can someone show me why the second equation is true ? I thought that $$f(y) = f(y+x)$$ is only true when $$x = 2pi$$ , since it’s a $$2pi$$-periodic function ?

## numerical integration – NIntegrate blowing up/behaving weirdly at the turning point of the integrand

I’m performing (what should be) a straightforward numerical integration (Fourier transform) of a function with no poles / singularities (at least in a particular parameter regime):

``````<< FourierSeries`

ClearAll("Global`*")

l = 1;
(Epsilon) = 10^-4;

p = 0;
P = 80;

Data1 = ConstantArray(Null, {P, 2});

While(p < P, p++;

w = -10 + 20 p/P;

(Sigma) = (
2 (-2 l^2 r^2 Sin(s/(2 Sqrt(l^-2 - r^2)))^2 + (1 - l^2 r^2) 2 Sinh(
s/(2 Sqrt(l^-2 - r^2)) - I (Epsilon))^2))/l^2;

RF = NInverseFourierTransform(-1/(4 (Pi)^2) E^(-I w s)/(Sigma), s,
w, Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0},
MinRecursion -> 6) // Chop;

Data1((p, 1)) = w;
Data1((p, 2)) = RF;

);

ListLinePlot(Data1, PlotRange -> All,
LabelStyle -> {FontSize -> 16, FontFamily -> "CMU Serif", Black},
PlotLabel -> StringForm("r = ``.", {r // N})) // Print;
``````

When the parameter `r` is small, the numerics seem to work fine/as expected, producing plots like this (`r = 0.1`):

(x-axis is the energy, y-axis is the transition rate of a two-level system – this shows a roughly thermal Planck spectrum).

When increasing `r` beyond some critical point, the numerics seem to blow up, giving a weird answer. For example, `r = 0.8` yields:

(see especially the magnitude of the y-axis).

I plotted the integrand for these two values below. The numerics seem to blow up when the function bifurcates from having one turning point to two turning points:

`r=0.1` (integrand plotted as a function of `s`)

and `r =0.8`

Why does `NIntegrate` not like this seemingly inconspicuous change in behaviour of the integrand? Any help would be greatly appreciated! (If there’s an analytic solution, even better :P)