SymPy, jacobians, and Euler-Lagrange?

New to Julia
1

I’m trying to use SymPy to check some hand calculations regarding the use of the Euler-Lagrange equation. However, I get stuck because I cannot find functions for Jacobians, etc. in SymPy…

First, I define generalized coordinates \mathsf{x}, generalized velocity \mathsf{v}, mass matrix \mathsf{M}(\mathsf{x}), kinetic and potential energies, K = \frac{1}{2}\mathsf{v}^TM\mathsf{v} and P, respectively, then Lagrangian \mathsf{L} \triangleq K-P, and generalized potential \mathsf{p} \triangleq \frac{\partial \mathsf{L}}{\partial \mathsf{v}}:

using SymPy, LinearAlgebra
#
@vars J_bc J_s m r_c h g T_θ T_ψ positive=true
@vars θ ψ ω_θ ω_ψ
#
# Mass matrix
M_11 = J_bc + m*(r_c^2+h^2)
M_22 = J_bc*cos(θ)^2 + m*(r_c*cos(θ)-h*sin(θ))^2 + J_s
M = [M_11 0;0 M_22]
#
# coordinate and velocity
x = [θ, ψ]
v = [ω_θ, ω_ψ]
#
# Kinetic energy K, potential energy P
K = (permutedims(v)*M*v)[1]
P = m*g*(r_c*sin(θ) + h*cos(θ))
# 
# Lagrangian L and momentum p
L = K-P
p = M*v

In addition to the trivial ODE \frac{d \mathsf{x}}{dt} = \mathsf{v}, the Euler-Lagrange equations include either a DAE:

\frac{d \mathsf{p}}{dt} = \frac{\partial \mathsf{L}}{\partial \mathsf{x}} + \mathsf{F} \\ \mathsf{p} = \mathsf{M}\mathsf{x}

where \mathsf{F} is the generalized force, or — by eliminating \mathsf{p} — an implicit ODE

\mathsf{M}(\mathsf{x})\frac{d\mathsf{v}}{dt} = -\frac{\mathsf{\partial p}}{\partial \mathsf{x}}\mathsf{v} + \frac{\partial \mathsf{L}}{\partial \mathsf{x}} + \mathsf{F}

OK – I can compute the gradient \frac{\partial \mathsf{L}}{\partial \mathsf{x}}:

julia> simplify.(diff(L,x)) |> show
SymMatrix(PyObject Matrix([
[g*m*(h*sin(θ) - r_c*cos(θ)) + ω_ψ**2*(-J_bc*sin(2*θ) + h**2*m*sin(2*θ) - 2*h*m*r_c*cos(2*θ) - m*r_c**2*sin(2*θ))],
[                                                                                                               0]]))

But how can I compute the Jacobian \frac{\mathsf{\partial p}}{\partial \mathsf{x}}\mathsf{v} ?

  • I cannot find a “jacobian” function to compute \frac{\mathsf{\partial p}}{\partial \mathsf{x}}
  • Instead, I try to compute the individual rows as \frac{\mathsf{\partial p}_i}{\partial \mathsf{x}}\mathsf{v} 
julia> diff(p[1],x)
SymMatrix(PyObject Matrix([
[0],
[0]]))

julia> diff(p[2],x) |> show
SymMatrix(PyObject Matrix([
[ω_ψ*(-2*J_bc*sin(θ)*cos(θ) + m*(-h*sin(θ) + r_c*cos(θ))*(-2*h*cos(θ) - 2*r_c*sin(θ)))],
[                                                                                    0]]))

which works, but …
Question: How can I compute the scalar product of \frac{\mathsf{\partial p}_i}{\partial \mathsf{x}} and \mathsf{v}

  • using permutedims on the gradient doesn’t work, and
  • using the dot function from LinearAlgebra doesn’t work

Any suggestions are appreciated!