Problem statement
In a problem of interacting particles, I encountered a type of geodesic equation in $\mathbb{R}^n$ with an additional rotation and dilation term $$ \ddot\gamma(t) + e^{t Q} \Lambda e^{-t Q} \gamma(t) = 0 , \tag{1} $$ where $Q=-Q^T$ is a given skew-symmetric matrix and $\Lambda = \Lambda^T \geq 0$ is a symmetric non-negative definite matrix. Both $Q$ and $\Lambda$ are fixed.
In addition, the initial and terminal position is fixed, i.e. $$ \gamma(0) = x \qquad\text{and}\qquad \gamma(1) = y . $$ In particular, I'm interested in the corresponding Neumann values, i.e. $\dot \gamma(0)$ and $\dot \gamma(1)$, (depending on the parameters $x,y, Q, \Lambda$). More specifically, the further investigation relies on the product-difference $$ \gamma(1)\cdot \dot\gamma(1) - \gamma(0)\cdot \dot \gamma(0) . $$
Solution approach so far
The case $Q=0$ or $\Lambda = \lambda \operatorname{Id}$ is easy solvable given in terms of elementary trigonometric functions.
The problem can be made time-homogeneous by the change of variable $\xi = e^{-t Q} \gamma$ resulting in the system $$ \ddot \xi + 2 Q \dot\xi + (\Lambda + Q^2) \xi = 0 . $$ The fundamental solution of this system can be written down (at least for $n=2$). However, after resolving for the boundary values, the calculations and parameter-dependencies become very untraceable, even in two dimensions.
Numerical tests
The ODE $(1)$ can be solved easily numerically. I include some pictures below for $n=2$ to explain the effect of $Q$ and $\Lambda$. Parameter is always $x=(0,0), y=(1,1)$, $Q=\begin{pmatrix} 0 & -q \\ q & 0 \end{pmatrix}$, $\Lambda=\begin{pmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{pmatrix}$
$q=\lambda_1=\lambda_2=0$ leads to Euclidean geodesic $q=\lambda_1=\lambda_2=0$ leads to Euclidean geodesic" />
$q=\lambda_2=0$, $\lambda_1=5$ favors excursion in $x_1$ direction $q=\lambda_2=0$, $\lambda_1=5$ favors excursion in $x_1$ direction" />
$q=0$, $\lambda_1=5$, $\lambda_2=6$ favors excursion in both $x_1, x_2$ direction $q=0$, $\lambda_1=5$, $\lambda_2=6$ favors excursion in both $x_1, x_2$ direction" />
$q=1$, $\lambda_1=5$, $\lambda_2=0$ gives a kind of spin introducing some wobbling effect $q=1$, $\lambda_1=5$, $\lambda_2=0$ gives a kind of spin introducing some wobbling effect" />
Mathematica
can produce the images via the commands and the parameter dependency can be easily explored
$Assumptions = l1 >= 0 && l2 >= 0 && q \[Element] Reals
Rot[q_, t_] = MatrixExp[{{0, -q}, {q, 0}} t]
sol = ParametricNDSolveValue[{{x1''[t], x2''[t]} +
Rot[q, t] .{{l1, 0}, {0, l2}}.Rot[-q, t].{x1[t], x2[t]} == {0,
0}, x1[0] == X1, x2[0] == X2, x1[1] == Y1, x2[1] == Y2}, {x1,
x2}, {t, 0, 1}, {l1, l2, q, X1, X2, Y1, Y2}]
Manipulate[
ParametricPlot[
Evaluate[#[t] & /@ sol[l1, l2, q, X1, X2, Y1, Y2]], {t, 0,
1}], {{l1, 0}, 0, 10}, {{l2, 0}, 0, 10}, {{q, 0}, -2 \[Pi],
2 \[Pi]}, {{X1, 0}, -3, 3}, {{X2, 0}, -3, 3}, {{Y1, 1}, -3,
3}, {{Y2, 1}, -3, 3}]
Question
I'm happy about any input on this problem or references, where similar issues are studied. Mainly, I wonder if the Dirichlet-to-Neumann map for second-order ODEs is somewhere investigates in the literature.