numerical investigation of stability of motion (confinement)

I am trying to find the required specifications of a RF trap, in which a proton can be confined.(trap dimensions,voltage frequency and amplitude used, etc). I have to solve the equations of motion numerically because the potential doesn't have a closed form for this specific geometry and equations of motion can not be solved analytically (indeed , only a series solution can be obtained for potential using Legendre polynomials and this series is present in the system of equations of motion) .This is not the case for hyperbolic Paul trap where the equations are happen to be Mathieu equations.

for this purpose, I solve this system of ODEs numerically (ode45--> 4th and 5th order Runge-Kutta) and see if the ion leaves the trap. you know, I can simulate the motion for a limited time span (typically about some milliseconds at most). So if the ion happen to leave the trap after a longer time span , I can't notice that.

The question is that how and when I can claim that the ion is confined in the trap?(using numerical methods)

the equations are:(in 2D polar coordinates with variables $r$ and $\theta$)(put some integer instead of $\infty$ in the series!)($C$ , $U$ , $V$ , $R$ and $\alpha$ are constants)

$$\ddot{r}-r\dot{\theta}^{2}=C(U+Vcos({\omega}t))\sum_{n=0}^{\infty}n(4n+1)\biggr[\sum_{m=0}^{n} (-1)^{m}\frac{(4n-2m)!(1-(cos(\alpha))^{2n-2m+1})}{4^{n}m!(2n-m)!(2n-2m+1)!}\biggr]\frac{r^{2n-1}}{R^{2n}}P_{2n}(cos(\theta))$$

$$r\ddot{\theta}+2\dot{r}\dot{\theta}=C(U+Vcos({\omega}t))\sum_{n=0}^{\infty}(4n+1)\biggr[\sum_{m=0}^{n} (-1)^{m}\frac{(4n-2m)!(1-(cos(\alpha))^{2n-2m+1})}{4^{n}m!(2n-m)!(2n-2m+1)!}\biggr]\frac{r^{2n-1}}{R^{2n}}{d \over d\theta}(P_{2n}(cos(\theta)))$$

which with this assignment:

$$X_1=r$$

$$X_2=\theta$$

$$X_3=\dot{r}$$

$$X_4=\dot{\theta}$$

becomes a simultaneous system of ODEs:

$\dot{X_1}=X_3$

$\dot{X_2}=X_4$

$\dot{X_3}=X_1X_4^{2}+C(U+Vcos({\omega}t))\sum_{n=0}^{\infty}n(4n+1)\biggr[\sum_{m=0}^{n} (-1)^{m}\frac{(4n-2m)!(1-(cos(\alpha))^{2n-2m+1})}{4^{n}m!(2n-m)!(2n-2m+1)!}\biggr]\frac{r^{2n-1}}{R^{2n}}P_{2n}(cos(X_2))$

$\dot{X_4}=\frac{-2X_3X_4}{X_1}+\frac{C(U+Vcos({\omega}t))\sum_{n=0}^{\infty}(4n+1)\biggr[\sum_{m=0}^{n} (-1)^{m} \frac{(4n-2m)!(1-(cos(\alpha))^{2n-2m+1})}{4^{n}m!(2n-m)!(2n-2m+1)!}\biggr]\frac{r^{2n-1}}{R^{2n}}{d \over dX_2}(P_{2n}(cos(X_2)))}{X_1}$

-
What about the problem keeps you from simulating for a long time span? Stability? Computational cost? –  Godric Seer Feb 13 '13 at 14:13
Hi @Peasant, and welcome to Scicomp! Could you provide the equations of motion that you are solving? –  Paul Feb 13 '13 at 15:53
@Godric: for some parameters(trap specifications and particle's initial conditions) the solution -for r variable (in 2D polar coordinates)- goes to infinity almost immediately (after about some nanoseconds) .But in general, I can't continue simulation for more than (at most)some milliseconds because of the increasing time needed for computations.(I have to solve this system for many different parameters , so time is important) –  Mostafa Feb 13 '13 at 23:06
@Paul:equations are too big so I can't write them here.(a simultaneous system of ODEs; in 2D polar coordinates, $r$, $\theta$) –  Mostafa Feb 13 '13 at 23:12

In your case, you are designing a system in which the ion will undergo a periodic motion. It is true that there are systems in which even autonomous ODEs (i.e., ODEs in which the forcing term is not time dependent -- i.e., like in your case where the field is fixed) can give rise to chaotic motion, but I would be very surprised to see that in your system. Rather, you are almost certainly going to see closed orbits of the ion.

For such systems, while it's not a guarantee that the orbit is stable, if you can observe a few hundred or thousand orbits and not see the ion leave the trap, then I imagine that that is good enough for practical purposes.

-
Thanks a lot. Do you mean an exactly periodic motion, or a motion around the center (in which the particle may never have the same status for two different $t$s)?(I have seen in some cases the $r$ variable has an almost periodic motion,but I've never seen an exactly periodic motion in the solutions). and, another question ; how many terms do you think I should keep in the series present in equations(written below in an answer (couldn't be written here)) to have a tolerable error(just for being able to distinguish trapped ones from non trapped ones)? Thanks. –  Mostafa Feb 13 '13 at 23:30
The motion will only be exactly periodic for very special potentials. Yours may simply not be that kind in which case you get a rotating ellipse like motion. In that case, $r$ is periodic, $\phi$ is too, but $\mathbf x$ is not. –  Wolfgang Bangerth Feb 14 '13 at 0:36
For the number of terms -- I don't know. –  Wolfgang Bangerth Feb 14 '13 at 0:36

If you do a finite simulation you can use two stopping criteria:

• It leaves the trap (It is obviously not trapped)
• The status of the system repeats itself (It is trapped)

The second condition requires that randomness does not play a significant role in your simulation.

With sufficient knowledge of the system you can add more stopping criteria:

• The status of the system is in a range of predefined 'safe' states
• The status of the system is 'between' previous states (does that mean it is safe?)

Of course you can add as many as you like, and may want to interpret the result as follows: a particle not leaving during a certain period, and also not meeting the stopping criteria is 'probably' trapped.

-