3
$\begingroup$

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!

Improve this question
$\endgroup$
3
$\begingroup$

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)"
Improve this answer
$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$ – keanhy14 Aug 6 '20 at 6:29
  • 1
    $\begingroup$ @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.. $\endgroup$ – Michael E2 Aug 6 '20 at 12:27
1
$\begingroup$

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

enter image description here

Improve this answer
$\endgroup$
8
  • 1
    $\begingroup$ What version are you using? With v12.1.1 on a Mac, the solutions blow up at around x == 7.2 $\endgroup$ – Bob Hanlon Aug 5 '20 at 14:12
  • $\begingroup$ @ BobHanlon My version is v12 Windows 64. Try the option MaxStepSize-> ... $\endgroup$ – Ulrich Neumann Aug 5 '20 at 14:15
  • $\begingroup$ With or without the MaxStepSize the InterpolatingFunction domains are shown as {{0., 7.2}} $\endgroup$ – Bob Hanlon Aug 5 '20 at 14:22
  • $\begingroup$ @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. $\endgroup$ – keanhy14 Aug 5 '20 at 15:39
  • $\begingroup$ @keanhy14 Don't know what's the reason for this unstability. Why do you need such a huge time range? $\endgroup$ – Ulrich Neumann Aug 6 '20 at 7:04
1
$\begingroup$

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]

WorkingPrecision->MachinePrecision

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

ParametricPlot of this first-order algebraic differential equation

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

enter image description here

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

Plot of the solution over the period using the rules

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

solution

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.

Improve this answer
$\endgroup$
2
  • 1
    $\begingroup$ 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? $\endgroup$ – keanhy14 Aug 5 '20 at 16:32
  • $\begingroup$ Does using the exact number slow down the numerical pde process? $\endgroup$ – Nam Nguyen Oct 18 '20 at 6:34

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.