Я представлю выполнение формализма Эйлера-Лагранжа для примера двойного маятника как нетривиальный, полный, стандартный пример, и это без использования специализированных функций в sympy.physics, только с использованием basi c средства дифференциации и написания кода от sympy. Надеемся, что он достаточно общий, так как вторая часть должна быть независимой от проблем и, следовательно, также непосредственно применимой к вашей ситуации.
from sympy import sin, cos, Symbol, symbols, solve
from sympy.utilities.lambdify import lambdify
Лагранжиан для физической модели
Сначала физическая установка для Лагранжиан использует два угла в качестве основных зависимых функций и строит кинетику c и потенциальную энергию через декартовы координаты.
# Symbols for the parameters of the problem
t,m1,m2,l1,l2,g = symbols("t,m_1 m_2 l_1 l_2 g")
# Variables of the problem
th1, th2 = Function("θ_1")(t), Function("θ_2")(t)
x1,y1 = l1*sin(th1), -l1*cos(th1)
x2,y2 = x1+l2*sin(th2), y1-l2*cos(th2)
# kinetic energy
vx1,vy1,vx2,vy2 = ( xx.diff(t) for xx in (x1,y1,x2,y2))
K1 = m1/2*(vx1**2+vy1**2)
K2 = m2/2*(vx2**2+vy2**2)
K = K1+K2
# potential energy
V = g*(m1*y1+m2*y2)
# Lagrangian
L = K - V
L = L.expand().simplify()
Чтобы получить абстрактную обработку, используйте массивы абстрактных параметров и векторы координат
params = [m1, l1, m2, l2, g]
q = [th1, th2]
dotq = [ qq.diff(t) for qq in q]
От Эйлера-Лагранжа до системы ОДУ первого порядка
В полученных выражениях становится довольно грязно, если перейти ко вторым производным от углов, можно получить более структурированный подход и, надеюсь, соответственно более быстрая оценка с использованием импульсных переменных
pk = diff(L, dotqk)
d(pk)/dt = diff(L, qk)
, где первое соотношение рассматривается как система уравнений для вычисления dotq из p.
Подготовьте новые переменные, приготовьтесь заменить функциональные символы простыми переменными.
N = len(q)
p = [ Symbol(f"p_{k+1}") for k in range(N)]
dotq_func, dotq = dotq, [ Symbol(f"Dq_{k+1}") for k in range(N)]
q_func, q = q, [ Symbol(f"q_{k+1}") for k in range(N)]
Теперь замените все функциональные члены простыми переменными
L=L.subs(list(zip(dotq_func, dotq))).subs(list(zip(q_func, q)))
Теперь устанавливается sh функция для вычисления dotq из q,p
p_eqns = [ p[k] - L.diff(dotq[k]) for k in range(N)]
dotq_expr = solve(p_eqns, dotq)
dotq_func = lambdify([*q, *p, *params],[ dotq_expr[dq] for dq in dotq])
Далее генерируется функция, которая вычисляет производные p
dotp_expr = [ L.diff(q[k]) for k in range(N)]
dotp_func = lambdify([*q,*dotq,*params],dotp_expr)
Теперь соберите сгенерированные функции в полную функцию ODE и решите с ней тестовую задачу, чтобы убедиться, что она работает
def odefunc(t,u,args):
q, p = u[:N], u[N:]
dotq = dotq_func(*q, *p, *args)
dotp = dotp_func(*q, *dotq, *args)
return [*dotq, *dotp]
Подтверждение с помощью численного решения
myparams = [1,10,5,5,9.81]
t = np.linspace(0,25,301)
u = odeint(odefunc,[2,1.2, 0,0],t,args=(myparams,), tfirst=True)
%matplotlib inline
plt.figure(figsize=(8,5))
plt.plot(t,u[:,0],t,u[:,1]);
plt.grid(); plt.show()
В результате получается следующий график для функций угла 
Outlook
Можно ожидать даже более высокой производительности, если использовать термин kineti c (как есть обычно в классической механике) гарантированно будет иметь квадратичные c скорости. Затем, непосредственно извлекая матрицу этой квадратичной c формы, можно было бы делегировать системное решение в преобразовании импульсов в скорость числовому линейному решателю системы, поэтому не сохраняйте выражения symboli c для обращения матрицы. Они могут быть большими в более высоких размерах.