I am trying to find the limit cycle of a nonlinear oscillator like a duffing or Van-der-Pol oscillator directly using the shooting method. I know how to use the Shooting method to solve a Boundary Value Problem. But I cannot apply the same to an initial value problem like finding the limit cycle of a nonlinear oscillator. I read some papers which discussed the required approach, but I didn't understand them because no examples were discussed. I would have shared them here, but they are not open access.
The whole point of using the Shooting method here is to get the limit cycle efficiently without much computation compared to integrating over the entire time interval.
I tried asking this question in Physics Stack Exchange, and I was advised to post it here as it is related to Numerical Methods. I want to get the limit cycle of a nonlinear oscillator in cases where it exists.
Best Answer
One of possible approaches. Assuming you have an autonomous system $$ \dot {\mathbf y} = \mathbf F(\mathbf y) $$ finding a limit cycle is finding a fixed point of its Poincare map (the wiki article is overcomplicated, the only thing you need there is the picture).
To begin with we need a Poincare section — that is a surface that definitely crosses the limit cycle. For Van-der-Pol oscillator that might be $y_1 = 0$ or $y_2 = 0$, each is a valid one. Let's use $y_1 = 0$.
The limit cycle we're searching is a solution to the following problem $$ \dot {\mathbf y} = \mathbf F(\mathbf y)\\ \mathbf y(0) = (0, z)\\ \mathbf y(T) = (0, z) $$ where unknowns are $z$ ($y_2$ coordinate of the intersection of the cycle and the Poincare section) and $T$ (cycle period).
Let's describe one shot procedure. Assume that we have some approximation $z_1$ to the solution $z$. Let's integrate the following ODE problem $$ \dot {\mathbf y} = \mathbf F(\mathbf y)\\ \mathbf y(0) = (0, z_1) $$ using some ODE solver. There is no fixed finish time, we are integrating until the following condition holds: $$ y_{n,1} < 0\\ y_{n+1,1} \geq 0 $$ Here $\mathbf y_n$ denotes numerical solution at $t_n$, and ${\mathbf y}_{n+1}$ is the solution at $t_{n+1}$. And $y_{n,1}$ denotes the first component of the $\mathbf y_n$ vector.
The condition tells that somewhere in $t \in (t_n, t_{n+1}]$ the numerical solution intersects the Poincare section $y_1 = 0$. Let $(t^*, \mathbf y^*)$ be the intersection time and point. Note that $\mathbf y^* = (0, z_2)$.
Let's do one more integration step and compute $\mathbf y_{n+2}$. Now we have the following table $$ \begin{array}{c|c||c||c|c} t_{n-1}&t_n&t^*&t_{n+1}&t_{n+2}\\\hline y_{n-1,1}&y_{n,1}&0&y_{n+1,1}&y_{n+2,1}\\\hline y_{n-1,2}&y_{n,2}&z_2&y_{n+1}&y_{n+2,2} \end{array} $$ The $t^*$ and $z_2$ values can be computed from the condition $y^*_1 = 0$ using the inverse interpolation.
The shot procedure describes a function $$ z_2 = G(z_1), $$ and the cycle corresponds to its fixed point $z = G(z)$. This equation can be solved numerically with many different methods, e.g. bisection, secant and many others. Since we are looking for a limiting cycle $z^{(n+1)} = G(z^{(n)})$ iterations will be convergent provided some sane $z^{(0)}$ is given.
Sample implementation in python