- #1
m4r35n357
- 654
- 148
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?
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, cosclass 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])
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) * 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.cosdef 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: