I'm trying to write a program to integrate the motion equations of the pendulum in the damped and forced case, that is, following this equation: $$ \frac{d^2\theta}{dt^2}=-\frac{g}{L}\sin(\theta)-\mu\frac{d\theta}{dt}+F\sin(\Omega t) $$ where $\mu$ is the friction coefficient, $F$ a multiplicative constant for the external force, which in turn is a periodic function with frequency $\Omega$.

I can define the usual system of ODEs without any difficult and use an integrator of the Runge Kutta type, but I encounter some difficulties while trying to use a symplectic integrator, say the Verlet integrators.

The position Verlet I used in different programs is: $$ x^{n+\frac{1}{2}}=x^n+\frac{h}{2}v^n\\ v^{n+1}=v^n+ha(x^{n+\frac{1}{2}})\\ x^{n+1}=x^{n+\frac{1}{2}}+\frac{h}{2}v^{n+1} $$ where $a(...)%$ is the acceleration, defined in the first equation of my question.

Just for the sake of clarity, here is the velocity Verlet: $$ v^{n+\frac{1}{2}}=v^n+\frac{h}{2}a(x^n)\\ x^{n+1}=x^n+hv^{n+\frac{1}{2}}\\ v^{n+1}=v^{n+\frac{1}{2}}+\frac{h}{2}a(x^{n+1}) $$

Here is my problem: in the expression of $a(...)$ the velocity appears in the damping term. But in both the Verlet algorithms I have no clue of how I can pass the velocity as an argument of the function.

In the position Verlet, I suppose I need to know $v^{n+\frac{1}{2}}$, while in the velocity Verlet i can use the initial conditions to kick the algorithm in first place, but in the third equation I think I need to use $v^{n+1}$, that is actually the result of the right hand side.

As I stated, I have no clue of how to resolve this issue, and actually, I'm not even sure if a solution is possible at all. I would be glad for any kind of answer, even a hint would be very kind.

  • 1
    $\begingroup$ This might be off the mark, but your system is non-Hamiltonian so it wouldn't be surprising if this impacted how a symplectic integrator works since they are designed for Hamiltonian systems. The damping term is the thing breaking the Hamiltonian-ness of the system so its maybe not surprising its the one causing you difficulties. $\endgroup$ – jacob1729 Jun 14 '19 at 16:46
  • 1
    $\begingroup$ Essentially a dupe of physics.stackexchange.com/q/239621/25301, as the answer there answers this question, despite the difference of original functions. $\endgroup$ – Kyle Kanos Jun 14 '19 at 19:03

I think that the answer and comments (including references to previous related questions), although correct and also interesting (I have discovered a kind of predictor-corrector variant of the velocity Verlet I was not aware of), are missing a point which could be relevant for your problem.

If the force depends linearly on velocity, the implicit step, which, for example, in the case of velocity Verlet would be $$ v_{n+1}=v_n + \frac{h}{2}\left(~~a_{n+1}(x_{n+1},v_{n+1},t_{n+1})+a_n(x_{n},v_{n},t_n)~~\right), $$ can be easily transformed into an explicit step, by solving a linear equation. Indeed, if $a_n = \tilde a_n(x_n,t_n) + w v_n$, the previous implicit equation for $v_{n+1}$ becomes: $$ v_{n+1}=v_n + \frac{h}{2}\left(~~a_n(x_{n},v_{n})+\tilde a_{n+1}(x_{n+1},t_{n+1}) + w v_{n+1}~~\right) $$ giving $$ v_{n+1} = \frac{v_n+\frac{h}{2}(a_n+\tilde a_{n+1})}{1-\frac{hw}{2} } $$

The resulting algorithm keeps most of the qualities of the usual velocity Verlet.

| cite | improve this answer | |

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.