Numerical Integration With Many Singularities

Question

Is possible a solution to this Numerical Integration?

The constant values are as follows:

  In[54]:= ClearAll[f1, f, P, c, α, V0, EE, II, k, m, X0, t, x, a]

In[236]:= P = 10.5; 
m = 48.2;
c = 1.;
α = 0.1;
V0 = 10.;
EE = 2.3;
II = 1.;
k = 68.9*10^3;

The function I want to integrate is :

In[130]:= 
f[x_, t_] := 
  1/(4*Pi^2) ((E^((I (V0 ξ + ω)^2)/(2 α ξ)) P Sqrt[2 π])/Sqrt[I α ξ]*
    E^(I*ξ*x)*E^(I*ω*t))/(EE II ξ^4 + k + I*ω*c - m ω^2);

The Numerical integration and the result is this:

    In[288]:= 
WDisp[x_, t_, ξmax_, ωmax_] := 
 NIntegrate[
  1/(4*Pi^2) ((E^((I (V0 ξ + ω)^2)/(2 α ξ)) P Sqrt[2 π])/Sqrt[I α ξ]*
    E^(I*ξ*x)*E^(I*ω*t))/(EE II ξ^4 + k + I*ω*c - m ω^2), {ξ, -ξmax, 
   ξmax}, {ω, -ωmax, ωmax}]

In[289]:= WDisp[0, 0, 20, 20] // Chop

During evaluation of In[289]:= NIntegrate::inumexpr: Expression {((0. +5. I) (10. ξ+ω)^2)/ξ} derived from integrand (0.666683 E^(((0. +5. I) (<<1>>+ω)^2)/ξ))/(Sqrt[(0. +0.1 I) ξ] (68900. +2.3 ξ^4+(0. +1. I) ω-48.2 ω^2)) is not numerical at {ξ,ω} = {0.,-19.8313}. >>

During evaluation of In[289]:= NIntegrate::mtdfb: Numerical integration with LevinRule failed. The integration continues with Method -> MultiDimensionalRule. >>

During evaluation of In[289]:= NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small. >>

During evaluation of In[289]:= NIntegrate::eincr: The global error of the strategy GlobalAdaptive has increased more than 2000 times. The global error is expected to decrease monotonically after a number of integrand evaluations. Suspect one of the following: the working precision is insufficient for the specified precision goal; the integrand is highly oscillatory or it is not a (piecewise) smooth function; or the true value of the integral is 0. Increasing the value of the GlobalAdaptive option MaxErrorIncreases might lead to a convergent numerical integration. NIntegrate obtained 0.0000941872 +1.05263*10^-7 I and 0.010442730338460396` for the integral and error estimates. >>

Out[289]= 0.0000941872 + 1.05263*10^-7 I

The Plot of The function Shows that it has high singularities:

Thank you in advance

Answer 1

The problem is, you get division by 0 at [Xi]==0.

Since NIntegrate starts integration a litte bit near the integration limits, split the [Xi] integration at [Xi]==0.

WDisp[x_, t_, \[Xi]max_, \[Omega]max_] := 
  NIntegrate[
    1/(4*Pi^2) ((E^((I (V0 \[Xi] + \[Omega])^2)/(2 \[Alpha] \[Xi])) P \
    Sqrt[2 \[Pi]])/Sqrt[I \[Alpha] \[Xi]]*E^(I*\[Xi]*x)*
  E^(I*\[Omega]*t))/(EE II \[Xi]^4 + k + I*\[Omega]*c - 
  m \[Omega]^2), {\[Xi], -\[Xi]max, 
   0, \[Xi]max}, {\[Omega], -\[Omega]max, \[Omega]max}, 
   MaxRecursion -> 100]

WDisp[0, 0, 20, 20]

Although convergence is very slowly because the integrand is highly oscillatory, you get a result without serious error messages

(*   0.000107936- 2.65831*10^-16 I    *)