1
$\begingroup$

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 For inipoin=3.

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!

$\endgroup$
2
$\begingroup$

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}]]

Figure 1 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}]

Figure 2

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.

Figure 3

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]

Fugure 4

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – keanhy14 Jun 7 at 3:33
0
$\begingroup$

@Alex Trounev 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 interval***)

sgpoints = 
  x0 /. NSolve[
    x0 (1 - 113 x0^4 + 2 x0^6) (1 - 58 x0^4 + 2 x0^6) == 0 && x0 > 0, 
    x0];
plist = Table[
   NIntegrate[f[\[Lambda]], {\[Lambda], x, sup[x]}, 
    Exclusions -> sgpoints, Method -> "PrincipalValue"], {x, 
    xlist}](*//Quiet*);
data = Transpose[{xlist, plist}];
ListLinePlot[data] 

The similar warning appear

warning: The integrand has evaluated to Overflow, Indeterminate, or Infinity with boundaries Then, what actually happened?

$\endgroup$
  • $\begingroup$ There is a suspicion that this integral is calculated incorrectly. Therefore, I replaced NIntegrate[] with Sum[] using GaussianQuadratureWeights. $\endgroup$ – Alex Trounev Jun 7 at 7:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.