Fundamental matrix solution of a differential equation $x'=A(t)x$

Question

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$?

Answer 1

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"}]

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}}

Answer 2

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

---

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} *)