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.)

interpolation – How to find the integrand by fitting data – spectrum reconstruction

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

enter image description here

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"})

enter image description here

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

enter image description here

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"})

enter image description here
enter image description here

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}));
   curvaturesList.Boole@Thread(values < (Mu))
   );
(*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})

enter image description here

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;
r = 0.1; (*radial coordinate*)

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):
enter image description here

(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:

enter image description here

(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)

enter image description here

and r =0.8

enter image description here

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)