Problem:
There is a set of coupled ODEs, where the time derivatives are equal to quartic polynomials. Here is an example: \begin{equation} \begin{aligned} \frac{dx_1}{dt} &= x_2 x_3 (y_1^2 + y_3^2 + 1) \\ \frac{dx_2}{dt} &= -x_1 x_3 (y_2^2 + y_3^2 + 1) \\ \frac{dx_3}{dt} &= x_1 x_2 (y_1^2 - y_2^2) \\ \end{aligned} \qquad \begin{aligned} \frac{dy_1}{dt} &= y_2 y_3 (x_1^2 + x_3^2 + 1) \\ \frac{dy_2}{dt} &= -y_1 y_3 (x_2^2 + x_3^2 + 1) \\ \frac{dy_3}{dt} &= y_1 y_2 (x_1^2 - x_2^2) \\ \end{aligned} \end{equation}
Periodic (numerical) solutions with period $T$ are wanted, i.e. $\vec{x}(0) = \vec{x}(T)$, $\vec{y}(0) = \vec{y}(T)$. The initial condition is not known, otherwise one would have easily solved it numerically via Rugge-Kutta. However, the boundary condition can be vaguely posed as follows. We start from the linear problem \begin{equation} \begin{aligned} \frac{dx_1}{dt} &= x_2 x_3 \\ \frac{dx_2}{dt} &= -x_1 x_3 \\ \frac{dx_3}{dt} &= 0 \\ \end{aligned} \qquad \begin{aligned} \frac{dy_1}{dt} &= y_2 y_3 \\ \frac{dy_2}{dt} &= -y_1 y_3 \\ \frac{dy_3}{dt} &= \\ \end{aligned} \end{equation} with constant $\vec{x}(t) = (x_1, 0, 0)$ and $\vec{y}(t) = (0, 1, 0)$. Then gradually turn on the non-linear terms. There may be an iterative process leading to the periodic solutions of the non-linear equations.
Analytic Attempt:
For the example problem, $\vec{y} = (0, 1, 0)$ solves the $y$ equations without knowing $x$. Then the $\vec{x}$ equation reduces to \begin{equation} \begin{aligned} \frac{dx_1}{dt} &= x_2 x_3 \\ \frac{dx_2}{dt} &= -2x_1 x_3 \\ \frac{dx_3}{dt} &= -x_1 x_2 \end{aligned} \end{equation}
Possible solutions for $T = \infty$ can be \begin{equation} \begin{aligned} x_1 &= \tanh ( \sqrt{2} t ) \\ x_2 &= \frac{\sqrt{2}}{\cosh( \sqrt{2} t )}\\ x_3 &= - \frac{1}{\cosh( \sqrt{2} t )} \end{aligned} \end{equation}
In this solution $x_1$ goes up from $0$ to $1$, I imagine in a periodic solution, it goes down symmetrically to $0$ on the other end.
Numerically calculation suggests that perturbing the initial conditions can produce periodic solutions, but how the initial condition controls the period is not clear.
Question:
This is probably an analytic solution for the example problem. However, the actual sets of ODEs contains $3n$ variables, so I seek to find a numerical approach.
Is there a variation method to find periodic solutions with unspecified boundary conditions? Suggestions of other methods are also very welcome.
Thank you!
Just a suggestion: there are two conserved quantities, $$x_1^2+x_2^2-x_3^2$$ and similarly for the $y$s. That brings the problem from 6 dimensions down to 4. If you could find 2 more, there would be a lot of theory available to help.