# NIntegrate warning: The integrand has evaluated to Overflow, Indeterminate, or Infinity with boundaries

The code

Clear["Global*"]
a = 110.; b = 55.; d = 1.; m1 = 18.; m2 = 42.; m = m2/m1;
n = 10;
deltap = .8;
inipoint = 3.;
tlength = 1.14;
w[lam_, xi_] := -((m1*a)/2) Log[1 -(lam^(-4) + 2*lam^2 - 3)/a] - (m2*b)/2 Log[1 - (lam^-4*xi^4 + 2 lam^2*xi^-2 - 3)/b]
dw[lam_, xi_] := D[w[lam, xi], lam]
f[lam_, xi_] := dw[lam, xi]/(1 - lam^3)
sup[x_] := ((d + x^3)/(1 + d))^(1/3)

intf[x_Real, xi_] :=
NIntegrate[f[lam, xi], {lam, x, sup[x]},
Exclusions -> (x0 /. NSolve[Denominator[Together@f[x0, xi]] == 0, x0]),
Method -> "PrincipalValue"]
eq1 := x[t]*x''[t] (1 - (1 + d/x[t]^3)^(-1/3)) +
1/2 x'[t]^2 (3 - d/x[t]^3 (1 + d/x[t]^3)^(-4/3) -
3 (1 + d/x[t]^3)^(-1/3)) + intf[x[t], xi[t]] == deltap;
eq2 := xi'[t]/
xi[t] == (m (x[t]^2*xi[t]^-2 -
x[t]^-4*xi[t]^4))/(3 n*(1 - (x[t]^-4*xi[t]^4 +
2 x[t]^2*xi[t]^-2 - 3)/b));

sol = Timing@
NDSolve[{eq1, eq2, xi[0] == 1, x'[0] == 0,
x[0] == inipoint}, {x[t], xi[t]}, {t, 0, tlength}];


when the "inipoint" is set to be 3, it gives the result

However, when "inipoint=4.", the warning is NIntegrate::inumri: The integrand has evaluated to Overflow, Indeterminate, or Infinity for all sampling points in the region with boundaries {{0.,1.68335*10^-22}}.

Someone has ideas about it? Any suggestions would be much appreciated!

First, it is possible to lower the integration time by three orders, using GaussianQuadratureWeights instead of NIntegrate. Check that the following code is executed in 0.046875 sec and has the same output like the original code that runs 51.0781 sec (on my laptop):

Clear["Global*"]
\[Alpha] = 110.; \[Beta] = 55.; \[Delta] = 1.; \[Mu]1 = 18.; \[Mu]2 = \
42.; \[Mu] = \[Mu]2/\[Mu]1;
\[Eta]b = 10;
deltap = .8;
inipoint = 3.;
tlength = 1.14;
w[\[Lambda]_, \[Xi]_] := -((\[Mu]1*\[Alpha])/2) Log[
1 - (\[Lambda]^(-4) + 2*\[Lambda]^2 -
3)/\[Alpha]] - (\[Mu]2*\[Beta])/2 Log[
1 - (\[Lambda]^-4*\[Xi]^4 + 2 \[Lambda]^2*\[Xi]^-2 - 3)/\[Beta]]
dw[\[Lambda]_, \[Xi]_] := D[w[\[Lambda], \[Xi]], \[Lambda]]
f[\[Lambda]_, \[Xi]_] := dw[\[Lambda], \[Xi]]/(1 - \[Lambda]^3)
sup[x_] := ((\[Delta] + x^3)/(1 + \[Delta]))^(1/3)
Get["NumericalDifferentialEquationAnalysis"];
np = 11; points = weights = Table[Null, {np}];

intf[x0_, \[Xi]0_] :=
Block[{y = x0, \[Xi]1 = \[Xi]0},
Do[points[[i]] =
GaussianQuadratureWeights[np, y, sup[y]][[i, 1]], {i, 1, np}];
Do[weights[[i]] =
GaussianQuadratureWeights[np, y, sup[y]][[i, 2]], {i, 1, np}];
int = Sum[(f[\[Lambda], \[Xi]1] /. \[Lambda] -> points[[i]])*
weights[[i]], {i, 1, np}]; int]
eq1 := x''[t] + (1/
2 x'[t]^2 (3 - \[Delta]/
x[t]^3 (1 + \[Delta]/x[t]^3)^(-4/3) -
3 (1 + \[Delta]/x[t]^3)^(-1/3)) + intf[x[t], \[Xi][t]] -
deltap)/x[t]/(1 - (1 + \[Delta]/x[t]^3)^(-1/3)) == 0;
eq2 := \[Xi]'[
t] == \[Xi][
t]*(\[Mu] (x[t]^2*\[Xi][t]^-2 -
x[t]^-4*\[Xi][t]^4))/(3 \[Eta]b*(1 - (x[t]^-4*\[Xi][t]^4 +
2 x[t]^2*\[Xi][t]^-2 - 3)/\[Beta]));

sol = Timing@
NDSolve[{eq1, eq2, \[Xi][0] == 1, x'[0] == 0,
x[0] == inipoint}, {x[t], \[Xi][t]}, {t, 0, tlength}]


Figure 1 shows x[t] (1) and $$\xi (t)$$ (2) for inipoint=3.

Plot[Evaluate[{x[t], \[Xi][t]} /. Last[sol]], {t, 0, tlength},
AxesLabel -> Automatic, PlotLegends -> Automatic,
PlotLabel -> Row[{"inipoint = ", inipoint}]]


Construct a parametric function

pfun = ParametricNDSolveValue[{eq1, eq2, \[Xi][0] == 1, x'[0] == 0,
x[0] == p}, {x[tlength], \[Xi][tlength]}, {t, 0, tlength}, {p}]


Then we have

pfun[3]

(*Out[]= {1.43111, 0.826233}*)


Figure 2 shows the dependence of functions {x[tlength], \[Xi][tlength]} on the parameter p. We see that when p>3.5095592, instability is observed.

Plot[pfun[p], {p, 2, 3.5}]


Compare the solution with p=3.5095592 and with p=3.5095595on fig.3. We see that the instability develops during the first rebound at x[t]-> 0.

To make sure that this is a numerical instability, we use an explicit method with a very small step. Figure 4 shows that it is possible to pass the first bounce at inipoint=3.6.

eq11 = {x'[t] == y[t],
y'[t] + (1/
2 y[t]^2 (3 - \[Delta]/
x[t]^3 (1 + \[Delta]/x[t]^3)^(-4/3) -
3 (1 + \[Delta]/x[t]^3)^(-1/3)) + intf[x[t], \[Xi][t]] -
deltap)/x[t]/(1 - (1 + \[Delta]/x[t]^3)^(-1/3)) ==
0}; eq21 = \[Xi]'[
t] == \[Xi][
t]*(\[Mu] (x[t]^2*\[Xi][t]^-2 -
x[t]^-4*\[Xi][t]^4))/(3 \[Eta]b*(1 - (x[t]^-4*\[Xi][t]^4 +
2 x[t]^2*\[Xi][t]^-2 - 3)/\[Beta]));

sol11 = Timing@
NDSolve[{eq11, eq2, \[Xi][0] == 1, y[0] == 0,
x[0] == 3.6}, {x[t], \[Xi][t], y[t]}, {t, 0, tlength},
Method -> "ExplicitEuler", MaxStepSize -> 1/5000000,
MaxSteps -> 10^7]

Plot[Evaluate[{x[t], \[Xi][t]} /. Last[sol11]], {t, 0, tlength},
AxesLabel -> Automatic, PlotLegends -> Automatic,
PlotLabel -> Row[{"inipoint = ", 3.6}], PlotRange -> All]


• That's a great job you've done! But when I use the NIntegrate function and change the interval xlist = Range[5, 8, 0.002] from 0.002 to 0.001 or other smaller values for the code below, Clear["Global*"] [Alpha] = 110.; [Beta] = 55.; [Delta] = 1; [Mu] = 42/18; f[[Lambda]_] := ((2*([Lambda] - [Lambda]^(-5)))/(1 - \ ([Lambda]^(-4) + 2*[Lambda]^2 - 3)/[Alpha]) + (2*[Mu] ([Lambda] - [Lambda]^(-5)))/(1 - ([Lambda]^(-4) + 2*[Lambda]^2 - 3)/[Beta]))/(1 - [Lambda]^3) sup[x_] := (([Delta] + x^3)/(1 + [Delta]))^(1/3) xlist = Range[5, 8, 0.002]; (***the change of inte Commented Jun 7, 2019 at 3:29
• Thank you for your help! I have the similar problem with codes listed below in the "answer part", would you mind to give some more suggestions? Commented Jun 7, 2019 at 3:33
• There is a suspicion that this integral is calculated incorrectly. Therefore, I replaced NIntegrate[] with Sum[] using GaussianQuadratureWeights. Commented Jun 7, 2019 at 7:25