Python example equivalent to: https://gist.github.com/geoffreygarrett/019d65059f0a2966f0db3158dfd2d03d

example_autodiff_stts.py

from sympy import symbols, Matrix, diff, lambdify
import numpy as np


# Compute the 1st-order State Transition Tensor
def dfdx(eom, wrt, at):
    u = eom.subs(at).evalf()
    m = len(u) if hasattr(u, "__len__") else 1
    n = len(wrt)
    t2 = np.zeros((m, n))

    # Loop over each variable to find partial derivatives
    for j, xj in enumerate(wrt):
        partial_diff = diff(eom, xj)
        evaluated_diff = partial_diff.subs(at).doit().evalf()
        t2[:, j] = np.array(evaluated_diff).flatten()

    return t2


# Compute the 2nd-order State Transition Tensor
def d2fdx2(eom, wrt, at):
    # Initialize
    m = np.size(eom)
    n = len(wrt)
    t3 = np.zeros((m, n, n))

    # Compute partial derivatives
    for j, xj in enumerate(wrt):
        for k, xk in enumerate(wrt):
            partial_diff = diff(diff(eom, xj), xk)
            evaluated_diff = partial_diff.subs(at).doit().evalf()
            t3[:, j, k] = np.array(evaluated_diff).flatten()

    return t3


# Compute the 3rd-order State Transition Tensor
def d3fdx3(eom, wrt, at):
    # Initialize
    m = np.size(eom)
    n = len(wrt)
    t4 = np.zeros((m, n, n, n))

    # Compute partial derivatives
    for j, xj in enumerate(wrt):
        for k, xk in enumerate(wrt):
            for l, xl in enumerate(wrt):
                partial_diff = diff(diff(diff(eom, xj), xk), xl)
                evaluated_diff = partial_diff.subs(at).doit().evalf()
                t4[:, j, k, l] = np.array(evaluated_diff).flatten()

    return t4


# Define physical parameters and symbols
x, y, z, vx, vy, vz, GM, CD, m, S, rho = symbols("x y z vx vy vz GM CD m S rho", real=True)
v = Matrix([vx, vy, vz])
r = Matrix([x, y, z])
accel = -GM * r / r.norm() ** 3 - 0.5 * rho * CD * S / m * v * v.norm()
eom = Matrix.vstack(v, accel)

# Define variables to differentiate and point of evaluation
wrt = [x, y, z, vx, vy, vz]
at = {x: 1, y: 0, z: 0, vx: 0, vy: 1, vz: 0, GM: 1, CD: 1, m: 1, S: 1, rho: 1}

# Compute and display tensors
dfdx_tensor = dfdx(eom, wrt, at)
d2fdx2_tensor = d2fdx2(eom, wrt, at)
d3fdx3_tensor = d3fdx3(eom, wrt, at)

# Define the perturbation
perturbation = np.array([0.1, 0.1, 0.1, 0.1, 0.1, 0.1])

# Compute the terms using tensordot
term1 = np.tensordot(dfdx_tensor, perturbation, axes=1)
term2 = np.tensordot(np.tensordot(d2fdx2_tensor, perturbation, axes=1), perturbation, axes=1) / 2
term3 = np.tensordot(np.tensordot(np.tensordot(d3fdx3_tensor, perturbation, axes=1), perturbation, axes=1),
                     perturbation,
                     axes=1) / 6

print("Term 1:", term1)
print("Term 2:", term2)
print("Term 3:", term3)

# Term 1: [ 0.1   0.1   0.1   0.15 -0.2  -0.15]
# Term 2: [ 0.     0.     0.    -0.005  0.02   0.025]
# Term 3: [ 0.      0.      0.     -0.0085 -0.003  -0.0035]

Full example seen here for an arbitrary ISS-like orbit.