Fehler im Lyapunov-Exponentencode

Question

Dank der Antwort von @Chris KI denke ich, dass ich die Frage richtig ausgedrückt habe.

Ich habe die folgenden Bewegungsgleichungen

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

Ich kann das oben genannte numerisch lösen, indem ich benutze

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]

Und aus der Handlung geht hervor, dass das System chaotisch ist

Dann habe ich versucht , mit Hilfe des Code geschrieben , die Lyapunov Exponenten berechnen hier . Dazu drücke ich natürlich die Gleichungen in einem Formalismus erster Ordnung erneut aus. Der Code ist

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

und dann

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

Und sobald ich den letzten Befehl ausführe, erzeugt Mma viele Fehler.

Nach dem lehrreichen Kommentar von ChrisK habe ich versucht, die Gleichungen erster Ordnung mit NDSolve auszuführen. Folgendes funktioniert

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]

Allerdings der Befehl

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

Gibt wieder Fehler zurück.

Irgendwelche Gedanken? Ich habe gelesen, dass der hier veröffentlichte Code kürzlich geändert wurde , aber ich habe das neueste Update verwendet. Ich denke, die wahre Frage ist, wie die Methoden für NDSolve im Originalcode übergeben werden, wobei die Befehle die Form NDSolveProcessEquations andNDSolve habenIterate

Answer 1

Ich denke, das Problem liegt in Ihren Gleichungen. In beiden Fällen wird NDSolvemit einem NDSolve::ndszFehler (Schrittgröße ist praktisch Null; Singularität oder steifes System vermutet) gestoppt. Wenn Sie sich die Lösung ansehen, in der Sie sich befinden sltn, werden Sie feststellen, dass sie nicht den gesamten Bereich von abdeckt$\tau$::

Die schnelle Variation, die Sie in Ihrem Diagramm über diesen Punkt hinaus sehen, ergibt sich aus Sin[z[\[Tau]]der Extrapolation z[\[Tau]]außerhalb des Bereichs, in dem sie definiert ist.