# Numerical Integration With Many Singularities

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:

• I don't think the integral is convergent, but I am an engineer, not a mathematician. You might ask the existence question on the mathematics SE. – MikeY Jun 16 '17 at 13:05
• Thank you, I can see that it does not converge – Edmond Muho Jun 16 '17 at 13:52

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    *)
`
• Thank you very much for your help, I appreciated – Edmond Muho Jun 16 '17 at 13:08