NDSolve with error that "Power::infy: Infinite expression 1/0.^0.4 encountered."

Question

Here is my code

b = 0.1;
a= 0.6; 
p0 = 1.5;
x0 = 1.3;

eqns = {x'[t] == y[t], 
   y'[t] == p0*x[t]^2 - 
     4 (1 + a b*
         x[t]^2 (x[t]^4 + 2 x[t]^-2 - 3)^(a - 1)) (x[t] - 
        x[t]^(-5))};
ics = {x[0] == x0, y[0] == 0};
sol=NDSolve[Join[eqns,ics],{x,y},{t,0,100}]

It solves with a dynamical system with a periodic solution, but the NDSolve gives an error that "Power::infy: Infinite expression 1/0.^0.4 encountered." Any suggestions would be much appreciated!

Answer 1

You should remove the pseudo-singularity at x[t] == 1 (x[t]-1 is a hidden factor in the numerator and denominator). It's probably better to do this by hand, but I went down the rabbit hole.

Clear[a, b, p0, x0, x, y, t];
eqns = {x'[t] == y[t], 
   y'[t] == 
    p0*x[t]^2 - 
     4 ((x[t] - x[t]^(-5)) + 
        a b*x[t]^2 PowerExpand[(x[t]^4 + 2 x[t]^-2 - 3 // 
              Factor)^(a - 1) (x[t] - x[t]^(-5) // Factor)])} /. (x[
       t] - 1)^p_ :> ((x[t] - 1)^2)^(p/2)
ics = {x[0] == x0, y[0] == 0}

Block[{b = 0.1, a = 0.6, p0 = 1.5, x0 = 1.3},
 sol = NDSolve[Join[eqns, ics], {x, y}, {t, 0, 100}]
 ]

ParametricPlot[Evaluate[{x[t], y[t]} /. First@sol], {t, 0, 100}, AspectRatio -> 1]

"SymplecticPartitionedRungeKutta" seems a better method for this problem but it does not make a great difference. To integrate out to t == 1000, it's about 60% faster than the default method (0.44s vs. 0.72s).

sol = NDSolve[Join[eqns, ics], {x, y}, {t, 0, 1000}, 
  Method -> {"SymplecticPartitionedRungeKutta", 
    "DifferenceOrder" -> 4, "PositionVariables" -> {x[t]}}]

Since the version seems to matter (I cannot reproduce the current answers), I'm using this:

$Version

"12.1.1 for Mac OS X x86 (64-bit) (June 19, 2020)"

Answer 2

Start with rational numbers

b = 0.1 // Rationalize;
a = 0.6 // Rationalize;
p0 = 1.5 // Rationalize;
x0 = 1.3 // Rationalize;
eqns = {x'[t] == y[t],y'[t] ==p0*x[t]^2 -4 (1 + a b*x[t]^2 (x[t]^4 + 2 x[t]^-2 - 3)^(a - 1)) (x[t] -x[t]^(-5))};
ics = {x[0] == x0, y[0] == 0};
sol = NDSolveValue[Join[eqns, ics], {x , y}, {t, 0, 100 } ]


Plot[Through[sol[t]], {t, 0, 100}]

Answer 3

Mathematica makes a difference between exact numbers and inexact ones. The boundary values You entered are inexact numbers.

The best practice is to use as often as possible exact numbers. This comes from the set of paradigms that Mathematica is built on.

There is a rule in Mathematica with exact numbers built-ins that are more powerful. That means inherently mathematical can appear to disappear with the use of exact numbers.

In this case, it means that Mathematica internally does not use really the number entered but an approximation. That suffices to avoid the warning General::infy: from which the waring in the question is derived. The use of Rationalize does the same trick. But it hides the truth behind.

So the most instructive advice on how to deal with these warnings is.

b = 0.1;
a = 0.6;
p0 = 1.5;
x0 = 1.3;
msol = NDSolve[{eqns, ics}, {x, y}, {t, 0, 100}, 
  WorkingPrecision -> MachinePrecision]

ParametricPlot[Evaluate[{x[t], y[t]} /. msol], {t, 0, 1.72}, 
 ColorFunction -> Hue, AspectRatio -> 1]

Between 1.72 and 1.72 is the periodicity of the motion in the x-y-plane.

Mind the small gap on the x-axis. This is due to the fact that 1.72 is slightly smaller than the period.

The set of boundary value will suffice to avoid the warning:

b = 1/10;
a = 3/5;
p0 = 3/2;
x0 = 13/10;

These are rational exact numbers in Mathematica. The exactness is caused by the use of Integers as a numerator and denominator.

b = 1/10;
a = 3/5;
p0 = 3/2;
x0 = 13/10;

eqns = {x'[t] == y[t], 
   y'[t] == 
    p0*x[t]^2 - 
     4*(1 + a*b*x[t]^2 (x[t]^4 + 2 x[t]^-5)^(a - 1)) (x[t] - 
        x[t]^(-5))};
ics = {x[0] == x0, y[0] == 0};
sol = NDSolve[{eqns, ics}, {x, y}, {t, 0, 100}]

Plot[{x[t], y[t]} /. sol, {t, 0, 1.73}]

The domain for which the interpolation functions are calculated shows the high numerical stability of the first order differential-algebraic equation system under question. The further the inexact number is put to zero the smaller the interval will get.

So inexact values for t0 in {t, t0, 100} make a big deal for those interested. Such an example is stable.

b = 0.1;
a = 0.6;
p0 = 1.5;
x0 = 1.3;
msol = NDSolve[{eqns, ics}, {x, y}, {t, 0.0000001, 100}]

But t0=0 will not be. There are many paths in Mathematica to avoid "Power::infy: Infinite expression 1/0.^0.4 encountered." but they are not explained in the documentation of the warning.