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

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!

## 3 Answers

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)"

• Thanks for your excellent answer! The problem of pseudo singularity is well solved and It's really what I want. Just one problem with your code, as in the document, The transformations made by PowerExpand are correct in general only if c is an integer or a and b are positive real numbers. So, the equations in your code will lead to a different solution. I just remove PowerExpand and it still works well. Aug 6, 2020 at 6:29
• @keanhy14 You can check by hand that PowerExpand does nothing illegal here. Removing PowerExpand does not remove the false singularity, which you can check with eqns /. x[t] -> 1.. With or without changing the original eqns, it should be unlikely that the function is evaluated at x[t] = 1., but somehow it usually hits on x[t] = 1., probably because of the many orbits. It misses it, however, with the original eqns with the option StartingStepSize -> 1/99.. Aug 6, 2020 at 12:27

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


• What version are you using? With v12.1.1 on a Mac, the solutions blow up at around x == 7.2 Aug 5, 2020 at 14:12
• @ BobHanlon My version is v12 Windows 64. Try the option MaxStepSize-> ... Aug 5, 2020 at 14:15
• With or without the MaxStepSize the InterpolatingFunction domains are shown as {{0., 7.2}} Aug 5, 2020 at 14:22
• @UlrichNeumann When t is increased 1000, i.e., {t,0,1000}, it doesn't work again. I think it is due to this part of the equation 4 (1 + a b*x[t]^2 (x[t]^4 + 2 x[t]^-2 - 3)^(a - 1)) (x[t] -x[t]^(-5)) , and x[t] can not be 1. Aug 5, 2020 at 15:39
• @keanhy14 Don't know what's the reason for this unstability. Why do you need such a huge time range? Aug 6, 2020 at 7:04

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.

• Thanks for your comments, it‘s instructive！There is something different in your code for the equation (x[t]^4 + 2 x[t]^-5)^(a - 1), and it should be (x[t]^4 + 2 x[t]^-2 - 3)^(a - 1). Moreover, when time range is {t,0,1000}, this method doesn't work. I wonder if there is stiffness in the equation due to this part (x[t]^4 + 2 x[t]^-2 - 3)^(a - 1) when x=1? Aug 5, 2020 at 16:32
• Does using the exact number slow down the numerical pde process? Oct 18, 2020 at 6:34