Assume $X(t)$ is a time dependent positive definite symmetric matrix which satisfies the matrix differential equation
$\dot{X}(t) = Q(X(t)), X(0) = 0, $
where $Q(X(t))$ is a symmetric positive definite matrix. To compute the solution $X(t)$, I employed various standard integration methods for explicit ordinary differential equations.
$X(t)$ being positive definite for all time values $t$, we can consider to compute instead $X(t)$, the upper triangular factor $U(t)$ of its Cholesky factorization $X(t) = U(t)^TU(t)$. The Cholesky factor satisfies the implicit matrix differential equation $\dot{U}(t)^TU(t)+U(t)^T\dot{U}(t) = Q(U(t)^TU(t)), U(0) = 0 $
and a smooth solution would be desirable to be computed.
My first question is: Was this problem already considered in the literature?
For the numerical integration of this implicit differential equation I tried to use an implicit method, as the IDA method available in the Sundials package. The main difficulty was to initialize the integrator with a consistent pair $(U(0),\dot{U}(0))$, which does not exist. By slightly perturbing $U(0) = 0 $ to $U(0) = \epsilon I $, the inegration can be sometimes be performed, but the solver reported quite frequently convergence issues of its corrector algorithm.
Using explicit solvers would involve the solution, at each time instance $t$, of the above Lyapunov-like matrix equation satisfied by $\dot{U}(t)$ for a given $U(t)$. Here, enforcing an upper triangular shape of the derivative is one issue which also needs to be addressed. The computation of the derivative may lead to the same difficulties as in the implicit case, at least when evaluating the derivative for $t = 0$.
My second question is: Is any known "trick" to manage the initialization issue?
An a third question: How it would be possible to enforce the smoothness of the solution, because of the inherent non-uniqueness of the Cholesky factorization (e.g., both $U(t)$ and $-U(t)$ are solutions)?