Erreurs dans le code de l'exposant Lyapunov

Question

Grâce à la réponse de @Chris K, je pense avoir ré-exprimé la question correctement.

J'ai les équations de mouvement suivantes

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

Je peux résoudre numériquement ce qui précède en utilisant

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]

Et d'après l'intrigue, il est clair que le système est chaotique

Ensuite, j'ai essayé de calculer l'exposant de Lyapunov en utilisant le code affiché ici . Bien entendu, pour ce faire, je ré-exprime les équations dans un formalisme du premier ordre. Le code est

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

et alors

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

Et dès que j'exécute la commande finale, Mma produit beaucoup d'erreurs.

Après le commentaire instructif de ChrisK, j'ai essayé d'exécuter les équations du premier ordre en utilisant NDSolve. Les oeuvres suivantes

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]

Cependant, la commande

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

Renvoie à nouveau des erreurs.

Des pensées? J'ai lu qu'il y avait une modification récente du code posté ici mais j'ai utilisé la dernière mise à jour. Je suppose que la vraie question est de savoir comment passer les méthodes pour NDSolve dans le code d'origine où les commandes sont de la forme NDSolveProcessEquations andNDSolveIterate

Answer 1

Je pense que le problème réside dans vos équations. Dans les deux cas, NDSolves'arrête avec une NDSolve::ndszerreur (la taille du pas est effectivement zéro; singularité ou système rigide suspecté). Si vous regardez la solution que vous avez sltn, vous verrez qu'elle ne couvre pas toute la gamme des$\tau$:

La variation rapide que vous voyez dans votre tracé au-delà de ce point provient de Sin[z[\[Tau]]son application à l'extrapolation de l' z[\[Tau]]extérieur de la plage où elle est définie.