In previous chapters we modeled robot dynamics using ordinary differential equations (ODEs) of the form

\[ \dot{\mathbf{x}}(t) = f\big(\mathbf{x}(t), \mathbf{u}(t), \boldsymbol{\theta}\big), \]

where \( \mathbf{x} \) collects joint positions and velocities, \( \mathbf{u} \) contains input torques/forces, and \( \boldsymbol{\theta} \) are physical parameters (masses, inertias, link lengths, friction coefficients, etc.). Differentiable physics means that both the dynamics \( f \) and the numerical integration scheme are treated as differentiable maps, so that trajectory-level functionals such as

\[ \mathcal{L}(\boldsymbol{\theta}, \mathbf{u}_{0:N-1}) = \sum_{k=0}^{N-1} \ell_k(\mathbf{x}_k, \mathbf{u}_k, \boldsymbol{\theta}) + \ell_N(\mathbf{x}_N, \boldsymbol{\theta}) \]

admit gradients \( \nabla_{\boldsymbol{\theta}} \mathcal{L} \) and \( \nabla_{\mathbf{u}_{0:N-1}} \mathcal{L} \). These gradients can be used for:

• parameter identification and model calibration (refining \( \boldsymbol{\theta} \)),
• co-design (simultaneously tuning mechanical parameters and trajectories),
• as a building block for trajectory optimization and learning-based methods (in later courses).

A typical computational loop for differentiable physics in robotics is:

flowchart TD
  P["Model parameters theta"] --> F["Dynamics map f(x,u,theta)"]
  F --> INT["Time integrator (discrete simulator)"]
  INT --> TRAJ["State trajectory x_0,...,x_N"]
  TRAJ --> LOSS["Trajectory loss L(theta,u_0,...,u_{N-1})"]
  LOSS --> GRAD["Gradients 'dL/dtheta', 'dL/du_k'"]
  GRAD --> OPT["Gradient-based update of theta or u_k"]
        

The key mathematical question of this lesson is: how do we compute these gradients efficiently and accurately, given a robot dynamics model and an integrator?

We operate in discrete time with sampling period \( \Delta t \). Let \( \mathbf{x}_k \approx \mathbf{x}(k\Delta t) \) and \( \mathbf{u}_k \approx \mathbf{u}(k\Delta t) \). A one-step numerical integrator defines a map

\[ \mathbf{x}_{k+1} = \Phi\big(\mathbf{x}_k, \mathbf{u}_k, \boldsymbol{\theta}\big), \quad k = 0,\dots,N-1, \]

where \( \Phi \) depends on the choice of integration scheme. For example, explicit Euler gives

\[ \mathbf{x}_{k+1} = \mathbf{x}_k + \Delta t \, f\big(\mathbf{x}_k,\mathbf{u}_k,\boldsymbol{\theta}\big). \]

For robot manipulators, a common continuous-time dynamics model in joint coordinates is

\[ \mathbf{M}(\mathbf{q},\boldsymbol{\theta}) \ddot{\mathbf{q}} + \mathbf{C}(\mathbf{q},\dot{\mathbf{q}},\boldsymbol{\theta}) \dot{\mathbf{q}} + \mathbf{g}(\mathbf{q},\boldsymbol{\theta}) + \mathbf{F}_d(\dot{\mathbf{q}},\boldsymbol{\theta}) = \boldsymbol{\tau}(t), \]

where \( \boldsymbol{\tau}(t) \) (written as \tau in LaTeX) collects joint torques. Defining \( \mathbf{x} = [\mathbf{q}^\top, \dot{\mathbf{q}}^\top]^\top \), we can write \( \dot{\mathbf{x}} = f(\mathbf{x},\boldsymbol{\tau},\boldsymbol{\theta}) \), and thus a differentiable \( \Phi \).

If \( f \) is continuously differentiable in its arguments, and the integrator uses a finite composition of additions, multiplications, and smooth functions, then \( \Phi \) is also continuously differentiable. This ensures that all Jacobians needed for gradient computation exist:

\[ \mathbf{A}_k := \frac{\partial \Phi}{\partial \mathbf{x}_k}(\mathbf{x}_k,\mathbf{u}_k,\boldsymbol{\theta}),\quad \mathbf{B}_k := \frac{\partial \Phi}{\partial \boldsymbol{\theta}}(\mathbf{x}_k,\mathbf{u}_k,\boldsymbol{\theta}),\quad \mathbf{G}_k := \frac{\partial \Phi}{\partial \mathbf{u}_k}(\mathbf{x}_k,\mathbf{u}_k,\boldsymbol{\theta}). \]

Consider a trajectory loss functional

\[ \mathcal{L}(\boldsymbol{\theta}) = \sum_{k=0}^{N-1} \ell_k(\mathbf{x}_k,\boldsymbol{\theta}) + \ell_N(\mathbf{x}_N,\boldsymbol{\theta}), \quad \mathbf{x}_{k+1} = \Phi(\mathbf{x}_k,\mathbf{u}_k,\boldsymbol{\theta}), \]

where we treat inputs \( \mathbf{u}_k \) as fixed. Define the sensitivity matrices with respect to parameters:

\[ \mathbf{S}_k := \frac{\partial \mathbf{x}_k}{\partial \boldsymbol{\theta}} \in \mathbb{R}^{n_x \times n_\theta}. \]

Differentiating the discrete dynamics recursion yields

\[ \mathbf{S}_{k+1} = \mathbf{A}_k \mathbf{S}_k + \mathbf{B}_k, \qquad \mathbf{S}_0 = \frac{\partial \mathbf{x}_0}{\partial \boldsymbol{\theta}}. \]

For explicit Euler,

\[ \mathbf{A}_k = \mathbf{I} + \Delta t \frac{\partial f}{\partial \mathbf{x}_k},\quad \mathbf{B}_k = \Delta t \frac{\partial f}{\partial \boldsymbol{\theta}}. \]

By the chain rule,

\[ \frac{\mathrm{d}\mathcal{L}}{\mathrm{d}\boldsymbol{\theta}} = \sum_{k=0}^{N-1} \left( \frac{\partial \ell_k}{\partial \mathbf{x}_k} \mathbf{S}_k + \frac{\partial \ell_k}{\partial \boldsymbol{\theta}} \right) + \frac{\partial \ell_N}{\partial \mathbf{x}_N} \mathbf{S}_N + \frac{\partial \ell_N}{\partial \boldsymbol{\theta}}. \]

Complexity. Forward sensitivity scales like \( \mathcal{O}(N n_x n_\theta) \), which is acceptable if the number of parameters \( n_\theta \) is small, but expensive when many parameters are calibrated. Differentiable physics often uses adjoint methods to reverse this scaling.

We derive the adjoint (discrete-time reverse-mode) method that computes \( \nabla_{\boldsymbol{\theta}} \mathcal{L} \) with cost \( \mathcal{O}(N n_x) \), independent of \( n_\theta \). Introduce Lagrange multipliers \( \boldsymbol{\lambda}_{k+1} \) for the dynamic constraints:

\[ \tilde{\mathcal{L}} = \sum_{k=0}^{N-1} \ell_k(\mathbf{x}_k,\boldsymbol{\theta}) + \ell_N(\mathbf{x}_N,\boldsymbol{\theta}) + \sum_{k=0}^{N-1} \boldsymbol{\lambda}_{k+1}^\top \big( \Phi(\mathbf{x}_k,\mathbf{u}_k,\boldsymbol{\theta}) - \mathbf{x}_{k+1} \big). \]

Taking variations with respect to \( \mathbf{x}_k \) and collecting terms yields the adjoint recursion. For \( k = N \) we have

\[ \frac{\partial \tilde{\mathcal{L}}}{\partial \mathbf{x}_N} = \frac{\partial \ell_N}{\partial \mathbf{x}_N} - \boldsymbol{\lambda}_N^\top = \mathbf{0}^\top \quad\Rightarrow\quad \boldsymbol{\lambda}_N^\top = \frac{\partial \ell_N}{\partial \mathbf{x}_N}. \]

For intermediate steps \( k = N-1,\dots,0 \), using the chain rule gives

\[ \frac{\partial \tilde{\mathcal{L}}}{\partial \mathbf{x}_k} = \frac{\partial \ell_k}{\partial \mathbf{x}_k} + \boldsymbol{\lambda}_{k+1}^\top \mathbf{A}_k - \boldsymbol{\lambda}_k^\top = \mathbf{0}^\top, \]

hence the backward recursion

\[ \boldsymbol{\lambda}_k^\top = \frac{\partial \ell_k}{\partial \mathbf{x}_k} + \boldsymbol{\lambda}_{k+1}^\top \mathbf{A}_k, \quad k = N-1,\dots,0. \]

Finally, differentiating \( \tilde{\mathcal{L}} \) with respect to \( \boldsymbol{\theta} \) and using the dynamics constraints gives

\[ \frac{\mathrm{d}\mathcal{L}}{\mathrm{d}\boldsymbol{\theta}} = \sum_{k=0}^{N-1} \left( \frac{\partial \ell_k}{\partial \boldsymbol{\theta}} + \boldsymbol{\lambda}_{k+1}^\top \mathbf{B}_k \right) + \frac{\partial \ell_N}{\partial \boldsymbol{\theta}}. \]

The computational cost is dominated by the forward rollout and the backward adjoint recursion, both scaling with \( N n_x^2 \) if \( \mathbf{A}_k \) is dense, but independent of \( n_\theta \).

flowchart TD
  X0["Forward pass: x_0"] --> FWD["Propagate x_{k+1} = Phi(x_k,u_k,theta)"]
  FWD --> TRAJ["Store trajectory x_0,...,x_N"]
  TRAJ --> LOSS["Compute loss terms l_k, l_N"]
  LOSS --> LAMBDA["Initialize lambda_N = dl_N/dx_N"]
  LAMBDA --> BWD["Backward pass: lambda_k = dl_k/dx_k + A_k^T lambda_{k+1}"]
  BWD --> GRAD["Accumulate gradient with B_k: dL/dtheta"]
        

Consider a single rotary joint with inertia parameter \( I \), viscous friction \( b \), stiffness \( k \), and input torque \( \tau(t) \). In continuous time:

\[ I \ddot{q} + b \dot{q} + k q = \tau(t). \]

Let the state be \( \mathbf{x} = [q,\dot{q}]^\top \). Then

\[ \dot{\mathbf{x}} = f(\mathbf{x},\tau,I) = \begin{bmatrix} \dot{q} \\ \dfrac{1}{I}\big(\tau - b \dot{q} - k q\big) \end{bmatrix}. \]

Using explicit Euler with step \( \Delta t \):

\[ \mathbf{x}_{k+1} = \mathbf{x}_k + \Delta t \, f(\mathbf{x}_k,\tau_k,I). \]

Suppose we have desired positions \( q_k^{\mathrm{ref}} \) and define the loss

\[ \mathcal{L}(I) = \frac{1}{2}\sum_{k=0}^N \big(q_k - q_k^{\mathrm{ref}}\big)^2 = \sum_{k=0}^N \ell_k(\mathbf{x}_k),\quad \ell_k(\mathbf{x}_k) = \frac{1}{2}(q_k - q_k^{\mathrm{ref}})^2. \]

We want \( \mathrm{d}\mathcal{L} / \mathrm{d}I \). The forward sensitivity matrix \( \mathbf{S}_k = \partial \mathbf{x}_k / \partial I \) obeys

\[ \mathbf{S}_{k+1} = \mathbf{S}_k + \Delta t \left( \frac{\partial f}{\partial \mathbf{x}_k} \mathbf{S}_k + \frac{\partial f}{\partial I} \right), \quad \mathbf{S}_0 = \mathbf{0}. \]

Here

\[ \frac{\partial f}{\partial I} = \begin{bmatrix} 0 \\ - \dfrac{1}{I^2}\big(\tau_k - b \dot{q}_k - k q_k\big) \end{bmatrix}. \]

The gradient is then

\[ \frac{\mathrm{d}\mathcal{L}}{\mathrm{d}I} = \sum_{k=0}^N \frac{\partial \ell_k}{\partial \mathbf{x}_k} \mathbf{S}_k = \sum_{k=0}^N (q_k - q_k^{\mathrm{ref}})\, [1\;\;0]\mathbf{S}_k. \]

Alternatively, the adjoint method introduces \( \lambda_k = \partial \mathcal{L} / \partial \mathbf{x}_k \), yielding a backward recursion and

\[ \frac{\mathrm{d}\mathcal{L}}{\mathrm{d}I} = \sum_{k=0}^{N-1} \boldsymbol{\lambda}_{k+1}^\top \Delta t \frac{\partial f}{\partial I}. \]

Python offers both manual implementations of the adjoint method and automatic differentiation via libraries such as jax, torch, or robotics-oriented tools like pinocchio (for dynamics and kinematics) that can be combined with autodiff backends. Below is a minimal manual adjoint implementation for the 1-DOF example:

import numpy as np

def f(x, tau, I, b=0.1, k=1.0):
    q, dq = x
    ddq = (tau - b * dq - k * q) / I
    return np.array([dq, ddq])

def step_euler(x, tau, I, dt):
    return x + dt * f(x, tau, I)

def rollout(x0, taus, I, dt):
    xs = [x0]
    for tau in taus:
        xs.append(step_euler(xs[-1], tau, I, dt))
    return np.array(xs)  # shape: (N+1, 2)

def loss(xs, q_ref):
    # xs[:,0] contains q
    return 0.5 * np.sum((xs[:, 0] - q_ref)**2)

def grad_I_adjoint(x0, taus, I, dt, q_ref):
    xs = rollout(x0, taus, I, dt)
    N = len(taus)
    # adjoint lambdas: lambda_k in R^2
    lambdas = np.zeros_like(xs)
    # terminal condition
    lambdas[N, 0] = xs[N, 0] - q_ref[N]  # d l_N / d q_N
    # precompute A_k and dPhi/dI for explicit Euler
    A = np.zeros((N, 2, 2))
    dPhi_dI = np.zeros((N, 2))
    for k in range(N):
        q, dq = xs[k]
        tau = taus[k]
        # partial derivatives of f wrt x
        # df/dq = [0, -(k/I)]
        # df/ddq = [1, -(b/I)]
        df_dq = np.array([0.0, -1.0 * k / I])
        df_ddq = np.array([1.0, -1.0 * 0.1 / I])
        Jx = np.column_stack([df_dq, df_ddq])  # 2x2
        A[k] = np.eye(2) + dt * Jx
        # df/dI
        ddq = (tau - 0.1 * dq - k * q) / I
        dddq_dI = -(tau - 0.1 * dq - k * q) / (I**2)
        df_dI = np.array([0.0, dddq_dI])
        dPhi_dI[k] = dt * df_dI
    # backward adjoint recursion
    for k in range(N - 1, -1, -1):
        # gradient of stage loss wrt x_k
        dL_dxk = np.zeros(2)
        dL_dxk[0] = xs[k, 0] - q_ref[k]
        lambdas[k] = dL_dxk + A[k].T @ lambdas[k + 1]
    # parameter gradient
    dL_dI = 0.0
    for k in range(N):
        dL_dI += lambdas[k + 1] @ dPhi_dI[k]
    return dL_dI

# Example usage
dt = 0.01
N = 100
taus = np.ones(N) * 1.0
t_grid = np.arange(N + 1) * dt
q_ref = 0.5 * np.sin(2.0 * np.pi * t_grid)
x0 = np.array([0.0, 0.0])
I = 0.5

xs = rollout(x0, taus, I, dt)
print("Loss:", loss(xs, q_ref))
print("dL/dI (adjoint):", grad_I_adjoint(x0, taus, I, dt, q_ref))
      

With jax, we can implement the same dynamics in a functional style and let autodiff compute gradients through the entire rollout, effectively performing the adjoint method automatically while we retain full control over the dynamics model.

In C++, differentiable physics is typically implemented using:

• Robotics dynamics libraries (e.g., Pinocchio, RBDL) to compute \( f \), and
• algorithmic differentiation (AD) tools (e.g., CppAD, adept) to obtain Jacobians and gradients.

The idea is to express one simulation step Phi(x,u,theta) using AD types, then invoke the AD library to obtain dPhi/dtheta or reverse-mode gradients. A skeletal example with CppAD for the 1-DOF system is:

#include <iostream>
#include <vector>
#include <cppad/cppad.hpp>

using CppAD::AD;

template <typename Scalar>
std::vector<Scalar> f(const std::vector<Scalar>& x,
                       Scalar tau, Scalar I,
                       Scalar b, Scalar k)
{
    Scalar q  = x[0];
    Scalar dq = x[1];
    Scalar ddq = (tau - b * dq - k * q) / I;
    std::vector<Scalar> dx(2);
    dx[0] = dq;
    dx[1] = ddq;
    return dx;
}

template <typename Scalar>
std::vector<Scalar> stepEuler(const std::vector<Scalar>& x,
                               Scalar tau, Scalar I,
                               Scalar b, Scalar k, Scalar dt)
{
    std::vector<Scalar> dx = f(x, tau, I, b, k);
    std::vector<Scalar> xnext(2);
    xnext[0] = x[0] + dt * dx[0];
    xnext[1] = x[1] + dt * dx[1];
    return xnext;
}

int main()
{
    using ADScalar = AD<double>;

    // Independent variable: I (inertia parameter)
    std::vector<ADScalar> I_vec(1);
    I_vec[0] = 0.5;
    CppAD::Independent(I_vec);

    ADScalar I = I_vec[0];
    ADScalar b = 0.1;
    ADScalar k = 1.0;
    ADScalar dt = 0.01;
    ADScalar tau = 1.0;

    std::vector<ADScalar> x(2);
    x[0] = 0.0; // q
    x[1] = 0.0; // dq

    // Simulate N steps and accumulate loss
    const int N = 100;
    ADScalar L = 0.0;
    for (int kstep = 0; kstep <= N; ++kstep) {
        ADScalar q_ref = 0.5 * CppAD::sin(2.0 * M_PI * dt * kstep);
        ADScalar err = x[0] - q_ref;
        L += 0.5 * err * err;
        if (kstep < N) {
            x = stepEuler(x, tau, I, b, k, dt);
        }
    }

    // Declare dependent variable
    std::vector<ADScalar> L_vec(1);
    L_vec[0] = L;
    CppAD::ADFun<double> tape(I_vec, L_vec);

    // Compute dL/dI at I = 0.5
    std::vector<double> I_val(1), grad(1);
    I_val[0] = 0.5;
    grad = tape.Jacobian(I_val); // since output dimension is 1
    std::cout << "Loss gradient dL/dI = " << grad[0] << std::endl;

    return 0;
}
      

For multi-DOF robots, one typically uses a library like Pinocchio to compute the dynamics and wraps the state and parameters in AD types, letting the AD tool perform reverse-mode differentiation through the entire simulation.

Java has mature linear algebra libraries such as EJML and numerical libraries that can be combined with:

• manual implementation of forward/adjoint recursions, or
• automatic differentiation frameworks (e.g., ND4J, custom reverse-mode).

The following simplified example illustrates manual forward sensitivity for the 1-DOF system using plain Java arrays (for clarity). In practice one would use matrix classes from EJML:

public class DifferentiablePendulum {

    static double[] f(double[] x, double tau, double I,
                      double b, double k) {
        double q = x[0];
        double dq = x[1];
        double ddq = (tau - b * dq - k * q) / I;
        return new double[]{dq, ddq};
    }

    static double[] stepEuler(double[] x, double tau, double I,
                              double b, double k, double dt) {
        double[] dx = f(x, tau, I, b, k);
        return new double[]{
            x[0] + dt * dx[0],
            x[1] + dt * dx[1]
        };
    }

    public static void main(String[] args) {
        double dt = 0.01;
        int N = 100;
        double I = 0.5;
        double b = 0.1;
        double k = 1.0;
        double tau = 1.0;

        double[] x = new double[]{0.0, 0.0};
        double[][] xTraj = new double[N + 1][2];
        xTraj[0] = x.clone();

        // Forward rollout
        for (int kstep = 0; kstep < N; ++kstep) {
            x = stepEuler(x, tau, I, b, k, dt);
            xTraj[kstep + 1] = x.clone();
        }

        // Simple reference trajectory
        double[] qRef = new double[N + 1];
        for (int kstep = 0; kstep <= N; ++kstep) {
            double t = dt * kstep;
            qRef[kstep] = 0.5 * Math.sin(2.0 * Math.PI * t);
        }

        // Forward sensitivity S_k = d x_k / dI (2x1 vector)
        double[][] S = new double[N + 1][2];
        S[0][0] = 0.0;
        S[0][1] = 0.0;

        for (int kstep = 0; kstep < N; ++kstep) {
            double q = xTraj[kstep][0];
            double dq = xTraj[kstep][1];

            // Jacobian of f wrt x
            double df_dq0 = 0.0;
            double df_dq1 = -k / I;
            double df_ddq0 = 1.0;
            double df_ddq1 = -b / I;

            // Jx * S_k
            double JxS0 = df_dq0 * S[kstep][0] + df_ddq0 * S[kstep][1];
            double JxS1 = df_dq1 * S[kstep][0] + df_ddq1 * S[kstep][1];

            // df/dI
            double num = tau - b * dq - k * q;
            double dddq_dI = -num / (I * I);
            double df_dI0 = 0.0;
            double df_dI1 = dddq_dI;

            S[kstep + 1][0] = S[kstep][0] + dt * (JxS0 + df_dI0);
            S[kstep + 1][1] = S[kstep][1] + dt * (JxS1 + df_dI1);
        }

        // Gradient dL/dI
        double dLdI = 0.0;
        for (int kstep = 0; kstep <= N; ++kstep) {
            double q = xTraj[kstep][0];
            double err = q - qRef[kstep];
            dLdI += err * S[kstep][0]; // only q-component
        }

        System.out.println("Gradient dL/dI = " + dLdI);
    }
}
      

This implementation mirrors the continuous-time sensitivity derivation and shows that Java can host differentiable physics computations using only basic language features plus a linear algebra package.

In MATLAB, differentiable physics can be approached in (at least) two complementary ways:

1. Symbolic differentiation of the dynamics using the Symbolic Math Toolbox, then generating efficient numeric functions for simulation and gradients.
2. Simulink-based sensitivity analysis, using Simscape Multibody for the physics and Simulink Design Optimization / linearization tools to compute sensitivities of outputs to physical parameters.

A symbolic example for the 1-DOF system:

syms q dq I b k tau real
x = [q; dq];

f = [dq;
     (tau - b*dq - k*q)/I];

% Jacobians needed for sensitivities
A = jacobian(f, x);   % df/dx
B = jacobian(f, I);   % df/dI

disp('A(q,dq,I) = ');
disp(A);
disp('B(q,dq,I) = ');
disp(B);

% Generate numeric functions for use in simulation loops
matlabFunction(f, 'File', 'f_pendulum', 'Vars', {q, dq, tau, I, b, k});
matlabFunction(A, 'File', 'A_pendulum', 'Vars', {q, dq, tau, I, b, k});
matlabFunction(B, 'File', 'B_pendulum', 'Vars', {q, dq, tau, I, b, k});
      

These generated functions can be called from a numeric script that unrolls the dynamics and the forward or adjoint recursions. For multi-body robots, one can use the Robotics System Toolbox or custom symbolic derivations of \( \mathbf{M}(\mathbf{q}) \), \( \mathbf{C}(\mathbf{q},\dot{\mathbf{q}}) \), and \( \mathbf{g}(\mathbf{q}) \), and differentiate them with respect to inertial parameters.

In Simulink, differentiable physics is realized by constructing a Simscape Multibody model of the robot, marking physical parameters (masses, inertias) as tunable, and using built-in tools to compute gradients of cost functionals (tracking error, energy, etc.) with respect to these parameters. This is conceptually equivalent to the discrete adjoint method derived earlier.

Wolfram Mathematica is well suited for differentiable physics because it treats symbolic and numeric computations in a unified framework. One can define the dynamics symbolically, solve or numerically integrate the system, and then differentiate a loss functional with respect to parameters using D or Grad.

(* 1-DOF rotary system *)
ClearAll["Global`*"];
q[t_]; dq[t_];

params = {I, b, k};
tau[t_] := 1.0;

eq = {
  q'[t] == dq[t],
  dq'[t] == (tau[t] - b*dq[t] - k*q[t])/I
};

ic = {q[0] == 0, dq[0] == 0};
tmax = 1.0;

sol = NDSolve[
  Join[eq, ic],
  {q, dq},
  {t, 0, tmax},
  Method -> "StiffnessSwitching"
][[1]];

qFun[t_, I_?NumericQ, b_?NumericQ, k_?NumericQ] :=
  Evaluate[q[t] /. sol /. {I -> I, b -> b, k -> k}];

(* Reference and loss functional *)
qRef[t_] := 0.5*Sin[2*Pi*t];

L[I_, b_, k_] := NIntegrate[
  0.5*(qFun[t, I, b, k] - qRef[t])^2,
  {t, 0, tmax}
];

(* Gradient with respect to parameters *)
gradL[I0_, b0_, k0_] :=
  Grad[L[I, b, k], {I, b, k}] /. {I -> I0, b -> b0, k -> k0};

(* Example evaluation *)
gradVal = gradL[0.5, 0.1, 1.0];
gradVal
      

For higher-dimensional robots, Mathematica can symbolically derive the manipulator equations via Lagrange's method, automatically compute Jacobians and sensitivity equations, and export optimized C or C++ code for embedded differentiable physics computations.

In this lesson we reframed robot dynamics and numerical integrators as differentiable maps, enabling analytic gradient computation of trajectory-level losses with respect to both physical parameters and inputs. We derived forward sensitivity equations and the discrete adjoint (reverse-mode) recursion, highlighting their computational trade-offs. Through a 1-DOF rotary example, we showed how these gradients can be implemented manually and via automatic differentiation in multiple programming languages and environments (Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica). Differentiable physics thus provides a rigorous bridge between classical robot modeling and modern gradient-based methods for identification and design.