What I have in mind is a mechanical system that is described by an implicit system of ODEs or a system of DAEs (*differential algebraic equations*). The system is asymptotically stable, meaning that there exists an attracting orbit that the system converges to for all initial values (except possibly for the trivial initial condition).


I am searching for a way to numerically calculate such a system's period when it has converged arbitrarly close to the attracting orbit. I have tried finding litterature on this subject, but find none that discuss how to solve this problem.


There is one pragmatic way to approach this problem, and that is of course to calculate the period by simulating the system (e.g. in Matlab) for a large amount of cycles and then divide the total time with the cycles: $T=\frac{t\_{total}}{N\_{cycles}}$, where $T$ is the period (or eigenfrequency) of the vibrating system.


Alas, there is a trade-off in precision doing so. Anyone know a more rigorous approach using some kind of numerical method?


Best regards



Nonsmooth dynamical system (DAE) - systematic way to calculate period numerically?
I would program this using an adaptive step size approach where the step size decreases as some power (square root works well, empirically) of the error associated with an assumed period. Specifically:


Start with a point $y\_0$ in the dynamics space that is "near" the asymptotic periodic orbit, some time step size $\Delta\_0$ which need not be as small as would be required for your eventual desired precisions, and some estimate $\tau\_0$ of the period. Divide that period by the step size to find $n\_1 = \left\lceil \frac{\tau\_0}{\Delta\_0} \right\rceil$, and simulate using $n\_1$ steps of size $\frac{\tau\_0}{n\_1}$. 


At this point the simulated trajectory will take you to $y\_1$ and let the previous state for one less and subsequent state for one more time step be $x\_1$ and $z\_1$, respectively. Now let $y\_1$ be the point of closest approach of the line $x\_1z\_1$ to the point $x\_0$, and let that point be at $\tilde{y}\_1 = \xi\_1 x\_1 + (1-\xi\_1 z\_1)$. Then your next guess for the period will be 
$$t\_1 = t\_0 + (2\xi\_1 - 1) \Delta\_0$$
and your degree of regard for the solution is $r\_0 = |y\_1-y\_0|$ for some suitably chosen norm function. 


Repeat this step, but every $k$-th iteration (for some suitably chosen $k$, balancing the cost of changing to a different time step size in the dynamic simulation against the desire to adjust the step size aggressively) modify your step size, doing
$$
\Delta\_{nk} = \Delta\_{(n-1)k} \sqrt{\frac{r\_{nk}}{r\_{(n-1)k}}}
$$
As I said, the square root might not be the best choice for how aggressively to decrease the step size; for well behaved systems one might do some theoretical work on what the best exponent should be, but $frac12$ seems like a good heuristic choice.


This scheme has the advantage that you can allow the trajectory to start some small distance away from the periodic orbit, regarding the imperfection in early periodicity to be comparable to the errors already inherent in the iterative approximations to a perfect time period. It also, once you have gotten close to the actual orbit and period, converges geometrically, which is much better than just going around for multiple periods.


If the Matlab routine does not let you recover the state for a step before and a step after the specified time interval, you may have to stop one interval early and move forward the last two intervals by explicit linear approximation to the equations.


