3. Advanced derivatives

3.1. Derivatives with respect to symbolic functions

Derivatives can be taken with respect to any of the following:

  • Variables: x = sym.symbols('x')

  • Symbolic function invocations: f_x = sym.Function('f')(x)

  • Derivatives of symbolic function invocations: df_x = f_x.diff(x)

For instance, consider the computation of the Lagrangian of a simple pendulum. We assume the generalized coordinate of the system is the pendulum angle \(\theta\). We take \(\theta = 0\) to mean that the pendulum extends rightwards along the x-axis. The potential energy of the system is given by:

\[U = m \cdot g \cdot r \cdot \sin{\theta}\]
from wrenfold import sym

m, g, r, t = sym.symbols("m, g, r, t")

# Theta as a function of time:
theta = sym.Function("theta")(t)
theta_dot = theta.diff(t)
print(theta_dot)  # prints: Derivative(theta(t), t)

# Potential energy of the pendulum mass:
T = m * g * r * sym.cos(theta)
print(T)  # prints: g*m*r*cos(theta(t))

Where \(r\) is the length of the pendulum, \(m\) the mass of the system, and \(g\) the acceleration due to gravity. The kinetic energy of the system is given by:

\[\begin{split}\begin{equation*} \begin{split} V &= \frac{1}{2} \cdot m \cdot \lVert \mathbf{v}\rVert^2 \\ \mathbf{v} &= r\cdot\frac{\partial}{\partial t}\left(\cos\theta\cdot\hat{\mathbf{i}} + \sin\theta\cdot\hat{\mathbf{j}}\right) \end{split} \end{equation*}\end{split}\]
# Position of the pendulum mass:
p_mass = r * sym.vector(sym.cos(theta), sym.sin(theta))
# Velocity of the pendulum mass:
v_mass = p_mass.diff(t)
# Kinetic energy:
V = sym.rational(1, 2) * m * v_mass.squared_norm()

# Simplify the expression by substituting: cos(theta)^2 + sin(theta)^2 --> 1
V = (
    V.distribute()
    .collect([m, r, theta_dot])
    .subs((sym.cos(theta) ** 2 + sym.sin(theta) ** 2) / 2, sym.rational(1, 2))
)
print(V)  # prints: m*r**2*Derivative(theta(t), t)**2/2

We can then compute the Lagrangian \(L\), and the Euler-Lagrange equation of the system:

\[\begin{split}\begin{equation*} \begin{split} L &= T - V \\[10pt] 0 &= \frac{\partial}{\partial t}\left(\frac{\partial}{\partial\dot{\theta}}L\right) - \frac{\partial}{\partial\theta}L \\ 0 &= g\cdot\sin\theta - r\cdot\ddot{\theta} \end{split} \end{equation*}\end{split}\]
# Compute canonical momentum.
# Note that we take the derivative wrt `Derivative(theta(t), t)`.
L = T - V
q_theta = L.diff(theta_dot)

# And the Euler-Lagrange equation, which involves taking the derivative wrt `theta(t)`.
euler_lagrange = q_theta.diff(t) - L.diff(theta)
print(euler_lagrange)  # prints: g*m*r*sin(theta(t)) - m*r**2*Derivative(theta(t), t, 2)

Note that while computing the equation of motion, we took the derivative with respect to both \(\theta\left(t\right)\) and \(\dot{\theta}\left(t\right)\).

3.2. Truncated derivatives

It is sometimes useful to truncate a derivative to prevent it from propagating through some terms. For instance, you might wish to compute a weighting function for a particular residual without the weight itself contributing to the Jacobian. To that end, wrenfold includes the stop_derivative() expression:

from wrenfold import sym

c, x = sym.symbols('c, x')
f_x = sym.Function('f')(x)

# Geman-McClure weight function:
w = (c**2) / (c**2 + f_x**2)
print(sym.stop_derivative(w).diff(x)) # prints: 0

# Weighted residual:
residual = sym.stop_derivative(w) * f_x
print(residual.diff(x)) # prints: Derivative(f(x), x)*StopDerivative(c**2/(c**2 + f(x)**2))