<< LinearAlgebra`Orthogonalization`

<< Graphics`MultipleListPlot`

SetPendParam[driveforce_, drivefreq_, dampen_, phase_] := 
        Module[{q = driveforce, b = dampen, p = phase, 
            w = drivefreq}, 
f1[t_, x_, y_, z_] = 
            y; 
f2[t_, x_, y_, z_] = 
            q * Cos[z + p] - Sin[x] - b * y; 
f3[t_, 
              x_, y_, z_] = 
            w; 
F[t_, {x_, y_, z_}] = {f1[t, x, y, z], 
              f2[t, x, y, z], f3[t, x, y, z]};];

Runge[initt_, initx_, inity_, stepsper_, 
          numb_] := 
Module[{it = initt, j, 
            stepper = stepsper, num = numb}, 
period = 
            2.0*N[π, 20]*(3/2); 
h = 
            (period/stepper); 
third  =  
            Table[{}, {num}]; 
T = 
            Table[(it + (j - 1)* h), {j, 1, 
                stepper}]; 
R = {initx, inity, 
              it}; 
Do[{For[j = 1, 
                j < stepper + 1, j++, 
k1 = 
                  h*F[T[[j]], R]; 
k2 = 
                  h*F[T[[j]] + h/2, R + k1/2]; 
k3 = 
                  h*F[T[[j]] + h/2, R + k2/2]; 
k4 = 
                  h*F[T[[j]] + h, R + k3]; 
Rnew = 
                  R + 1/6 (k1 + 2 k2 + 2 k3 + 
                          k4); 
R = 
                  Rnew;
]; 
second = 
                R; 
third[[k]] = 
                second}
, {k, 1, 
              num}]
];

SetPendParam[0.9, 2/3, 0.5, 0]
fourth1  =  Table[{}, {1}];
Runge[0, 0, 1, 1000, 1];
Do[fourth1[[k]]  =  {Mod[third[[k]][[1]] - π, 
            2  π], third[[k]][[2]]}, {k, 1, 
        1}]
fourth2  =  Table[{}, {1}];
Runge[0, 0.001, 1, 1000, 1];
Do[fourth2[[k]]  =  {Mod[third[[k]][[1]] - π, 
            2  π], third[[k]][[2]]}, {k, 1, 
        1}]
temp1  =  Last[fourth1];
temp2  =  Last[fourth2];
dθτ  =  
        temp2[[1]] - temp1[[1]];
dωτ  =  
        temp2[[2]] - temp1[[2]];
q1θ  =  dθ(τ/0.001);
q1ω  =  dω(τ/0.001);
fourth1  =  Table[{}, {1}];
Runge[0, 0, 1, 1000, 1];
Do[fourth1[[k]]  =  {Mod[third[[k]][[1]] - π, 
            2  π], third[[k]][[2]]}, {k, 1, 1}]
fourth2  =  Table[{}, {1}];
Runge[0, 0, 1.001, 1000, 1];
Do[fourth2[[k]]  =  {Mod[third[[k]][[1]] - π, 
            2  π], third[[k]][[2]]}, {k, 1, 1}]
temp1  =  Last[fourth1];
temp2  =  Last[fourth2];
dθτ  =  temp2[[1]] - temp1[[1]];
dωτ  =  temp2[[2]] - temp1[[2]];
q2θ  =  dθ(τ/0.001);
q2ω  =  dω(τ/0.001);
monodromy  =  {{q1θ, q2θ}, {q1ω, 
            q2ω}};
monodromy // MatrixForm
Eigensystem[monodromy]
Inverse[Eigenvectors[monodromy]] . monodromy . 
      Eigenvectors[monodromy]

SetPendParam[0.9, 2/3, 0.5, 0]
fourth1  =  Table[{}, {1}];
Runge[0, 1, 0, 10000, 1];
Do[fourth1[[k]]  =  {Mod[third[[k]][[1]] - π, 
            2  π], third[[k]][[2]]}, {k, 1, 
        1}]
fourth2  =  Table[{}, {1}];
Runge[0, 1.0001, 0, 10000, 1];
Do[fourth2[[k]]  =  {Mod[third[[k]][[1]] - π, 
            2  π], third[[k]][[2]]}, {k, 1, 
        1}]
temp1  =  Last[fourth1];
temp2  =  Last[fourth2];
dθτ  =  temp2[[1]] - temp1[[1]];
dωτ  =  temp2[[2]] - temp1[[2]];
q1θ  =  
        dθ(τ/0.0001);
q1ω  =  
        dω(τ/0.0001);
fourth1  =  Table[{}, {1}];
Runge[0, 1, 0, 10000, 1];
Do[fourth1[[k]]  =  {Mod[third[[k]][[1]] - π, 
            2  π], third[[k]][[2]]}, {k, 1, 1}]
fourth2  =  Table[{}, {1}];
Runge[0, 1, 0.0001, 10000, 1];
Do[fourth2[[k]]  =  {Mod[third[[k]][[1]] - π, 
            2  π], third[[k]][[2]]}, {k, 1, 1}]
temp1  =  Last[fourth1];
temp2  =  Last[fourth2];
dθτ  =  temp2[[1]] - temp1[[1]];
dωτ  =  temp2[[2]] - temp1[[2]];
q2θ  =  dθ(τ/0.0001);
q2ω  =  
        dω(τ/0.0001);
monodromy  =  {{q1θ, q2θ}, {q1ω, 
            q2ω}};
Print["The monodromy matrix is  ", 
      monodromy // MatrixForm]
Eigensystem[monodromy];
basis  =  {{Eigenvectors[
                  monodromy][[1]][[1]], Eigenvectors[
                  monodromy][[2]][[1]]}, {Eigenvectors[
                  monodromy][[1]][[2]], Eigenvectors[
                  monodromy][[2]][[2]]}};
Print[basis // MatrixForm];
Print["The diagonalized monodromy matrix is  ", 
      Chop[Inverse[basis] . monodromy . basis] // 
        MatrixForm]
\(Print[\*""Giving Eigenvalues(Lyapunov Exponents) of  (Subscript[λ, 1]) = "", \ \(basis[\([1]\)]\)[\([1]\)], \*"" and (Subscript[λ, 2]) = "", \ \(basis[\([2]\)]\)[\([2]\)]]\)
vector1  =  {Eigenvectors[
                monodromy][[1]][[1]], Eigenvectors[
                monodromy][[2]][[1]]};
vector2  =  {Eigenvectors[
                monodromy][[1]][[2]], Eigenvectors[
                monodromy][[2]][[2]]};
n1  =  (Sqrt[Eigenvectors[monodromy][[1]][[1]]^2 + Eigenvectors[monodromy][[2]][[1]]^2])
n2  =  (Sqrt[Eigenvectors[monodromy][[1]][[2]]^2 + Eigenvectors[monodromy][[2]][[2]]^2])

Monodromy[d_, spacing_, fvalue_, 
        cycles_] := 
Module[{dd  =  d,  
          space  =  spacing, ff  =  fvalue, 
          cyc = cycles}, 
a = 1; aa = 0; c = 0; cc = 1; 
        n11 = 0; n21 = 0; nnew = 0; 
        nnew2 = 0; 
periodrange  =  
          3*π; 
Do[{SetPendParam[ff, 2/3, 0.5, 
              0]; 
fourth1  =  
              Table[{}, {1}]; 
newthing  =  {0, 
                0}; 
Runge[0 + γ*periodrange, a, aa, 
              space, 1]; 
Do[
              fourth1[[k]]  =  {third[[k]][[1]], third[[k]][[2]]}, {k, 1, 1}]; 
fourth2  =  
              Table[{}, {1}]; 
Runge[
              0 + γ*periodrange, a +  dd, aa, space, 
              1]; 
Do[
              fourth2[[k]]  =  {third[[k]][[1]], third[[k]][[2]]}, {k, 1, 1}]; 
temp1  =  
              Last[fourth1]; 
temp2  =  
              Last[fourth2]; 
dθτ  =  
              temp1[[1]] - 
                temp2[[1]]; 
dωτ  =  
              temp1[[2]] - 
                temp2[[2]]; 
q1θ  =  
              dθ(τ/dd); 
q1ω  = 
              dω(τ/dd); 
fourth1  =  
              Table[{}, {1}]; 
Runge[0, c, cc, space, 
              1]; 
Do[
              fourth1[[k]]  =  {third[[k]][[1]], third[[k]][[2]]}, {k, 1, 1}]; 
fourth2  =  
              Table[{}, {1}]; 
Runge[0, c, cc + dd, space, 
              1]; 
Do[
              fourth2[[k]]  =  {third[[k]][[1]], third[[k]][[2]]}, {k, 1, 1}]; 
temp1  =  
              Last[fourth1]; 
temp2  =  
              Last[fourth2]; 
dθτ  =  
              temp1[[1]] - 
                temp2[[1]]; 
dωτ  =  
              temp1[[2]] - 
                temp2[[2]]; 
q2θ  =  
              dθ(τ/dd); 
q2ω  =  
              dω(τ/dd); 
monodromy  =  {{q1θ, q2θ}, {q1ω, 
                  q2ω}}; 
mono  =  
              GramSchmidt[
                Transpose[
                  monodromy]]; 
n1  =  √(Eigenvectors[mono][[1]][[1]]^2 + Eigenvectors[mono][[2]][[1]]^2); 
n2  =  √(Eigenvectors[mono][[1]][[2]]^2 + Eigenvectors[mono][[2]][[2]]^2); 
n11  =  Log[n1]; 
n21  =  Log[n2]; 
nnew = 
              Re[nnew + n11]; 
nnew2 = 
              Re[nnew2 + 
                  n21]; 
a = mono[[1]][[1]]; 
            aa = mono[[1]][[2]]; c = mono[[2]][[1]]; 
            cc = mono[[2]][[2]];}, {γ, 1, 
            cyc}]; 
nnew  =  ((nnew )/(cyc*3  π)); 
nnew2  =  
          (nnew2/(cyc*3  π));
]

lyapunov1 = Table[{0, 0}, {1}];
lyapunov2 = Table[{0, 0}, {1}];
\(\(Do[{
Monodromy[0.0001, 1000, α, 
            250]; 
β = 
            nnew; 
χ = 
            nnew2; 
newthing2\  = \ {{α, β}}; 
newthing3\  = \ {{α, χ}}; 
lyapunov1\  = \ 
            Join[newthing2, lyapunov1]; 
lyapunov2\  = \ 
            Join[newthing3, 
              lyapunov2]; 
Print[\*""The two Lyapunov Exponents are (Subscript[λ, 1]) = "", \ 
            nnew, \ \*"" and (Subscript[λ, 2]) = "", \ 
            nnew2, \ " for a driving force value of ", \ α]}, {α, 1.0, 1.025,  .001}];\)\)
lyapunov1 =  
        Delete[lyapunov1, Length[lyapunov1]];
lyapunov2 =  
        Delete[lyapunov2, Length[lyapunov2]];
MultipleListPlot[lyapunov1, lyapunov2];

lyapunov1 = Table[{0, 0}, {1}];
lyapunov2 = Table[{0, 0}, {1}];
\(\(Do[{
Monodromy[0.0001, 1000, α, 
            500]; 
β = 
            nnew; 
χ = 
            nnew2; 
newthing2\  = \ {{α, β}}; 
newthing3\  = \ {{α, χ}}; 
lyapunov1\  = \ 
            Join[newthing2, lyapunov1]; 
lyapunov2\  = \ 
            Join[newthing3, 
              lyapunov2]; 
Print[\*""The two Lyapunov Exponents are (Subscript[λ, 1]) = "", \ 
            nnew, \ \*"" and (Subscript[λ, 2]) = "", \ 
            nnew2, \ " for a driving force value of ", \ α]}, {α, 1.025, 1.035,  .001}];\)\)
lyapunov1 =  
        Delete[lyapunov1, Length[lyapunov1]];
lyapunov2 =  
        Delete[lyapunov2, Length[lyapunov2]];
MultipleListPlot[lyapunov1, lyapunov2];

lyapunov1 = Table[{0, 0}, {1}];
lyapunov2 = Table[{0, 0}, {1}];
\(\(Do[{
Monodromy[0.0001, 100, α, 
            50]; 
β = 
            nnew; 
χ = 
            nnew2; 
newthing2\  = \ {{α, β}}; 
newthing3\  = \ {{α, χ}}; 
lyapunov1\  = \ 
            Join[newthing2, lyapunov1]; 
lyapunov2\  = \ 
            Join[newthing3, 
              lyapunov2]; 
Print[\*""The two Lyapunov Exponents are (Subscript[λ, 1]) = "", \ 
            nnew, \ \*"" and (Subscript[λ, 2]) = "", \ 
            nnew2, \ " for a driving force value of ", \ α]}, {α, 1.0, 1.025,  .0001}];\)\)
lyapunov1 =  
        Delete[lyapunov1, Length[lyapunov1]];
lyapunov2 =  
        Delete[lyapunov2, Length[lyapunov2]];
MultipleListPlot[lyapunov1, lyapunov2];

GramSchmidt[monodromy]

GramSchmidt[Transpose[monodromy]]

n1
n2

nnew

nnew2

Eigenvalues[monodromy]

Clear[α]

lyapunov1 = Table[{0, 0}, {1}];
lyapunov2 = Table[{0, 0}, {1}];
Do[{
α = 
            Random[Real, {0.8, 1.0}]; 
Clear[
            newthing]; 
Monodromy[0.0001, 
            1000, α]; 
β = 
            newthing[[1]]; 
χ = 
            newthing[[2]]; 
newthing2  =  {{α, β}}; 
newthing3  =  {{α, χ}}; 
lyapunov1  =  
            Join[newthing2, lyapunov1]; 
lyapunov2  =  
            Join[newthing3, lyapunov2];}, {2}];
lyapunov1 =  
        Delete[lyapunov1, Length[lyapunov1]];
lyapunov2 =  
        Delete[lyapunov2, Length[lyapunov2]];
ListPlot[lyapunov1];
ListPlot[lyapunov2];

lyapunov1 =  Delete[lyapunov1, Length[lyapunov1]];

ojhdsf  =  Delete[lyapunov1, Length[lyapunov1]]

ojhdsf

lyapunov1