# SymPy, jacobians, and Euler-Lagrange?

**URL:** https://discourse.julialang.org/t/sympy-jacobians-and-euler-lagrange/53666
**Category:** New to Julia
**Created:** [January 20, 2021, 1:03pm UTC](https://discourse.julialang.org/t/sympy-jacobians-and-euler-lagrange/53666 "2021-01-20T13:03:16Z")
**Posts on this page:** 1
**Page:** 1

<div class="post-metadata">

### Author: ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)
#### Post date: [January 20, 2021, 1:03pm UTC](https://discourse.julialang.org/t/sympy-jacobians-and-euler-lagrange/53666/1 "2021-01-20T13:03:16Z")

</div>

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}}:

```julia
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
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
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!