I Symplectic integrator, non-separable Hamiltonian

654
148
I have been attempting to modify a symplectic integrator that I wrote a while ago. It works very well for "separable" hamiltonians, but I want to use it to simulate a double pendulum.

I am using the Stormer-Verlet equation (3) from this source. From the article "Even order 2 follows from its symmetry."

In the case of a separable hamiltonian, ##\nabla q## is a function only of ##q##, and similarly for ##p##, so that the equations form a symmetrical sequence of three function calls.

For a non-separable hamiltonian, this is no longer true, and it is necessary to use the full equations, but these are no longer symmetrical (the first two are implicit whilst the last is explicit).

Anyhow, I have implemented the full non symmetrical equations, and while they are "symplectic" in the sense that there is no systematic "energy creep" in the output, they are only first order WRT step size, and my attempts to increase the order via composition are ineffective (the composition is still first order).

So, my question is this: is it possible to compose these implicit equations, or does their asymmetry prevent this? In other words, have I just made an error somewhere in my implementation?
 
1,416
9
Thanks for the thread! This is an automated courtesy bump. Sorry you aren't generating responses at the moment. Do you have any further information, come to any new conclusions or is it possible to reword the post? The more details the better.
 
654
148
Thanks for the thread! This is an automated courtesy bump. Sorry you aren't generating responses at the moment. Do you have any further information, come to any new conclusions or is it possible to reword the post? The more details the better.
np ;) I've checked this over, stepped away for a few days, and still can't find any mistakes in my implementation.
This is an invocation that I use a lot!
Python:
time -p ./DoublePendulum.py <initial-conditions.double-pendulum.json >/tmp/data-Python
Here is a data file (rename to initial-conditions.double-pendulum.json):
JavaScript:
{
  "Simulator": "Double Pendulum",
  "IC": {
    "g": 1.0,
    "m": 1.0,
    "l1": 1.0,
    "l2": 1.0,
    "m1": 1.0,
    "m2": 1.0,
    "th1": 0.4,
    "th2": 0.8,
    "pth1": 0.0,
    "pth2": 0.0,
    "start": 0.0,
    "end": 100.0,
    "step": 0.01,
    "integrator": "b4",
    "scheme" : "suzuki",
    "tol" : 1.0e-15,
    "plotratio": 1
  }
}
and here is the Python code (rename and chmod +x DoublePendulum.py):
Python:
#!/usr/bin/env python3

from sys import stdin, stderr, argv
from math import sin, cos, log10
from json import loads


class Dual:

    def __init__(self, jet):
        self.val = jet[0]
        self.der = jet[1]

    def __pos__(self):
        return Dual([self.val, self.der])

    def __neg__(self):
        return Dual([- self.val, - self.der])

    def __add__(self, other):
        if isinstance(other, Dual):
            return Dual([self.val + other.val, self.der + other.der])
        else:
            return Dual([self.val + other, self.der])

    def __radd__(self, other):
        return self.__add__(other)

    def __sub__(self, other):
        if isinstance(other, Dual):
            return Dual([self.val - other.val, self.der - other.der])
        else:
            return Dual([self.val - other, self.der])

    def __rsub__(self, other):
        return Dual([- self.val + other, - self.der])

    def __mul__(self, other):
        if isinstance(other, Dual):
            return Dual([self.val * other.val, self.der * other.val + self.val * other.der])
        else:
            return Dual([self.val * other, self.der * other])

    def __rmul__(self, other):
        return self.__mul__(other)

    def __truediv__(self, other):
        if isinstance(other, Dual):
            return Dual([self.val / other.val, (self.der * other.val - self.val * other.der) / other.val**2])
        else:
            return Dual([self.val / other, self.der / other])

    def __rtruediv__(self, other):
        return Dual([other / self.val, - other * self.der / self.val**2])

    def __pow__(self, a):
        return Dual([self.val**a, a * self.val**(a - 1.0) * self.der])

    @property
    def var(self):
        return Dual([self.val, 1.0])

    @property
    def sqr(self):
        return Dual([self.val**2, 2 * self.der * self.val])

    @property
    def sin(self):
        return Dual([sin(self.val), self.der * cos(self.val)])

    @property
    def cos(self):
        return Dual([cos(self.val), - self.der * sin(self.val)])


class DoublePendulum(object):
    def __init__(self, l1, m1, l2, m2, th1_0, pth1_0, th2_0, pth2_0, tol):
        self.l1 = l1
        self.m1 = m1
        self.l2 = l2
        self.m2 = m2
        self.th1 = Dual([th1_0, 0.0])
        self.pth1 = Dual([pth1_0, 0.0])
        self.th2 = Dual([th2_0, 0.0])
        self.pth2 = Dual([pth2_0, 0.0])
        self.h0 = self.hamiltonian(self.th1, self.pth1, self.th2, self.pth2).val
        self.tol = tol

    def hamiltonian(self, th1, pth1, th2, pth2):
        return (self.l2**2 * self.m2 * pth1.sqr + self.l1**2 * (self.m1 + self.m2) * pth2.sqr
                - 2.0 * self.m2 * self.l1 * self.l2 * pth1 * pth2 * (th1 - th2).cos) \
               / (2.0 * self.l1**2 * self.l2**2 * self.m2 * (self.m1 + self.m2 * (th1 - th2).sin.sqr)) \
               - (self.m1 + self.m2) * self.l1 * th1.cos - self.m2 * self.l2 * th2.cos

    # def hamiltonian_single_pendulum(self, th1, pth1, th2, pth2):
    #     return 0.5 * pth1.sqr / (self.m1 * self.l1**2) - self.m1 * self.l1 * th1.cos

    def qth1_update(self, c, th1):
        return self.th1.val - th1 + 0.5 * c * self.hamiltonian(Dual([th1, 0.0]), self.pth1.var, self.th2, self.pth2).der

    def qth2_update(self, c, th2):
        return self.th2.val - th2 + 0.5 * c * self.hamiltonian(self.th1, self.pth1, Dual([th2, 0.0]), self.pth2.var).der

    def q_update_1(self, c):
        th1 = secant(self.qth1_update, self.th1.val, c, self.tol)
        th2 = secant(self.qth2_update, self.th2.val, c, self.tol)
        self.th1 = Dual([th1, 0.0])
        self.th2 = Dual([th2, 0.0])

    def q_update_2(self, c):
        th1 = self.th1.val + 0.5 * c * self.hamiltonian(self.th1, self.pth1.var, self.th2, self.pth2).der
        th2 = self.th2.val + 0.5 * c * self.hamiltonian(self.th1, self.pth1, self.th2, self.pth2.var).der
        self.th1 = Dual([th1, 0.0])
        self.th2 = Dual([th2, 0.0])

    def pth1_update(self, d, pth1):
        th1_var = self.th1.var
        return self.pth1.val - pth1 - 0.5 * d * (self.hamiltonian(th1_var, self.pth1, self.th2, self.pth2).der
                                               + self.hamiltonian(th1_var, Dual([pth1, 0.0]), self.th2, self.pth2).der)

    def pth2_update(self, d, pth2):
        th2_var = self.th2.var
        return self.pth2.val - pth2 - 0.5 * d * (self.hamiltonian(self.th1, self.pth1, th2_var, self.pth2).der
                                               + self.hamiltonian(self.th1, self.pth1, th2_var, Dual([pth2, 0.0])).der)

    def p_update(self, d):
        pth1 = secant(self.pth1_update, self.pth1.val, d, self.tol)
        pth2 = secant(self.pth2_update, self.pth2.val, d, self.tol)
        self.pth1 = Dual([pth1, 0.0])
        self.pth2 = Dual([pth2, 0.0])

    def stormer_verlet(self, h):
        self.q_update_1(h)
        self.p_update(h)
        self.q_update_2(h)

    def solve(self, h, start, end, tr):
        t = 0.0
        i = 0
        while t < end:
            if t >= start and i % tr == 0:
                self.plot(t)
            self.stormer_verlet(h)
            i += 1
            t = h * i
        self.plot(t)

    def plot(self, time):
        x1 = self.l1 * sin(self.th1.val)
        y1 = - self.l1 * cos(self.th1.val)
        x2 = x1 + self.l2 * sin(self.th2.val)
        y2 = y1 - self.l2 * cos(self.th2.val)
        error = abs(self.hamiltonian(self.th1, self.pth1, self.th2, self.pth2).val - self.h0)
        print("{:.9e} {:.9e} {:.9e} {:.9e} {:.5e} {:.9e}".format(
            x1, y1, x2, y2, time, 10 * log10(error if error > 1.0e-18 else 1.0e-18)))


def secant(func, parameter, cd, tol, spread=0.001):
    a = parameter * (1 - spread) + (-1 if parameter * -1 > 0 else 1) * spread
    b = parameter * (1 + spread) + (1 if parameter * -1 < 0 else -1) * spread
    f_a = func(cd, a)
    f_b = func(cd, b)
    counter = c = f_c = 1
    while abs(f_c) > tol:
        c = (b * f_a - a * f_b) / (f_a - f_b)
        f_c = func(cd, c)
        b = a
        f_b = f_a
        a = c
        f_a = f_c
        counter += 1
        if counter == 1000:
            raise RuntimeError("Giving Up!")
    return c


if __name__ == "__main__":
    print("Simulator: {}".format(argv[0]), file=stderr)
    input_data = stdin.read()
    ic = loads(input_data)['IC']
    print(input_data, file=stderr)
    DoublePendulum(ic['l1'], ic['m1'], ic['l2'], ic['m2'], ic['th1'], ic['pth1'], ic['th2'], ic['pth2'], ic['tol']).solve(
        ic['step'], ic['start'], ic['end'], ic['plotratio'])
else:
    print(__name__ + " module loaded", file=stderr)
 
654
148
OK, no response from stack exchange.

I have backed away for a while and returned to experimenting with this, and my conclusions remain the same: non-separable (implicit) symplectic integrators (being non-symmetrical) are not composable i.e. the Stormer-Verlet method is only first-order for non-separable Hamiltonians.

The examples I have been using are the double pendulum, and the Kerr metric (MTW p.900 equation 33.35).

It would be really nice to get confirmation (or refutation!) of this - surely someone on Physics Forums knows the answer?
 

Want to reply to this thread?

"Symplectic integrator, non-separable Hamiltonian" You must log in or register to reply here.

Related Threads for: Symplectic integrator, non-separable Hamiltonian

Replies
3
Views
3K
Replies
4
Views
6K
  • Posted
Replies
0
Views
1K
Replies
4
Views
1K
  • Posted
Replies
4
Views
5K
  • Posted
Replies
0
Views
797
  • Posted
Replies
0
Views
3K
Replies
0
Views
3K

Physics Forums Values

We Value Quality
• Topics based on mainstream science
• Proper English grammar and spelling
We Value Civility
• Positive and compassionate attitudes
• Patience while debating
We Value Productivity
• Disciplined to remain on-topic
• Recognition of own weaknesses
• Solo and co-op problem solving

Hot Threads

Top