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

Question

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!

Answer 1

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]