It is difficult in general to find a solution to this equation. One way would be to use Kronecker calculus, but in this case that would be of no help.
First note that the period of the system is $\pi$ and so the stability condition is given by the existence of positive definite matrix-valued functions $P,Q:[0,2\pi\mapsto\mathbb{R}^{2\times 2}$, $P$ differentiable, $Q$ continuous, verifying $P(0)=P(k\pi)$ and $Q(0)=Q(k\pi)$, for some $k\in\mathbb{N}$, $k\ne 0$, such that
$$
\dot{P}(t)+P(t) A(t) + A(t)^T P(t) + Q(t)=0
$$ holds for all $t\in[0,k\pi]$.
This is not analytically solvable in most cases and one often has to rely on algorithms. However, in the present case, one can exploit the nice structure of the matrix $A(t)$ to parametrize all the possible solutions.
To this aim, define
$$
Z(t)=\begin{bmatrix}
\sin(t) & \cos(t)\\ -\cos(t) & \sin(t)
\end{bmatrix}
$$
we see that $Z(t)^TZ(t)=I$. Define also
$$
J=\begin{bmatrix}
0 & 1\\ -1 & 0
\end{bmatrix}
$$
and observe that $\dot{Z}(t)=Z(t)J^T$.
We now observe that if we use the change of variables $z(t)=Z(t)^Tx(t)$ we get that
$$\dot{z}(t)=\dot{Z}(t)^Tx(t)+Z(t)^T\dot{x}(t)=(\dot{Z}(t)^TZ(t)+Z(t)^TA(t)Z(t))z(t).$$
Using now the fact that $\dot{Z}(t)=Z(t)J^T$, we get that $\dot{Z}(t)^TZ(t)=JZ(t)^TZ(t)=J$ and, so,
$$\dot{z}(t)=J_0z(t)$$
where
$$
J_0=J+Z(t)^TA(t)Z(t)=\begin{bmatrix}
0 & 1\\ -1 & -a
\end{bmatrix}.
$$
So, this shows that by a proper change of variables, one can express the system into LTI form, for which stability is easy to assess. In fact, the eigenvalues of $J_0$ have all negative real part for all $a>0$.
Therefore, there exist some positive definite matrices $P_0$ and $Q_0$ such that
$$J_0^TP_0+P_0J_0+Q_0=0$$
holds.
From that, it is possible to show that $(P_0,Q_0)$ is a solution of the above Lyapunov equation if and only if $(P(t),Q(t))=(Z(t)P_0Z(t)^T,Z(t)Q_0Z(t)^T)$ is a solution to the time-varying Lyapunov equation.
The proof of this result follows from direct substitution of the expressions and very tedious algebra.
This allows us to state that the following statements are equivalent
- The system $\dot{x}(t)=A(t)x(t)$ is asympotically stable.
- The system $\dot{z}(t)=J_0z(t)$ is asympotically stable.
- There exist some positive definite matrices $P_0$ and $Q_0$ such that
$$J_0^TP_0+P_0J_0+Q_0=0$$ holds.
- There exist some positive definite matrix-valued functions $P,Q:[0,2\pi\mapsto\mathbb{R}^{2\times 2}$, $P$ differentiable, $Q$ continuous, verifying $P(0)=P(2\pi)$ and $Q(0)=Q(2\pi)$ such that
$$
\dot{P}(t)+P(t) A(t) + A(t)^T P(t) + Q(t)=0
$$ holds for all $t\in[0,2\pi]$.
Best Answer
The "normal" solution is only for time-invariant systems (as you said), but the idea is the same. So, find the solution at time $t+T$ as $$\begin{align} x(t+T) &= \phi(t+T,t_0) x(t_0) + \int_{t_0}^{t+T} \phi(t+T,\tau) B(\tau) u(\tau) d\tau \\ &= \phi(t+T,t) \phi(t,t_0) x(t_0) + \phi(t+T,t) \int_{t_0}^{t} \phi(t,\tau) B(\tau) u(\tau) d\tau + \int_{t}^{t+T} \phi(t+T,\tau) B(\tau) u(\tau) d\tau \\ &= \phi(t+T,t) x(t) + \int_{t}^{t+T} \phi(t+T,\tau) B(\tau) u(\tau) d\tau \end{align}$$ Assuming Zero-Order Hold and since $B$ is time invariant, we can rewrite these as $$x_{k+1} = A_k x_k + B_k u_k$$ where $$A_k := \phi(t_k+T,t_k) ~~ \text{and} ~~ B_k := \left( \int_{t_k}^{t_k+T} \phi(t_k+T, \tau) d\tau \right) B$$ In general, we cannot go further than that, but in your specific case you can obtain $\phi(\cdot)$ explicitly.