Consider a bunch of interacting oscillators (e.g., a chain of atoms), interacting due to anharmonicity in the potential energy. You can Taylor expand the force on each oscillator about equilibrium positions, and get a set of coupled 2nd order ODEs describing the dynamics of this system:

$$ \ddot{q_i}=-\sum_i\omega_i^2q_i-\sum_{ijk}K_{ijk}q_jq_k-\sum_{ijkl}K_{ijkl}q_jq_kq_l-\cdots$$

where indices $i,j,k$ refer to particular modes in the system, $q_i$ are the normal mode coordinates of mode $i$, $K_{ijk}$ are the coupling constants for modes $i,j,k$, and $K_{ijkl}$ are the coupling constants between modes $i,j,k,l$, and so forth. The harmonic part of the force includes a $\omega_i$ term, which is the natural frequency of mode $i$.

I numerically integrated these coupled 2nd order ODEs throughout time and noticed the following peculiarities:

  1. The purely harmonic case is always stable in the sense that the system potential energy fluctuates about equilibrium, and the chain remains intact, provided $\omega_i$ is real. This makes sense according to conventional wisdom on stability and positive 2nd derivatives of the potential energy.
  2. Including cubic terms $K_{ijk}$ causes "instability" in the sense that the system structurally falls apart during numerical time integration. In other words, the chain of oscillators breaks. Why is this?
  3. I imagine there are some combination of anharmonic coupling constants to realize a stable system... Is there a way to find these?
  4. What areas of math/physics study this sort of instability? I don't think the theories of Lyapunov apply here... This sort of instability seems more structural, but I can't find anyone studying the "stability" of coupled nonlinear ODEs like this.
  • 3
    $\begingroup$ Are you sure the instability is physical and not numerical? Have you looked at what the stable time steps need to be as the number of terms increases? What is the stiffness of the equation as you add more terms? $\endgroup$ – tpg2114 Jan 23 '20 at 5:15
  • $\begingroup$ @tpg2114 The instability is physical, because if I make the higher order constants small enough, stability will be achieved. It's not clear to me how these constants affect stability, though. I haven't tested the "stiffness", that's something I just learned about and looks useful. $\endgroup$ – Drew Jan 23 '20 at 5:21
  • 2
    $\begingroup$ Not sure if I agree that small constants -> stability means it is physical. It could be that small constants reduce the stiffness of the system and so the numerical method becomes stable again. There might be some ground work to do first to verify you are using the right numerical method to integrate -- stiff systems need to be treated differently than non-stiff ones, and positive eigenvalues need different methods than negative eigenvalues. $\endgroup$ – tpg2114 Jan 23 '20 at 5:26
  • 1
    $\begingroup$ @tpg2114 Interesting, I'll look into the integration scheme then. The constants $K_i>0$ always, but this isn't enough to ensure stability for some nonzero higher order terms. Right now I'm using Verlet integration - I'll have to look into others. $\endgroup$ – Drew Jan 23 '20 at 5:29
  • 1
    $\begingroup$ You might want to check out my question on our sister site, SciComp.SE about the stability limits of the velocity-Verlet integration method. Might help you get started -- re-reading it, it looks like some values of the spring constant I was using create an unstable scheme, where if the spring constant gets too big relative to the time step, the scheme breaks. $\endgroup$ – tpg2114 Jan 23 '20 at 5:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.