моделирование — ODE решатель из лагранжиана / вариационные методы в переполнении стека

У меня есть общий вопрос, который я сформулирую в контексте более конкретной ситуации.

Если кто-то хочет найти динамику двойного маятника, он может математически вывести уравнения движения, переписать ODE, чтобы иметь специальную форму, полезную для численных расчетов, и решить ODE, используя, скажем, odeint в c ++ (см. Пример этого переполнение стека https://stackoverflow.com/a/30582741).

Теперь представьте, что мы хотим сделать то же самое для n связанных маятников (n известно во время выполнения). Это требует, чтобы мы написали так называемый лагранжиан (кинетическая энергия — потенциальная энергия), и различные производные этой функции будут теми ОДУ, которые нам нужно решить. Кроме того, эти ODE должны быть переписаны в форме, подходящей для odeint. Это может быть трудно сделать вручную.

В таких программах, как Mathematica и Maple, это на самом деле довольно просто. Можно явным образом вывести необходимые ОДУ из лагранжиана, и для решения ОДУ не требуется, чтобы мы приводили уравнения в какой-либо специальной форме (см. Пример в Mathematica здесь https://mathematica.stackexchange.com/a/84279).

Можно ли сделать это в C ++, не переживая слишком много проблем?

Потенциальный подход:

Потенциальным подходом может быть использование пакета с ++ ginac. Это может помочь нам в аналитическом выводе ODE. Но я не знаю, как переписать выражения, поступающие из ginac, в форму, подходящую для численных расчетов в odeint. Есть идеи?