Search code examples
pythonsympyphysicsdifferentiationorbital-mechanics

Sympy implicit differentiation for euler-lagrange


I'm trying to perform implicit differentiation on some expressions using sympy's idiff function.

In this case, rdot is dr/ds where s is an affine parameter. I want to perform implicit differentiation on Ltdot,Lphidot and Lrdot with respect to the same affine parameter s.

import numpy as np
from sympy import *
from sympy.physics.mechanics import *

#definition of variables
s = dynamicsymbols('s')
r = Function('r')(s)
rdot = Function('rdot')(s)
t = Function('t')(s)
tdot = Function('tdot')(s)
phi = Function('phi')(s)
phidot = Function('phidot')(s)


def F(x):
    return 1-(1/x)


# Largangian
def L(a,b,c, adot, bdot, cdot, photon = true): #r,t,phi
    return F(a)*(bdot)**2 - adot**2/F(a) - (a*cdot)**2


L = L(r, t, phi, rdot, tdot, phidot, photon = True)
Lt = diff(L, t)
Ltdot = diff(L, tdot)
Lphi = diff(L, phi)
Lphidot = diff(L, phidot)
Lr = diff(L, r)
Lrdot = diff(L, rdot)


#E-L equations printed to be used to solve equations
print('d/ds(', Ltdot, ') =', Lt) #EL1
print('d/ds(', Lphidot, ') =', Lphi) #EL2
print('d/ds(', Lrdot, ') =', Lr) #EL3


#FIX THIIISSSSS------------------------------------------------
LHS_EL1 = idiff(Ltdot, [t, tdot], s)
LHS_EL2 = idiff(Lphidot, [phi, phidot], s)
LHS_EL3 = idiff(Lrdot, [r, rdot], s)
#i want to do implicit differentiation wrt to affine parameter s, same that r is differentiated by to make rdot!!


print('d/ds(', LHS_EL1, ') =', Lt) #EL1 finalised
print('d/ds(', LHS_EL2, ') =', Lphi) #EL2 finalised
print('d/ds(', LHS_EL3, ') =', Lr) #EL3 finalised

I get the following error msg:

Traceback (most recent call last):
File "/Users/myname/PycharmProjects/untitled/.idea/14.1.py", line 53, in <module>
LHS_EL1 = idiff(Ltdot, [t, tdot], s)
File "/Users/myname/PycharmProjects/untitled/venv/lib/python3.6/site-packages/sympy/geometry/util.py", line 589, in idiff
yp = solve(eq.diff(x), dydx)[0].subs(derivs)
IndexError: list index out of range

Any thoughts on how I could achieve what I want or any help debugging would be greatly appreciated!


Solution

  • It's a little confusing to have the implicit 't' as the "time" variable for s and t be the function t(s). When you differentiate wrt t do you mean "Function('t')" or "s.args[0]"? If the latter, then if T = s.args[0] then

    >>> diff(L, T)
    2*(1 - 1/r(s(t)))*tdot(s(t))*Derivative(s(t),
    t)*Derivative(tdot(s(t)), s(t)) -
    2*phidot(s(t))**2*r(s(t))*Derivative(r(s(t)), s(t))*Derivative(s(t),
    t) - 2*phidot(s(t))*r(s(t))**2*Derivative(phidot(s(t)),
    s(t))*Derivative(s(t), t) + tdot(s(t))**2*Derivative(r(s(t)),
    s(t))*Derivative(s(t), t)/r(s(t))**2 -
    2*rdot(s(t))*Derivative(rdot(s(t)), s(t))*Derivative(s(t), t)/(1 -
    1/r(s(t))) + rdot(s(t))**2*Derivative(r(s(t)), s(t))*Derivative(s(t),
    t)/((1 - 1/r(s(t)))**2*r(s(t))**2)