3
$\begingroup$

Thanks to the answer by @Chris K I think I have re-expressed the question properly.

I have the following equations of motion

eqnx = (1/4)*(-((α1*α2*μ*Sin[2*z[τ]]^2)/(α*β*(α*Cos[z[τ]]^2 + β*Sin[z[τ]]^2 - x[τ])^2)) - 
     (α1^2*Sin[z[τ]]^2*(3*α^3 + 2*α^2*β + 3*α*β^2 + 4*μ + 4*(α^3 - α*β^2 - μ)*Cos[2*z[τ]] + α^3*Cos[4*z[τ]] - 2*α^2*β*Cos[4*z[τ]] + 
        α*β^2*Cos[4*z[τ]] - 8*α*(α + β + (α - β)*Cos[2*z[τ]])*x[τ] + 8*α*x[τ]^2))/(2*α^2*(α*Cos[z[τ]]^2 + β*Sin[z[τ]]^2 - x[τ])^2) - 
     (α2^2*Cos[z[τ]]^2*(3*α^2*β + 2*α*β^2 + 3*β^3 + 4*μ + 4*(α^2*β - β^3 + μ)*Cos[2*z[τ]] + α^2*β*Cos[4*z[τ]] - 2*α*β^2*Cos[4*z[τ]] + 
        β^3*Cos[4*z[τ]] - 8*β*(α + β + (α - β)*Cos[2*z[τ]])*x[τ] + 8*β*x[τ]^2))/(2*β^2*(α*Cos[z[τ]]^2 + β*Sin[z[τ]]^2 - x[τ])^2) + 
     Derivative[1][x][τ]^2/(-μ + (α - x[τ])*(β - x[τ])*x[τ]) - ((α*Cos[z[τ]]^2 + β*Sin[z[τ]]^2 - x[τ])*(α*β - 2*(α + β)*x[τ] + 3*x[τ]^2)*
       Derivative[1][x][τ]^2)/(μ - α*β*x[τ] + (α + β)*x[τ]^2 - x[τ]^3)^2 + (4*Derivative[1][z][τ]^2)/(α*Cos[z[τ]]^2 + β*Sin[z[τ]]^2) + 
     (2*Derivative[1][x][τ]*(Derivative[1][x][τ] + (α - β)*Sin[2*z[τ]]*Derivative[1][z][τ]))/(μ - α*β*x[τ] + (α + β)*x[τ]^2 - x[τ]^3) + 
     (2*(α*Cos[z[τ]]^2 + β*Sin[z[τ]]^2 - x[τ])*Derivative[2][x][τ])/(-μ + (α - x[τ])*(β - x[τ])*x[τ])); 
eqnz = (1/4)*(-((α1*α2*(α - β)*μ*Sin[2*z[τ]]^3)/(α*β*(α*Cos[z[τ]]^2 + β*Sin[z[τ]]^2 - x[τ])^2)) + (8*α2^2*μ*Cos[z[τ]]^3*Sin[z[τ]]*(β - x[τ]))/
      (β^2*(α*Cos[z[τ]]^2 + β*Sin[z[τ]]^2 - x[τ])^2) - (2*α1*α2*μ*Sin[4*z[τ]])/(α*β*(α*Cos[z[τ]]^2 + β*Sin[z[τ]]^2 - x[τ])) + 
     (8*α1^2*μ*Cos[z[τ]]*Sin[z[τ]]^3*(-α + x[τ]))/(α^2*(α*Cos[z[τ]]^2 + β*Sin[z[τ]]^2 - x[τ])^2) - 
     (1/β^2)*8*α2^2*Cos[z[τ]]*Sin[z[τ]]*(Cos[z[τ]]^2*(β - x[τ])^2 + (Sin[z[τ]]^2*(α*Cos[z[τ]]^2 + β*Sin[z[τ]]^2)*(β - x[τ])^2 + 
         Cos[z[τ]]^2*(-μ + (α - x[τ])*(β - x[τ])*x[τ]))/(α*Cos[z[τ]]^2 + β*Sin[z[τ]]^2 - x[τ])) + 
     (1/α^2)*8*α1^2*Cos[z[τ]]*Sin[z[τ]]*(Sin[z[τ]]^2*(α - x[τ])^2 + (Cos[z[τ]]^2*(α*Cos[z[τ]]^2 + β*Sin[z[τ]]^2)*(α - x[τ])^2 + 
         Sin[z[τ]]^2*(-μ + (α - x[τ])*(β - x[τ])*x[τ]))/(α*Cos[z[τ]]^2 + β*Sin[z[τ]]^2 - x[τ])) + ((α - β)*Sin[2*z[τ]]*Derivative[1][x][τ]^2)/
      (-μ + α*β*x[τ] - (α + β)*x[τ]^2 + x[τ]^3) + (4*(α - β)*Sin[2*z[τ]]*Derivative[1][z][τ]^2)/(α*Cos[z[τ]]^2 + β*Sin[z[τ]]^2) + 
     (2*(α - β)*Sin[2*z[τ]]*(α + β + (α - β)*Cos[2*z[τ]] - 2*x[τ])*Derivative[1][z][τ]^2)/(α*Cos[z[τ]]^2 + β*Sin[z[τ]]^2)^2 - 
     (8*Derivative[1][z][τ]*(Derivative[1][x][τ] + (α - β)*Sin[2*z[τ]]*Derivative[1][z][τ]))/(α*Cos[z[τ]]^2 + β*Sin[z[τ]]^2) + 
     (8*(α*Cos[z[τ]]^2 + β*Sin[z[τ]]^2 - x[τ])*Derivative[2][z][τ])/(α*Cos[z[τ]]^2 + β*Sin[z[τ]]^2));

I can solve numerically the above by using

intl = 0; 
lim = 10^4; 
x0 = 2.; 
rule = {α1 -> 1, α2 -> 2, μ -> 1, α -> 2.85383, β -> 3.18783}; 
sltn = First[NDSolve[{{(eqnx /. rule) == 0, (eqnz /. rule) == 0}, x[intl] == x0, Derivative[1][x][intl] == 0, z[intl] == 0.403, Derivative[1][z][intl] == 0.1}, z, 
     {τ, intl, lim}, Method -> {"BDF"}]]; 
z[τ] /. sltn
plot5 = Show[Plot[Sin[z[τ] /. sltn], {τ, intl, lim}, BaseStyle -> {17, FontFamily -> "Times New Roman"}, AxesLabel -> {"τ", "Sin(θ(τ))"}, AxesStyle -> Thick, 
    PlotRange -> {{0, lim}, {-1.05, 1.05}}, PlotStyle -> Red], ImageSize -> Large]

And from the plot it is clear that the system is chaotic enter image description here

Then I tried to compute the Lyapunov exponent by using the code posted here. Of course, in order to do so, I re-express the equations in a first-order formalism. The code is

eqnxx = {Derivative[1][x][t] == y[t], Derivative[1][y][t] == 
     (y[t] - 0.334*Sin[2*z[t]]*g[t])/(2.85*Cos[z[t]]^2 + 3.19*Sin[z[t]]^2 - x[t]) - 
      (2*(-1 + (2.85 - x[t])*(3.19 - x[t])*x[t])*g[t]^2)/((2.85*Cos[z[t]]^2 + 3.19*Sin[z[t]]^2)*
        (2.85*Cos[z[t]]^2 + 3.19*Sin[z[t]]^2 - x[t])) - 0.5*(1/(2.85*Cos[z[t]]^2 + 3.19*Sin[z[t]]^2 - x[t]) - 
        (9.1 - 12.1*x[t] + 3*x[t]^2)/(-1 + (2.85 - x[t])*(3.19 - x[t])*x[t]))*y[t]^2 + 
      (0.25*(-1 + (2.85 - x[t])*(3.19 - x[t])*x[t])*(0.44*Sin[2*z[t]]^2 + 0.394*Cos[z[t]]^2*
          (237 - 21.7*Cos[2*z[t]] + 0.356*Cos[4*z[t]] - 25.5*(6.04 - 0.334*Cos[2*z[t]] - x[t])*x[t]) + 
         0.123*Sin[z[t]]^2*(213 - 27*Cos[2*z[t]] + 0.318*Cos[4*z[t]] - 22.8*(6.04 - 0.334*Cos[2*z[t]] - x[t])*x[t])))/
       (2.85*Cos[z[t]]^2 + 3.19*Sin[z[t]]^2 - x[t])^3, Derivative[1][z][t] == g[t], 
    Derivative[1][g][t] == (0.0418*Sin[2*z[t]]*(y[t]^2/(-1 + (2.85 - x[t])*(3.19 - x[t])*x[t]) - 
         (4*x[t]*g[t]^2)/(2.85*Cos[z[t]]^2 + 3.19*Sin[z[t]]^2)^2))/(1 - (2*x[t])/(6.04 - 0.334*Cos[2*z[t]])) + 
      (2*y[t]*g[t])/(6.04 - 0.334*Cos[2*z[t]] - 2*x[t]) - (1.5*Sin[2*z[t]])/(1 - (2*x[t])/(6.04 - 0.334*Cos[2*z[t]])) - 
      (0.211 - 0.22*(-1 + (4*(2.85 - x[t])*(3.19 - x[t]))/(6.04 - 0.334*Cos[2*z[t]] - 2*x[t])^2) + 
       (0.491*(2.85 - x[t])^2)/(6.04 - 0.334*Cos[2*z[t]] - 2*x[t])^2 + 
       4*(-0.432 + (0.394*(3.19 - x[t])^2)/(6.04 - 0.334*Cos[2*z[t]] - 2*x[t])^2 + 0.105*x[t]) - 0.117*x[t])};

and then

LyapunovExponents[eqnxx, {x -> 2, y -> 0, z -> 0.403, g -> 0.1}, 
 ShowPlot -> True]

And as soon as I run the final command Mma produces a lot of errors.

After the instructive comment by ChrisK I tried to run the first-order equations using NDSolve. The following works

x0 = 2.3;
aa = First[
  NDSolve[{eqnxx, x[intl] == x0, y[intl] == 0.1, z[intl] == 0, 
    g[intl] == 0}, z, {t, intl, lim}, 
   Method -> {"StiffnessSwitching", "NonstiffTest" -> False, 
     Method -> {"ExplicitRungeKutta", Automatic}}, AccuracyGoal -> 5, 
   PrecisionGoal -> 5], MaxSteps -> Infinity]

However, the command

LyapunovExponents[eqnxx, {x -> 2.3, y -> 0.1, z -> 0, g -> 0}, 
 ShowPlot -> True]

Again returns errors.

Any thoughts? I read that there was a recent modification to the code posted here but I have used the latest update. I guess the true question is how to pass the Methods for NDSolve in the original code where the commands are of the form NDSolveProcessEquationsandNDSolveIterate

Improve this question
$\endgroup$
5
  • $\begingroup$ I tried running your new code, but it didn't work (maybe transcribed wrong because there's a mix of \[Theta] and z in the first two code blocks). Running NDSolve on the bottom code (eqnxx) gives me an Power::infy -- Infinite expression 1/0. encountered error with the default solver and an NDSolve::ndsz -- step size is effectively zero; singularity or stiff system suspected error with Method -> {"BDF"}. $\endgroup$
    – Chris K
    Commented Sep 11, 2020 at 2:05
  • $\begingroup$ If there is a \theta it is a transcription error. I cannot find it now. It should be a z. I have not tried to perform the NDSolve for the eqnxx but I believe you, which I think it is peculiar as eqnxx is the same as the first equations but written in first order formalism $\endgroup$
    – user49048
    Commented Sep 11, 2020 at 9:00
  • $\begingroup$ @ChrisK thanks for your time. I fixed the typo you pointed out. I also now provide an example where NDSolve with some methods manage to solve eqnxx but still the Lyapunov exponent returns errors. I phrased in the main post my true question. How do I pass the NDSolve methods and Precision and so on in the original code that you wrote where you used NDSolve`ProcessEquations, etc? $\endgroup$
    – user49048
    Commented Sep 11, 2020 at 9:50
  • $\begingroup$ You can pass options to the NDSolve inside LyapunovExponents with NDSolveOpts->{your options here}. Unfortunately the options you found to work on the system alone don't work here (Power::infy error). The LyapunovExponents code adds some linearized equations to your equations, so it's solving a larger system. You could experiment to find other NDSolve options that might work, but there's still the risk of division by zero lurking in your equations for y'[t] and g'[t]. It feels kind of touchy to me. $\endgroup$
    – Chris K
    Commented Sep 12, 2020 at 20:43
  • $\begingroup$ @ChrisK Thanks again for taking time to explain the commands in your original code. I have been playing all day today so I think I have a better feel now. It seems that you are right though and I might end up having to consider the dynamics after expanding the equations near $0$ and/or infinity and see if I can get an interesting behaviour. Anyway, thanks for the comments and the help $\endgroup$
    – user49048
    Commented Sep 12, 2020 at 20:47

1 Answer 1

Reset to default
4
$\begingroup$

I think the problem is in your equations. In both cases, NDSolve stops with an NDSolve::ndsz (step size is effectively zero; singularity or stiff system suspected) error. If you look at the solution you have in sltn, you'll see it does not cover the whole range of $\tau$: enter image description here

The rapid variation you see in your plot beyond that point comes from Sin[z[\[Tau]] being applied to the extrapolation of z[\[Tau]] outside the range where it is defined.

Improve this answer
$\endgroup$
3
  • $\begingroup$ Ok. Thanks a lot for that. I just wanted to make sure that it was my mistake and not some compatibility issue. $\endgroup$
    – user49048
    Commented Sep 9, 2020 at 20:39
  • 1
    $\begingroup$ No prob, good to put that LyapunovExponents function to the test! $\endgroup$
    – Chris K
    Commented Sep 9, 2020 at 20:40
  • $\begingroup$ Can you have a look at the example I gave after your answer? I managed to solve the equations of motion without errors using the BDF method but then the Lyapunov exponent still spits errors. Do you have any recommendations? $\endgroup$
    – user49048
    Commented Sep 10, 2020 at 18:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.