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

In this question about Floquet theory the author asked about the fundamental matrix solution $$X(t)$$ of the following $$2\pi$$-periodic differential equation $${\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$$?

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 == X0}, x, {t, 0, tend}];
sol = NDSolveValue[{x'[t] == A[t].x[t], x == 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 == X0}, x, {t, 0, tend}];
sol = NDSolveValue[{x'[t] == A[t].x[t], x == X0}, x, {t, 0, tend}];


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

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[
{x==x0,y==y0}],{x,y},{t,0,t1},opts][];

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