# More cool automatic differentiation (Hamiltonian)

• Python
Well I think it is cool anyhow ;)

Here is a dependency-free variable-order double pendulum simulation in about 120 lines of Python. Have you seen the equations of motion for this system?

As usual, this is based on code that I have provided here, but trimmed for maximum impact. Can you see how it works?

Code:
./double.py 8 0.1 1000 1 1 1 1 .1 0 .1 0
Python:
#!/usr/bin/env python3

from sys import argv
from math import sin, cos

class Dual:

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

def __str__(self):
return "{:.6e} {:.6e}".format(self.val, self.der)

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

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

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

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) * self.der])

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

@property
def sqr(self):
return self * self

@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)])

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

def main():
n = int(argv[1])
h = float(argv[2])
steps = int(argv[3])

g = 1.0
l1, m1, l2, m2 = float(argv[4]), float(argv[5]), float(argv[6]), float(argv[7])
th1_0, pth1_0, th2_0, pth2_0 = float(argv[8]), float(argv[9]), float(argv[10]), float(argv[11])
th1 = [1.0 for _ in range(n + 1)]
pth1 = [1.0 for _ in range(n + 1)]
th2 = [1.0 for _ in range(n + 1)]
pth2 = [1.0 for _ in range(n + 1)]
for step in range(1, steps + 1):
th1[0], pth1[0], th2[0], pth2[0] = th1_0, pth1_0, th2_0, pth2_0
th1_dual = Dual([th1[0], 0.0])
pth1_dual = Dual([pth1[0], 0.0])
th2_dual = Dual([th2[0], 0.0])
pth2_dual = Dual([pth2[0], 0.0])
for k in range(n):
th1[k + 1] = hamiltonian(g, l1, m1, l2, m2, th1_dual, pth1_dual.var, th2_dual, pth2_dual).der / (k + 1)
pth1[k + 1] = - hamiltonian(g, l1, m1, l2, m2, th1_dual.var, pth1_dual, th2_dual, pth2_dual).der / (k + 1)
th2[k + 1] = hamiltonian(g, l1, m1, l2, m2, th1_dual, pth1_dual, th2_dual, pth2_dual.var).der / (k + 1)
pth2[k + 1] = - hamiltonian(g, l1, m1, l2, m2, th1_dual, pth1_dual, th2_dual.var, pth2_dual).der / (k + 1)
result_th1, result_pth1, result_th2, result_pth2 = th1[n], pth1[n], th2[n], pth2[n]
for i in range(n - 1, -1, -1):
result_th1 = result_th1 * h + th1[i]
result_pth1 = result_pth1 * h + pth1[i]
result_th2 = result_th2 * h + th2[i]
result_pth2 = result_pth2 * h + pth2[i]
th1_0, pth1_0, th2_0, pth2_0 = result_th1, result_pth1, result_th2, result_pth2
x1 = l1 * sin(th1_0)
y1 = - l1 * cos(th1_0)
x2 = x1 + l2 * sin(th2_0)
y2 = y1 + - l2 * cos(th2_0)
print("{:.9e} {:.9e} {:.9e} {:.9e} {:.5e}".format(x1, y1, x2, y2, step * h))

main()

Last edited:
scottdave

## Answers and Replies

Related Programming and Computer Science News on Phys.org
jedishrfu
Mentor
Thanks for sharing!

One thing I noticed is the call to main() versus a better a python convention as follows:

Python:
if __name__ == '__main__':
main()
This allows you to run your program directly or use it in the context of another program ala library of functions.

https://stackoverflow.com/questions/419163/what-does-if-name-main-do

An interesting aside here is that some folks at Cornell did a learning system that discovers the equations of motion for this kind of system:

https://www.wired.com/2009/04/newtonai/

Thanks for sharing!
Cheers, I'm trying to drum up a bit of awareness of these techniques, they are not talked about much on the internet! This snippet uses automatic differentiation in two different ways, and I have never seen an example of this approach before, it was just a little experiment.
One thing I noticed is the call to main() versus a better a python convention as follows:

Python:
if __name__ == '__main__':
main()
This allows you to run your program directly or use it in the context of another program ala library of functions.

https://stackoverflow.com/questions/419163/what-does-if-name-main-do
Yep, this is a cut'n'paste from "real" code. I usually start off with this:
Python:
if __name__ == "__main__":
main()
else:
print(__name__ + " module loaded", file=stderr)
and cut bits out if I decide that "dual use" is not appropriate.
An interesting aside here is that some folks at Cornell did a learning system that discovers the equations of motion for this kind of system:

https://www.wired.com/2009/04/newtonai/
A bit scant on details, reminds me a bit of the Nelder-Mead approach (or maybe Monte-Carlo) ;) I wonder if they have made any progress since then?

jedishrfu
Mentor
My feeling was that they used a genetic algorithm with each vector in the set of solutions containing an expression to mutate and to test.

There's also a website http://www.compadre.org/osp Open Source Physics that implement many of these kinds of things in Java. They provide an ODE solver (Euler, RungeKutta...) and you get to code the details of how it steps.

Its great though that you've provided one in Python.

jedishrfu
Mentor
Have you considered writing an Insight article on this? @Greg Bernhardt might be interested.

Greg Bernhardt
I have considered it, and haven't written off the idea, but you should know that I am more than 30 years out of an academic environment, and I didn't do Maths or Physics (Electrical & Electronic Engineering)!

I have some background material in the form of a tex file in the source, which I have used to capture the idea, but if you read it you will see that I am not great at mathematical rigour ;) Are there any mathematicians out there that could proof-read and correct my mistakes?

Also I need to add something about the hamiltonian solver I mentioned in the OP as I only just got it going (I wasn't even sure until today that it was doable).

Last edited:
Full gmpy2 double pendulum with plotter now on GitHub (see README).

There's also a website http://www.compadre.org/osp Open Source Physics that implement many of these kinds of things in Java. They provide an ODE solver (Euler, RungeKutta...) and you get to code the details of how it steps.

Its great though that you've provided one in Python.
I've skimmed the docs and I'm not convinced they are using the same methods. They might be using the Taylor Series Method to solve ODEs (the docs I found were not detailed enough), but so far I'm pretty sure they aren't using dual numbers to differentiate hamiltonians. I'll have to find the sources to be sure . . .

I'll report back when I find out more. [UPDATE] I really don't like their site, I haven't managed to find any sources so far!

Last edited:
jedishrfu
Mentor
The OSP site java code is freely available. In addition, there are numerous examples. I've used it in some projects. Keep in mind that its primarily a teaching tool, not a high-performance tool and so they may not implement the most efficient methods.

Here's a writeup on their uses as a teaching tool:

I can't say whether they've implemented the Taylor Series method or not. However, if its a well-known method used in numerical methods then I'm sure someone somewhere has implemented it before.

Don't let that stop you though. It appears you have some worthy code that people should know about and a PF Insight article would help promote its use.

The OSP site java code is freely available. In addition, there are numerous examples. I've used it in some projects. Keep in mind that its primarily a teaching tool, not a high-performance tool and so they may not implement the most efficient methods.

Here's a writeup on their uses as a teaching tool:

I can't say whether they've implemented the Taylor Series method or not. However, if its a well-known method used in numerical methods then I'm sure someone somewhere has implemented it before.
There are already many implementations in many languages, going back half a century now! This is the one that persuaded me to do it myself. It uses code generation to produce efficient but ugly and verbose c code from equations in a config file. But there is also TIDES, ATOMFT, and many others, all slightly different approaches.
My approach is just "code it by hand, it's simpler".
Don't let that stop you though. It appears you have some worthy code that people should know about and a PF Insight article would help promote its use.
I'll beef up the technical note this weekend to cover the hamiltonian solver, and see how I feel then.

BTW the code I posted can be made much more efficient by storing the hamiltonian derivatives outside the "derivative jet" loop, like this.

Last edited:
One thing I forgot to mention in the OP, I have previously done my hamiltonian simulations using symplectic integrators, but the only types simple enough for me to program myself (Stormer-Verlet, Yoshida, Suzuki and derivatives) were only effective on "separable" hamiltonians. The method in this thread does not have this restriction, as should be clear from the example in question ;)

Last edited:
scottdave
Homework Helper
This is all awesome. I just started to learn Python, last year (for data science). It's cool to see use in dynamic systems..

This is all awesome. I just started to learn Python, last year (for data science). It's cool to see use in dynamic systems..
To me it is surprising ;)
I first started in automatic differentiation by reverse engineering the Lorenz equation Taylor code in c. I originally intended to do a scripted version in Julia, but the one-based indexing, runtime startup delay, and the general ""difficulty" of the language meant I did it in Python instead.
Now I have expanded beyond looking at chaotic ODEs, the operator overloading features of Python have now persuaded me to freeze development of the c version.

jedishrfu
Mentor
You should still consider doing it in Julia. Its becoming a major force in computer simulation because of its similarity to Matlab and because its free - free free free (turbo tax commercial)

Yes, Julia has a penalty for compiling code at the start but benchmarks indicate it runs much faster than many other competing languages.

Also you might find more interest for your work too.

You should still consider doing it in Julia. Its becoming a major force in computer simulation because of its similarity to Matlab and because its free - free free free (turbo tax commercial)
Well, the work is now completed (twice), no point doing it for a third time. In any case, there is a lot of overlap with this Julia package. My intention was to produce code that was short enough for me and others to tinker with, and learn something new.

The Taylor Series recurrences are sheer madness with one-based indexing, believe me! I never finished them, and it wasn't any fun trying (I am not a professional!).

BTW the python implementation is already faster than the c one for high (say 2000th) order Taylor Series ;)

I would much rather attempt an insights article than punish myself further with Julia!

Last edited:
OK, I'm posting a working outline of the technical description of my method, hopefully this will also notify @Greg Bernhardt, but I would like to canvas any bored mathematicians to point out anything that I can change to make this more rigorous.

One aspect of this work that I do not intend to document myself is the way the Taylor Coefficient recurrence relations are defined, this is done much better elsewhere (open access, look for the download link).

So now you know how it all comes together!

For clarification, this is based on a very old technique. My "contribution", such as it is, is limited to:
• Making a simple Python "playground" for learning about and experimenting with the method in a tiny codebase - also it runs quite nicely on a Raspberry Pi class computer
• I think the Hamiltonian analysis that combines the two forms of Automatic Differentiation (Taylor Series Method and Dual Numbers) is probably obscure, I haven't seen it anywhere else.

#### Attachments

• 115.3 KB Views: 83
Last edited:
Hmmm, energy creep, wow :(
Seems to be related only to the step size, and not the integrator order.
Looks like I need to abandon this idea and use a symplectic integrator instead. I already have one, but I need to make it implicit because the double pendulum is a non-separable hamiltonian (I suppose I can at least then use it for metrics/geodesic analysis in General Relativity).
Is there no end to this? ;)