1
$\begingroup$

In this question about Floquet theory the author asked about the fundamental matrix solution $X(t)$ of the following $2\pi$-periodic differential equation $${\displaystyle {\dot {x}}=A(t)x}$$ with

$$A(\text{t})\text{=}\left( \begin{array}{cc} \sin (t) & \sin ^2(t) \cos (t) \\ \sin ^2(t) \cos (t) & \cos (t) \\ \end{array} \right);$$

In the answers, the matrix (monodromy matrix) $X(2\pi) = B$ was calculated using Maple's numerical DE solver: $$ B = \pmatrix{0.736606947094663 & 0.310166738881922\cr -3.21321753950662 &0.00457171990219575\cr}$$ Its eigenvalues are $e^{\pm ic}$ where $c/\pi \approx 0.3791557842561098$ is not rational because it does not have a small denominator: the continued fraction starts $[0; 2, 1, 1, 1, 3, 7, 2, 2, 1, 1, 10, 3]$.

How can we do such computations in Mathematica, especially $B=X(2\pi)$ and $c$?

$\endgroup$

2 Answers 2

2
$\begingroup$

As in the other answer, if you don't need high precision, just use NDSolve.

ClearAll["Global`*"];
A[t_] := {{Sin[t], Sin[t]^2*Cos[t]}, {Sin[t]^2*Cos[t], Cos[t]}}
X0 = {0, 1};
tend = 60;
s = x[t] /. NDSolve[{x'[t] == A[t].x[t], x[0] == X0}, x, {t, 0, tend}];
sol = NDSolveValue[{x'[t] == A[t].x[t], x[0] == X0}, x, {t, 0, tend}];

That you could plot:

Plot[{s[[1, 1]], s[[1, 2]]}, {t, 0, 60}, PlotStyle -> {"Blue", "Red"}]

result

If you are looking for fundamental solution for the given matrix, numerically, you can also use NDSolve:

ClearAll["Global`*"];
A[t_] := {{Sin[t], Sin[t]^2*Cos[t]}, {Sin[t]^2*Cos[t], Cos[t]}}
X0 = {{1, 0}, {0, 1}};
tend = 60;
s = x[t] /. NDSolve[{x'[t] == A[t].x[t], x[0] == X0}, x, {t, 0, tend}];
sol = NDSolveValue[{x'[t] == A[t].x[t], x[0] == X0}, x, {t, 0, tend}];

With the result sol[2 Pi]: {{0.73660731, 0.31016668}, {-3.213216, 0.0045717184}}

$\endgroup$
1
$\begingroup$

One can use

A[t_] := {{Sin[t],Sin[t]^2*Cos[t]},
          {Sin[t]^2*Cos[t],Cos[t]}};

sol[{x0_,y0_},t1_,opts___] := {x[t1],y[t1]} /. NDSolve[Join[
            Thread[{x'[t],y'[t]}==A[t].{x[t],y[t]}],
               {x[0]==x0,y[0]==y0}],{x,y},{t,0,t1},opts][[1]];

X[t_,opts___] := {sol[{1,0},t,opts],
                  sol[{0,1},t,opts]} // Transpose;

X[2*Pi] // MatrixForm

enter image description here


Higher precision. One can use opts. It seems that NDSolve throws errors when increasing precision too much, not sure why exactly. The following settings seem to be ok:

X2Pi = N[X[2*Pi,PrecisionGoal->45,WorkingPrecision->50],22]
(* {{ 0.7366066589412074195214, 0.3101667201208725536489},
    {-3.213215459433480170506,  0.004571639239150046552961}} *)

Det[X2Pi]
(* 1.000000000000000000000, it is known that this must be equal to 1 *)

c = Log[Eigenvalues[X2Pi]]/Pi // Im // First
(* 0.3791557988001068280320 *)

ContinuedFraction[c]
(* {0,2,1,1,1,3,7,2,2,1,1,4,1,1,42,18,1,1,7,2,9,1,2,2} *)
$\endgroup$

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.