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

$\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
    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$
    – kcr
    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$
    – kcr
    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
    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$
    – kcr
    Sep 12, 2020 at 20:47

1 Answer 1

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.

$\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$
    – kcr
    Sep 9, 2020 at 20:39
  • 1
    $\begingroup$ No prob, good to put that LyapunovExponents function to the test! $\endgroup$
    – Chris K
    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$
    – kcr
    Sep 10, 2020 at 18:18

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.