Pendulum-on-a-Cart Physics

import matplotlib.pyplot as plt
import numpy as np
import sympy as sp
from IPython.display import display

There are two masses in the system: the cart, and the pendulum bob (we'll assume the pendulum arm is massless).

We could treat the pendulum arm as a stiff spring, in which case Hooke's Law and Newton's Laws would give us the equations of motion of the system. But if we want the pendulum arm to be perfectly rigid, we have to treat it as a constraint. Lagrangian mechanics gives us a way to derive the equations of motion for a system with constraints.

Let the cart mass be $m_1$ kg, and the pendulum bob $m_2$ kg. We'll describe the position of the cart with a variable $x$ (positive $x$ to the right), and the angle of the pendulum with a variable $\theta$ (positive $\theta$ counter-clockwise). Let's define $v = \dot{x}$, and $\omega = \dot{\theta}$. Let $\theta = 0$ occur when the pendulum is hanging straight down. Finally, let's say the length of the pendulum arm is $l$.

def plot_pendulum_cart(x, theta, length = 1.0):
    """
    Generates a plot of the cart and pendulum.
    """
    fig, ax = plt.subplots()
    ax.set_xlim(-4.0, 4.0)
    ax.set_ylim(-1.25 * length, 1.25 * length)
    ax.set_aspect("equal")
    cart_pos = np.array([x, 0.0])
    pendulum_pos = cart_pos + length * np.array([np.sin(theta), -np.cos(theta)])
    ax.plot([cart_pos[0], pendulum_pos[0]], [cart_pos[1], pendulum_pos[1]], color="black")
    ax.scatter(cart_pos[0], cart_pos[1], color="black")
    ax.scatter(pendulum_pos[0], pendulum_pos[1], color="blue")

plot_pendulum_cart(-.5, 0.75 * np.pi)

To construct the Lagrangian, we need to know the potential and kinetic energy of the cart and pendulum bob.

Since the cart is constrained to move horizontally, its potential energy never changes. So we can call it zero.

$$ \begin{equation} V_\text{cart} = 0 \end{equation} $$

Its kinetic energy is almost as simple.

$$ \begin{equation} T_\text{cart} = \frac{1}{2} m_1 v^2 \end{equation} $$

Relative to the cart, the position of the pendulum is $(l \sin(\theta), -l \cos(\theta))$. So as the cart moves, the pendulum's coordinates are $(x + l \sin(\theta), -l \cos(\theta))$. Its potential energy is $-mgh$, where $g$ is the acceleration due to gravity (9.8 meters per second squared), and $h$ is its height.

$$ \begin{equation} V_\text{bob} = -m_2 gl \cos(\theta) \end{equation} $$

To compute the pendulum's kinetic energy, we need to know the magnitude of its velocity. Differentiating its coordinates with respect to time, we find that the velocity of the pendulum bob is $(v + \omega l \cos(\theta), \omega l \sin(\theta))$.

$$ \begin{equation} \begin{aligned} |v_\text{bob}|^2 &= (v + \omega l \cos(\theta))^2 + (\omega l \sin(\theta))^2 \\ &= v^2 + 2 v \omega l \cos(\theta) + \omega^2 l^2 \cos(\theta)^2 + \omega^2 l^2 \sin(\theta)^2 \\ &= v^2 + 2 v \omega l \cos(\theta) + \omega^2 l^2 \\ \end{aligned} \end{equation} $$$$ \begin{equation} T_\text{bob} = \frac{1}{2} m_2 |v_\text{bob}|^2 = \frac{1}{2} m_2 v^2 + m_2 v \omega l \cos(\theta) + \frac{1}{2} m_2 \omega^2 l^2 \end{equation} $$

From this we can construct the Lagrangian of our system: $T - V$.

$$ \begin{equation} \begin{aligned} \mathcal{L} &= T_\text{cart} + T_\text{bob} - V_\text{cart} - V_\text{bob} \\ &= \frac{1}{2} m_1 v^2 + \frac{1}{2} m_2 v^2 + m_2 v \omega l \cos(\theta) + \frac{1}{2} m_2 \omega^2 l^2 + m_2 gl \cos(\theta) \\ &= \frac{1}{2} (m_1 + m_2) v^2 + \frac{1}{2} m_2 \omega^2 l^2 + m_2 l (v \omega + g) \cos(\theta) \end{aligned} \end{equation} $$

The Euler-Lagrange equations give us the equations of motion in terms of the Lagrangian.

$$ \begin{equation} \begin{aligned} \frac{d}{dt} \frac{\partial \mathcal{L}}{\partial v} &= \frac{\partial \mathcal{L}}{\partial x} \\ \frac{d}{dt} \frac{\partial \mathcal{L}}{\partial \omega} &= \frac{\partial \mathcal{L}}{\partial \theta} \end{aligned} \end{equation} $$

If we just wanted to model the pendulum acting under the force of gravity, these equations would suffice. But to incorporate external forces, we need to add additional terms to the right hand side. We will need an external force $f$ that acts on the cart. And to make our simulation a little more realistic, we can also include damping terms: $-c_1 v$ for the cart, and $-c_2 \omega$ for the pendulum. With these terms, our modified Euler-Lagrange equations end up as follows.

$$ \begin{equation} \begin{aligned} \frac{d}{dt} \frac{\partial \mathcal{L}}{\partial v} &= \frac{\partial \mathcal{L}}{\partial x} + f -c_1 v \\ \frac{d}{dt} \frac{\partial \mathcal{L}}{\partial \omega} &= \frac{\partial \mathcal{L}}{\partial \theta} -c_2 \omega \end{aligned} \end{equation} $$

The algebra gets a little involved here, so let's use sympy.

# Constants
m1 = sp.symbols("m1", real=True)
m2 = sp.symbols("m2", real=True)
length = sp.symbols("l", real=True)
g = sp.symbols("g", real=True)
c1 = sp.symbols("c1", real=True) # cart damping coefficient
c2 = sp.symbols("c2", real=True) # pendulum damping coefficient

# Variables/Functions
t = sp.symbols("t", real=True)
x = sp.Function('x', real=True)(t)
theta = sp.Function('θ', real=True)(t)
f = sp.Function('f', real=True)(t)

# Derivatives
v = x.diff(t)
a = v.diff(t)
omega = theta.diff(t)
alpha = omega.diff(t)

# Lagrangian
V_cart = 0
T_cart = m1 * v**2 / 2
V_bob = -m2 * g * length * sp.cos(theta)
T_bob = m2 * v**2 / 2 + m2 * v * omega * length * sp.cos(theta) + m2 * omega**2 * length**2 / 2
L = T_cart + T_bob - V_cart - V_bob

# Euler-Lagrange Equations
cart_force_term = f
cart_damping_term = -c1 * v
pendulum_damping_term = -c2 * omega
x_equation = L.diff(v, t) - L.diff(x) - cart_force_term - cart_damping_term
theta_equation = L.diff(omega, t) - L.diff(theta) - pendulum_damping_term

# Solve
solution = sp.solve([x_equation, theta_equation], [a, alpha])
a_solution = solution[a]
alpha_solution = solution[alpha]
print("cart acceleration (x double dot)")
display(sp.simplify(a_solution))
print("pendulum angular acceleration (theta double dot)")
display(sp.simplify(alpha_solution))

cart acceleration (x double dot)

$\displaystyle \frac{- c_{1} l \frac{d}{d t} x{\left(t \right)} + c_{2} \cos{\left(θ{\left(t \right)} \right)} \frac{d}{d t} θ{\left(t \right)} + \frac{g l m_{2} \sin{\left(2 θ{\left(t \right)} \right)}}{2} + l^{2} m_{2} \sin{\left(θ{\left(t \right)} \right)} \left(\frac{d}{d t} θ{\left(t \right)}\right)^{2} + l f{\left(t \right)}}{l \left(m_{1} + m_{2} \sin^{2}{\left(θ{\left(t \right)} \right)}\right)}$

pendulum angular acceleration (theta double dot)

$\displaystyle \frac{c_{1} l m_{2} \cos{\left(θ{\left(t \right)} \right)} \frac{d}{d t} x{\left(t \right)} - c_{2} m_{1} \frac{d}{d t} θ{\left(t \right)} - c_{2} m_{2} \frac{d}{d t} θ{\left(t \right)} - g l m_{1} m_{2} \sin{\left(θ{\left(t \right)} \right)} - g l m_{2}^{2} \sin{\left(θ{\left(t \right)} \right)} - \frac{l^{2} m_{2}^{2} \sin{\left(2 θ{\left(t \right)} \right)} \left(\frac{d}{d t} θ{\left(t \right)}\right)^{2}}{2} - l m_{2} f{\left(t \right)} \cos{\left(θ{\left(t \right)} \right)}}{l^{2} m_{2} \left(m_{1} + m_{2} \sin^{2}{\left(θ{\left(t \right)} \right)}\right)}$

Let's code this up.

# Constants
m1 = 1.0      # mass of cart (kg)
m2 = 1.0      # mass of pendulum bob (kg)
length = 1.0  # length of pendulum (m)
g = 9.81      # acceleration due to gravity (m/s^2)
c1 = 0.02     # damping coefficient of cart (N/(m/s))
c2 = 0.02     # damping coefficient of pendulum (N/(rad/s))


def derivatives(state, f):
    """
    Given a state vector y = (x, v, theta, omega) and the external force (f), this method computes
    the derivative dy/dt = (dx/dt, dv/dt, d theta/dt, d omega/dt).

    This lets us express the dynamics of the system as a first-order ODE: dy/dt = derivatives(y, f).
    """
    v = state[1]
    theta = state[2]
    omega = state[3]
    ct = np.cos(theta)
    st = np.sin(theta)

    denominator = m1 + m2 * st**2  # both equations include this term

    x_acceleration = m2 * (g * ct + length * omega**2) * st  # Euler-Lagrange (conservative forces)
    x_acceleration += f - c1 * v + c2 * omega * ct / length  # external forces
    x_acceleration /= denominator

    theta_acceleration = -(m2 * length * omega**2 * ct + g * (m1 + m2)) * st  # Euler-Lagrange
    theta_acceleration += - f * ct + c1 * v * ct - c2 * (m1 + m2) * omega     # external forces
    theta_acceleration /= length * denominator

    return np.array([v, x_acceleration, omega, theta_acceleration])


def rk4_step(state, f, dt):
    """
    Computes one step of the Runge-Kutta 4th order integrator.

    The first argument is the current state of the system (x, v, theta, omega), the second is the
    external force (which we assume is constant over the integration interval), and the third
    argument is the time step (in seconds).

    It returns a new state vector (x, v, theta, omega).
    """
    k1 = derivatives(state, f)
    k2 = derivatives(state + 0.5 * dt * k1, f)
    k3 = derivatives(state + 0.5 * dt * k2, f)
    k4 = derivatives(state + dt * k3, f)
    return state + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)


def generate_trajectory(initial_conditions, dt, steps, f = None):
    """
    Integrates forward from a set of initial conditions.

    The first argument specifies (x, v, theta, omega) at time t = 0. The second argument is the
    timestep in seconds, and the third is the number of steps to integrate for. If provided, f is
    used as the external force (in Newtons) applied to the cart. It should be an array, and len(f)
    == steps. Otherwise a force of zero is assumed for all time.

    Returns an array of shape (steps + 1, 6) where each row is a state vector (t, x, v, theta,
    omega, f).
    """
    # We store states as (t, x, v, theta, omega, f). The first row holds the initial conditions, and
    # row i holds the state after i steps (i.e. at time t = i * dt).
    states = np.zeros((steps + 1, 6))
    states[0, 1:5] = initial_conditions
    states[:, 0] = np.array([step * dt for step in range(steps + 1)])
    # Apply a sinusoidal external force.
    if f is not None:
        states[0:steps, 5] = f

    # Simulate.
    for step in range(steps):
        states[step + 1, 1:5] = rk4_step(states[step, 1:5], states[step, 5], dt)

    return states


def plot_trajectory(trajectory):
    fig, axs = plt.subplots(nrows=4, ncols=1)

    axs[3].set_xlabel("time (s)")
    axs[0].set_ylabel("x (m)")
    axs[1].set_ylabel("v (m/s)")
    axs[2].set_ylabel("theta (rad)")
    axs[3].set_ylabel("omega (rad/s)")

    axs[0].plot(trajectory[:, 0], trajectory[:, 1])
    axs[1].plot(trajectory[:, 0], trajectory[:, 2])
    axs[2].plot(trajectory[:, 0], trajectory[:, 3])
    axs[3].plot(trajectory[:, 0], trajectory[:, 4])

    fig.tight_layout()


# Define simulation parameters.
dt = 1e-2     # integration time step (s)
t_max = 10.0  # simulation duration (s)
steps = int(t_max / dt)

# Let the pendulum swing freely.
x = 0.0
v = 0.0
theta = 0.5 * np.pi
omega = 0.0
trajectory = generate_trajectory(np.array([x, v, theta, omega]), dt, steps)
plot_trajectory(trajectory)

# Apply a sinusoidal force.
theta = 0.0
f = 5.0 * np.array([1.2 * np.sin(np.pi * t * dt) for t in range(steps)])
trajectory = generate_trajectory(np.array([x, v, theta, omega]), dt, steps, f)
plot_trajectory(trajectory)

Double Pendulum

We can derive the equations of motion in the same way, but they get pretty complicated. I'm going to use sympy all the way through.

from sympy.vector import CoordSys3D
import sympy as sp

# Constants
m1 = sp.symbols("m1", real=True) # mass of cart (kg)
m2 = sp.symbols("m2", real=True) # mass of first arm (kg)
m3 = sp.symbols("m3", real=True) # mass of second arm (kg)
l1 = sp.symbols("l1", real=True) # length of first arm (m)
l2 = sp.symbols("l2", real=True) # length of second arm (m)
g = sp.symbols("g", real=True)   # acceleration due to gravity (m/s^2)
c1 = sp.symbols("c1", real=True) # cart damping coefficient
c2 = sp.symbols("c2", real=True) # arm 1 damping coefficient
c3 = sp.symbols("c3", real=True) # arm 2 damping coefficient

# Variables/Functions
t = sp.symbols("t", real=True)
x = sp.Function('x', real=True)(t)
theta1 = sp.Function('theta1', real=True)(t)
theta2 = sp.Function('theta2', real=True)(t)
f = sp.Function('f', real=True)(t)

# Derivatives
v = x.diff(t)
a = v.diff(t)
omega1 = theta1.diff(t)
alpha1 = omega1.diff(t)
omega2 = theta2.diff(t)
alpha2 = omega2.diff(t)

# Positions and Velocities
R3 = CoordSys3D('R3')
x1 = x * R3.i
x2_displacement = l1 * sp.sin(theta1) * R3.i - l1 * sp.cos(theta1) * R3.j
x2 = x1 + x2_displacement / 2  # center of mass of the first arm
x3_displacement = l2 * sp.sin(theta2) * R3.i - l2 * sp.cos(theta2) * R3.j
x3 = x1 + x2_displacement + x3_displacement / 2  # center of mass of the second arm
v1 = x1.diff(t)
v2 = x2.diff(t)
v3 = x3.diff(t)
v1_mag_sqrd = sp.simplify(v1.dot(v1))
v2_mag_sqrd = sp.simplify(v2.dot(v2))
v3_mag_sqrd = sp.simplify(v3.dot(v3))

# Lagrangian
V1 = 0
V2 = m2 * g * x2.dot(R3.j)  # mgh (potential energy due to gravity)
V3 = m3 * g * x3.dot(R3.j)
# The cart only has linear kinetic energy.
T1 = m1 / 2 * v1_mag_sqrd   # linear kinetic energy
# The arms have linear and rotational kinetic energy. The moment of inertia of a rod of mass m and
# length l about its center of mass (and perpendicular to the rod) is 1/12 * m * l^2.
T2 = m2 / 2 * v2_mag_sqrd + m2 * l1**2 / 12 * omega1**2
T3 = m3 / 2 * v3_mag_sqrd + m3 * l2**2 / 12 * omega2**2
L = sp.simplify(T1 + T2 + T3 - V1 - V2 - V3)
print("Lagrangian")
display(L)

# Euler-Lagrange Equations
damping1 = -c1 * v
damping2 = -c2 * omega1
damping3 = -c3 * (omega2 - omega1)
# Sympy struggles to solve the equations with damping on the rotary joints.
# (Ran for 10+ minutes with no result.)
x_equation = sp.simplify(L.diff(v, t) - L.diff(x) - f - damping1)
theta1_equation = sp.simplify(L.diff(omega1, t) - L.diff(theta1)) # - damping2)
theta2_equation = sp.simplify(L.diff(omega2, t) - L.diff(theta2)) # - damping3)

# Solving and simplifying can take some time (2 minutes, 40 seconds on a 2023 MacBook Pro).
print("Solve...")
solution = sp.solve([x_equation, theta1_equation, theta2_equation], [a, alpha1, alpha2])
a1_sol = solution[a]
a2_sol = solution[alpha1]
a3_sol = solution[alpha2]

print("Simplify...")
a1_sol = sp.simplify(a1_sol)
a2_sol = sp.simplify(a2_sol)
a3_sol = sp.simplify(a3_sol)
# print("x double dot")
display(a1_sol)
print("theta1 double dot")
display(a2_sol)
print("theta2 double dot")
display(a3_sol)

# Replace symbols with text.
# Instead of x(t), we want x. Instead of Derivative(x(t), t), we want v, etc.
substitutions = [(v, "v"), (x, "x"), (omega1, "omega1"), (theta1, "theta1"), (omega2, "omega2"), (theta2, "theta2"), (f, "f")]
a1_sol = a1_sol.subs(substitutions)
a2_sol = a2_sol.subs(substitutions)
a3_sol = a3_sol.subs(substitutions)
print(f"a1 =", str(a1_sol).replace("cos", "np.cos").replace("sin", "np.sin"))
print(f"a2 =", str(a2_sol).replace("cos", "np.cos").replace("sin", "np.sin"))
print(f"a3 =", str(a3_sol).replace("cos", "np.cos").replace("sin", "np.sin"))

# This should break out subexpressions that are used multiple times.
# So instead of a single extremely long line for each acceleration, we'll have many short lines.
# But there seems to be a bug -- sometimes I get subexpressions that aren't used, and the physics
# is obviously wrong (the cart accelerates off to infinity, for example).
"""
print("Common subexpression elimination...")
subexpressions, expressions = sp.cse([a1_sol, a2_sol, a3_sol], optimizations='basic')
for name, value in subexpressions:
    value = str(value).replace("cos", "np.cos").replace("sin", "np.sin")
    print(f"{name} = {value}")
print(f"a1 =", str(expressions[0]).replace("cos", "np.cos").replace("sin", "np.sin"))
print(f"a2 =", str(expressions[0]).replace("cos", "np.cos").replace("sin", "np.sin"))
print(f"a3 =", str(expressions[1]).replace("cos", "np.cos").replace("sin", "np.sin"))
"""

Lagrangian

$\displaystyle \frac{g l_{1} m_{2} \cos{\left(\theta_{1}{\left(t \right)} \right)}}{2} + g l_{1} m_{3} \cos{\left(\theta_{1}{\left(t \right)} \right)} + \frac{g l_{2} m_{3} \cos{\left(\theta_{2}{\left(t \right)} \right)}}{2} + \frac{5 l_{1}^{2} m_{2} \left(\frac{d}{d t} \theta_{1}{\left(t \right)}\right)^{2}}{24} + \frac{l_{1}^{2} m_{3} \left(\frac{d}{d t} \theta_{1}{\left(t \right)}\right)^{2}}{2} + \frac{l_{1} l_{2} m_{3} \cos{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \frac{d}{d t} \theta_{1}{\left(t \right)} \frac{d}{d t} \theta_{2}{\left(t \right)}}{2} + \frac{l_{1} m_{2} \cos{\left(\theta_{1}{\left(t \right)} \right)} \frac{d}{d t} \theta_{1}{\left(t \right)} \frac{d}{d t} x{\left(t \right)}}{2} + l_{1} m_{3} \cos{\left(\theta_{1}{\left(t \right)} \right)} \frac{d}{d t} \theta_{1}{\left(t \right)} \frac{d}{d t} x{\left(t \right)} + \frac{5 l_{2}^{2} m_{3} \left(\frac{d}{d t} \theta_{2}{\left(t \right)}\right)^{2}}{24} + \frac{l_{2} m_{3} \cos{\left(\theta_{2}{\left(t \right)} \right)} \frac{d}{d t} \theta_{2}{\left(t \right)} \frac{d}{d t} x{\left(t \right)}}{2} + \frac{m_{1} \left(\frac{d}{d t} x{\left(t \right)}\right)^{2}}{2} + \frac{m_{2} \left(\frac{d}{d t} x{\left(t \right)}\right)^{2}}{2} + \frac{m_{3} \left(\frac{d}{d t} x{\left(t \right)}\right)^{2}}{2}$

Solve...
Simplify...

$\displaystyle \frac{50 c_{1} m_{2} \frac{d}{d t} x{\left(t \right)} - 72 c_{1} m_{3} \cos^{2}{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \frac{d}{d t} x{\left(t \right)} + 120 c_{1} m_{3} \frac{d}{d t} x{\left(t \right)} - 15 g m_{2}^{2} \sin{\left(2 \theta_{1}{\left(t \right)} \right)} + 36 g m_{2} m_{3} \sin{\left(\theta_{1}{\left(t \right)} \right)} \cos{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \cos{\left(\theta_{2}{\left(t \right)} \right)} - 60 g m_{2} m_{3} \sin{\left(2 \theta_{1}{\left(t \right)} \right)} + 36 g m_{2} m_{3} \sin{\left(\theta_{2}{\left(t \right)} \right)} \cos{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \cos{\left(\theta_{1}{\left(t \right)} \right)} - 15 g m_{2} m_{3} \sin{\left(2 \theta_{2}{\left(t \right)} \right)} + 72 g m_{3}^{2} \sin{\left(\theta_{1}{\left(t \right)} \right)} \cos{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \cos{\left(\theta_{2}{\left(t \right)} \right)} - 60 g m_{3}^{2} \sin{\left(2 \theta_{1}{\left(t \right)} \right)} + 72 g m_{3}^{2} \sin{\left(\theta_{2}{\left(t \right)} \right)} \cos{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \cos{\left(\theta_{1}{\left(t \right)} \right)} - 36 g m_{3}^{2} \sin{\left(2 \theta_{2}{\left(t \right)} \right)} - 25 l_{1} m_{2}^{2} \sin{\left(\theta_{1}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{1}{\left(t \right)}\right)^{2} - 9 l_{1} m_{2} m_{3} \left(\sin{\left(\theta_{1}{\left(t \right)} - 2 \theta_{2}{\left(t \right)} \right)} + \sin{\left(3 \theta_{1}{\left(t \right)} - 2 \theta_{2}{\left(t \right)} \right)}\right) \left(\frac{d}{d t} \theta_{1}{\left(t \right)}\right)^{2} + 30 l_{1} m_{2} m_{3} \sin{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \cos{\left(\theta_{2}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{1}{\left(t \right)}\right)^{2} + 36 l_{1} m_{2} m_{3} \sin{\left(\theta_{1}{\left(t \right)} \right)} \cos^{2}{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{1}{\left(t \right)}\right)^{2} - 110 l_{1} m_{2} m_{3} \sin{\left(\theta_{1}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{1}{\left(t \right)}\right)^{2} - 18 l_{1} m_{3}^{2} \left(\sin{\left(\theta_{1}{\left(t \right)} - 2 \theta_{2}{\left(t \right)} \right)} + \sin{\left(3 \theta_{1}{\left(t \right)} - 2 \theta_{2}{\left(t \right)} \right)}\right) \left(\frac{d}{d t} \theta_{1}{\left(t \right)}\right)^{2} + 72 l_{1} m_{3}^{2} \sin{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \cos{\left(\theta_{2}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{1}{\left(t \right)}\right)^{2} + 72 l_{1} m_{3}^{2} \sin{\left(\theta_{1}{\left(t \right)} \right)} \cos^{2}{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{1}{\left(t \right)}\right)^{2} - 120 l_{1} m_{3}^{2} \sin{\left(\theta_{1}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{1}{\left(t \right)}\right)^{2} - 30 l_{2} m_{2} m_{3} \sin{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \cos{\left(\theta_{1}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{2}{\left(t \right)}\right)^{2} - 25 l_{2} m_{2} m_{3} \sin{\left(\theta_{2}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{2}{\left(t \right)}\right)^{2} + 9 l_{2} m_{3}^{2} \left(\sin{\left(2 \theta_{1}{\left(t \right)} - 3 \theta_{2}{\left(t \right)} \right)} + \sin{\left(2 \theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)}\right) \left(\frac{d}{d t} \theta_{2}{\left(t \right)}\right)^{2} - 60 l_{2} m_{3}^{2} \sin{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \cos{\left(\theta_{1}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{2}{\left(t \right)}\right)^{2} + 36 l_{2} m_{3}^{2} \sin{\left(\theta_{2}{\left(t \right)} \right)} \cos^{2}{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{2}{\left(t \right)}\right)^{2} - 60 l_{2} m_{3}^{2} \sin{\left(\theta_{2}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{2}{\left(t \right)}\right)^{2} - 50 m_{2} f{\left(t \right)} + 72 m_{3} f{\left(t \right)} \cos^{2}{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} - 120 m_{3} f{\left(t \right)}}{2 \left(- 25 m_{1} m_{2} + 36 m_{1} m_{3} \cos^{2}{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} - 60 m_{1} m_{3} + 15 m_{2}^{2} \cos^{2}{\left(\theta_{1}{\left(t \right)} \right)} - 25 m_{2}^{2} + 36 m_{2} m_{3} \cos^{2}{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} - 36 m_{2} m_{3} \cos{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \cos{\left(\theta_{1}{\left(t \right)} \right)} \cos{\left(\theta_{2}{\left(t \right)} \right)} + 60 m_{2} m_{3} \cos^{2}{\left(\theta_{1}{\left(t \right)} \right)} + 15 m_{2} m_{3} \cos^{2}{\left(\theta_{2}{\left(t \right)} \right)} - 85 m_{2} m_{3} + 36 m_{3}^{2} \cos^{2}{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} - 72 m_{3}^{2} \cos{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \cos{\left(\theta_{1}{\left(t \right)} \right)} \cos{\left(\theta_{2}{\left(t \right)} \right)} + 60 m_{3}^{2} \cos^{2}{\left(\theta_{1}{\left(t \right)} \right)} + 36 m_{3}^{2} \cos^{2}{\left(\theta_{2}{\left(t \right)} \right)} - 60 m_{3}^{2}\right)}$

theta1 double dot

$\displaystyle \frac{- 30 c_{1} m_{2} \cos{\left(\theta_{1}{\left(t \right)} \right)} \frac{d}{d t} x{\left(t \right)} + 18 c_{1} m_{3} \cos{\left(\theta_{1}{\left(t \right)} - 2 \theta_{2}{\left(t \right)} \right)} \frac{d}{d t} x{\left(t \right)} - 42 c_{1} m_{3} \cos{\left(\theta_{1}{\left(t \right)} \right)} \frac{d}{d t} x{\left(t \right)} + 30 g m_{1} m_{2} \sin{\left(\theta_{1}{\left(t \right)} \right)} + 18 g m_{1} m_{3} \sin{\left(\theta_{1}{\left(t \right)} - 2 \theta_{2}{\left(t \right)} \right)} + 42 g m_{1} m_{3} \sin{\left(\theta_{1}{\left(t \right)} \right)} + 30 g m_{2}^{2} \sin{\left(\theta_{1}{\left(t \right)} \right)} + 9 g m_{2} m_{3} \sin{\left(\theta_{1}{\left(t \right)} - 2 \theta_{2}{\left(t \right)} \right)} + 63 g m_{2} m_{3} \sin{\left(\theta_{1}{\left(t \right)} \right)} + 24 g m_{3}^{2} \sin{\left(\theta_{1}{\left(t \right)} \right)} + 18 l_{1} m_{1} m_{3} \sin{\left(2 \theta_{1}{\left(t \right)} - 2 \theta_{2}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{1}{\left(t \right)}\right)^{2} + \frac{15 l_{1} m_{2}^{2} \sin{\left(2 \theta_{1}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{1}{\left(t \right)}\right)^{2}}{2} + 9 l_{1} m_{2} m_{3} \sin{\left(2 \theta_{1}{\left(t \right)} - 2 \theta_{2}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{1}{\left(t \right)}\right)^{2} + 21 l_{1} m_{2} m_{3} \sin{\left(2 \theta_{1}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{1}{\left(t \right)}\right)^{2} + 12 l_{1} m_{3}^{2} \sin{\left(2 \theta_{1}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{1}{\left(t \right)}\right)^{2} + 30 l_{2} m_{1} m_{3} \sin{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{2}{\left(t \right)}\right)^{2} + \frac{45 l_{2} m_{2} m_{3} \sin{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{2}{\left(t \right)}\right)^{2}}{2} + \frac{15 l_{2} m_{2} m_{3} \sin{\left(\theta_{1}{\left(t \right)} + \theta_{2}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{2}{\left(t \right)}\right)^{2}}{2} + 6 l_{2} m_{3}^{2} \sin{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{2}{\left(t \right)}\right)^{2} + 6 l_{2} m_{3}^{2} \sin{\left(\theta_{1}{\left(t \right)} + \theta_{2}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{2}{\left(t \right)}\right)^{2} + 30 m_{2} f{\left(t \right)} \cos{\left(\theta_{1}{\left(t \right)} \right)} - 18 m_{3} f{\left(t \right)} \cos{\left(\theta_{1}{\left(t \right)} - 2 \theta_{2}{\left(t \right)} \right)} + 42 m_{3} f{\left(t \right)} \cos{\left(\theta_{1}{\left(t \right)} \right)}}{l_{1} \left(- 25 m_{1} m_{2} + 36 m_{1} m_{3} \cos^{2}{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} - 60 m_{1} m_{3} + 15 m_{2}^{2} \cos^{2}{\left(\theta_{1}{\left(t \right)} \right)} - 25 m_{2}^{2} + 36 m_{2} m_{3} \cos^{2}{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} - 36 m_{2} m_{3} \cos{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \cos{\left(\theta_{1}{\left(t \right)} \right)} \cos{\left(\theta_{2}{\left(t \right)} \right)} + 60 m_{2} m_{3} \cos^{2}{\left(\theta_{1}{\left(t \right)} \right)} + 15 m_{2} m_{3} \cos^{2}{\left(\theta_{2}{\left(t \right)} \right)} - 85 m_{2} m_{3} + 36 m_{3}^{2} \cos^{2}{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} - 72 m_{3}^{2} \cos{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \cos{\left(\theta_{1}{\left(t \right)} \right)} \cos{\left(\theta_{2}{\left(t \right)} \right)} + 60 m_{3}^{2} \cos^{2}{\left(\theta_{1}{\left(t \right)} \right)} + 36 m_{3}^{2} \cos^{2}{\left(\theta_{2}{\left(t \right)} \right)} - 60 m_{3}^{2}\right)}$

theta2 double dot

$\displaystyle \frac{18 c_{1} m_{2} \cos{\left(2 \theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \frac{d}{d t} x{\left(t \right)} - 12 c_{1} m_{2} \cos{\left(\theta_{2}{\left(t \right)} \right)} \frac{d}{d t} x{\left(t \right)} + 36 c_{1} m_{3} \cos{\left(2 \theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \frac{d}{d t} x{\left(t \right)} - 36 c_{1} m_{3} \cos{\left(\theta_{2}{\left(t \right)} \right)} \frac{d}{d t} x{\left(t \right)} - 18 g m_{1} m_{2} \sin{\left(2 \theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} + 12 g m_{1} m_{2} \sin{\left(\theta_{2}{\left(t \right)} \right)} - 36 g m_{1} m_{3} \sin{\left(2 \theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} + 36 g m_{1} m_{3} \sin{\left(\theta_{2}{\left(t \right)} \right)} - 9 g m_{2}^{2} \sin{\left(2 \theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} + 3 g m_{2}^{2} \sin{\left(\theta_{2}{\left(t \right)} \right)} - 18 g m_{2} m_{3} \sin{\left(2 \theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} + 12 g m_{2} m_{3} \sin{\left(\theta_{2}{\left(t \right)} \right)} - 30 l_{1} m_{1} m_{2} \sin{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{1}{\left(t \right)}\right)^{2} - 72 l_{1} m_{1} m_{3} \sin{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{1}{\left(t \right)}\right)^{2} - \frac{27 l_{1} m_{2}^{2} \sin{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{1}{\left(t \right)}\right)^{2}}{2} - \frac{3 l_{1} m_{2}^{2} \sin{\left(\theta_{1}{\left(t \right)} + \theta_{2}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{1}{\left(t \right)}\right)^{2}}{2} - 33 l_{1} m_{2} m_{3} \sin{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{1}{\left(t \right)}\right)^{2} - 3 l_{1} m_{2} m_{3} \sin{\left(\theta_{1}{\left(t \right)} + \theta_{2}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{1}{\left(t \right)}\right)^{2} - 18 l_{2} m_{1} m_{3} \sin{\left(2 \theta_{1}{\left(t \right)} - 2 \theta_{2}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{2}{\left(t \right)}\right)^{2} - 9 l_{2} m_{2} m_{3} \sin{\left(2 \theta_{1}{\left(t \right)} - 2 \theta_{2}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{2}{\left(t \right)}\right)^{2} - \frac{3 l_{2} m_{2} m_{3} \sin{\left(2 \theta_{2}{\left(t \right)} \right)} \left(\frac{d}{d t} \theta_{2}{\left(t \right)}\right)^{2}}{2} - 18 m_{2} f{\left(t \right)} \cos{\left(2 \theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} + 12 m_{2} f{\left(t \right)} \cos{\left(\theta_{2}{\left(t \right)} \right)} - 36 m_{3} f{\left(t \right)} \cos{\left(2 \theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} + 36 m_{3} f{\left(t \right)} \cos{\left(\theta_{2}{\left(t \right)} \right)}}{l_{2} \left(- 25 m_{1} m_{2} + 36 m_{1} m_{3} \cos^{2}{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} - 60 m_{1} m_{3} + 15 m_{2}^{2} \cos^{2}{\left(\theta_{1}{\left(t \right)} \right)} - 25 m_{2}^{2} + 36 m_{2} m_{3} \cos^{2}{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} - 36 m_{2} m_{3} \cos{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \cos{\left(\theta_{1}{\left(t \right)} \right)} \cos{\left(\theta_{2}{\left(t \right)} \right)} + 60 m_{2} m_{3} \cos^{2}{\left(\theta_{1}{\left(t \right)} \right)} + 15 m_{2} m_{3} \cos^{2}{\left(\theta_{2}{\left(t \right)} \right)} - 85 m_{2} m_{3} + 36 m_{3}^{2} \cos^{2}{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} - 72 m_{3}^{2} \cos{\left(\theta_{1}{\left(t \right)} - \theta_{2}{\left(t \right)} \right)} \cos{\left(\theta_{1}{\left(t \right)} \right)} \cos{\left(\theta_{2}{\left(t \right)} \right)} + 60 m_{3}^{2} \cos^{2}{\left(\theta_{1}{\left(t \right)} \right)} + 36 m_{3}^{2} \cos^{2}{\left(\theta_{2}{\left(t \right)} \right)} - 60 m_{3}^{2}\right)}$

a1 = (50*c1*m2*v - 72*c1*m3*v*np.cos(theta1 - theta2)**2 + 120*c1*m3*v - 50*f*m2 + 72*f*m3*np.cos(theta1 - theta2)**2 - 120*f*m3 - 15*g*m2**2*np.sin(2*theta1) + 36*g*m2*m3*np.sin(theta1)*np.cos(theta2)*np.cos(theta1 - theta2) - 60*g*m2*m3*np.sin(2*theta1) + 36*g*m2*m3*np.sin(theta2)*np.cos(theta1)*np.cos(theta1 - theta2) - 15*g*m2*m3*np.sin(2*theta2) + 72*g*m3**2*np.sin(theta1)*np.cos(theta2)*np.cos(theta1 - theta2) - 60*g*m3**2*np.sin(2*theta1) + 72*g*m3**2*np.sin(theta2)*np.cos(theta1)*np.cos(theta1 - theta2) - 36*g*m3**2*np.sin(2*theta2) - 25*l1*m2**2*omega1**2*np.sin(theta1) - 9*l1*m2*m3*omega1**2*(np.sin(theta1 - 2*theta2) + np.sin(3*theta1 - 2*theta2)) + 36*l1*m2*m3*omega1**2*np.sin(theta1)*np.cos(theta1 - theta2)**2 - 110*l1*m2*m3*omega1**2*np.sin(theta1) + 30*l1*m2*m3*omega1**2*np.sin(theta1 - theta2)*np.cos(theta2) - 18*l1*m3**2*omega1**2*(np.sin(theta1 - 2*theta2) + np.sin(3*theta1 - 2*theta2)) + 72*l1*m3**2*omega1**2*np.sin(theta1)*np.cos(theta1 - theta2)**2 - 120*l1*m3**2*omega1**2*np.sin(theta1) + 72*l1*m3**2*omega1**2*np.sin(theta1 - theta2)*np.cos(theta2) - 25*l2*m2*m3*omega2**2*np.sin(theta2) - 30*l2*m2*m3*omega2**2*np.sin(theta1 - theta2)*np.cos(theta1) + 9*l2*m3**2*omega2**2*(np.sin(2*theta1 - 3*theta2) + np.sin(2*theta1 - theta2)) + 36*l2*m3**2*omega2**2*np.sin(theta2)*np.cos(theta1 - theta2)**2 - 60*l2*m3**2*omega2**2*np.sin(theta2) - 60*l2*m3**2*omega2**2*np.sin(theta1 - theta2)*np.cos(theta1))/(2*(-25*m1*m2 + 36*m1*m3*np.cos(theta1 - theta2)**2 - 60*m1*m3 + 15*m2**2*np.cos(theta1)**2 - 25*m2**2 + 60*m2*m3*np.cos(theta1)**2 - 36*m2*m3*np.cos(theta1)*np.cos(theta2)*np.cos(theta1 - theta2) + 15*m2*m3*np.cos(theta2)**2 + 36*m2*m3*np.cos(theta1 - theta2)**2 - 85*m2*m3 + 60*m3**2*np.cos(theta1)**2 - 72*m3**2*np.cos(theta1)*np.cos(theta2)*np.cos(theta1 - theta2) + 36*m3**2*np.cos(theta2)**2 + 36*m3**2*np.cos(theta1 - theta2)**2 - 60*m3**2))
a2 = (-30*c1*m2*v*np.cos(theta1) - 42*c1*m3*v*np.cos(theta1) + 18*c1*m3*v*np.cos(theta1 - 2*theta2) + 30*f*m2*np.cos(theta1) + 42*f*m3*np.cos(theta1) - 18*f*m3*np.cos(theta1 - 2*theta2) + 30*g*m1*m2*np.sin(theta1) + 42*g*m1*m3*np.sin(theta1) + 18*g*m1*m3*np.sin(theta1 - 2*theta2) + 30*g*m2**2*np.sin(theta1) + 63*g*m2*m3*np.sin(theta1) + 9*g*m2*m3*np.sin(theta1 - 2*theta2) + 24*g*m3**2*np.sin(theta1) + 18*l1*m1*m3*omega1**2*np.sin(2*theta1 - 2*theta2) + 15*l1*m2**2*omega1**2*np.sin(2*theta1)/2 + 21*l1*m2*m3*omega1**2*np.sin(2*theta1) + 9*l1*m2*m3*omega1**2*np.sin(2*theta1 - 2*theta2) + 12*l1*m3**2*omega1**2*np.sin(2*theta1) + 30*l2*m1*m3*omega2**2*np.sin(theta1 - theta2) + 45*l2*m2*m3*omega2**2*np.sin(theta1 - theta2)/2 + 15*l2*m2*m3*omega2**2*np.sin(theta1 + theta2)/2 + 6*l2*m3**2*omega2**2*np.sin(theta1 - theta2) + 6*l2*m3**2*omega2**2*np.sin(theta1 + theta2))/(l1*(-25*m1*m2 + 36*m1*m3*np.cos(theta1 - theta2)**2 - 60*m1*m3 + 15*m2**2*np.cos(theta1)**2 - 25*m2**2 + 60*m2*m3*np.cos(theta1)**2 - 36*m2*m3*np.cos(theta1)*np.cos(theta2)*np.cos(theta1 - theta2) + 15*m2*m3*np.cos(theta2)**2 + 36*m2*m3*np.cos(theta1 - theta2)**2 - 85*m2*m3 + 60*m3**2*np.cos(theta1)**2 - 72*m3**2*np.cos(theta1)*np.cos(theta2)*np.cos(theta1 - theta2) + 36*m3**2*np.cos(theta2)**2 + 36*m3**2*np.cos(theta1 - theta2)**2 - 60*m3**2))
a3 = (-12*c1*m2*v*np.cos(theta2) + 18*c1*m2*v*np.cos(2*theta1 - theta2) - 36*c1*m3*v*np.cos(theta2) + 36*c1*m3*v*np.cos(2*theta1 - theta2) + 12*f*m2*np.cos(theta2) - 18*f*m2*np.cos(2*theta1 - theta2) + 36*f*m3*np.cos(theta2) - 36*f*m3*np.cos(2*theta1 - theta2) + 12*g*m1*m2*np.sin(theta2) - 18*g*m1*m2*np.sin(2*theta1 - theta2) + 36*g*m1*m3*np.sin(theta2) - 36*g*m1*m3*np.sin(2*theta1 - theta2) + 3*g*m2**2*np.sin(theta2) - 9*g*m2**2*np.sin(2*theta1 - theta2) + 12*g*m2*m3*np.sin(theta2) - 18*g*m2*m3*np.sin(2*theta1 - theta2) - 30*l1*m1*m2*omega1**2*np.sin(theta1 - theta2) - 72*l1*m1*m3*omega1**2*np.sin(theta1 - theta2) - 27*l1*m2**2*omega1**2*np.sin(theta1 - theta2)/2 - 3*l1*m2**2*omega1**2*np.sin(theta1 + theta2)/2 - 33*l1*m2*m3*omega1**2*np.sin(theta1 - theta2) - 3*l1*m2*m3*omega1**2*np.sin(theta1 + theta2) - 18*l2*m1*m3*omega2**2*np.sin(2*theta1 - 2*theta2) - 3*l2*m2*m3*omega2**2*np.sin(2*theta2)/2 - 9*l2*m2*m3*omega2**2*np.sin(2*theta1 - 2*theta2))/(l2*(-25*m1*m2 + 36*m1*m3*np.cos(theta1 - theta2)**2 - 60*m1*m3 + 15*m2**2*np.cos(theta1)**2 - 25*m2**2 + 60*m2*m3*np.cos(theta1)**2 - 36*m2*m3*np.cos(theta1)*np.cos(theta2)*np.cos(theta1 - theta2) + 15*m2*m3*np.cos(theta2)**2 + 36*m2*m3*np.cos(theta1 - theta2)**2 - 85*m2*m3 + 60*m3**2*np.cos(theta1)**2 - 72*m3**2*np.cos(theta1)*np.cos(theta2)*np.cos(theta1 - theta2) + 36*m3**2*np.cos(theta2)**2 + 36*m3**2*np.cos(theta1 - theta2)**2 - 60*m3**2))

'\nprint("Common subexpression elimination...")\nsubexpressions, expressions = sp.cse([a1_sol, a2_sol, a3_sol], optimizations=\'basic\')\nfor name, value in subexpressions:\n    value = str(value).replace("cos", "np.cos").replace("sin", "np.sin")\n    print(f"{name} = {value}")\nprint(f"a1 =", str(expressions[0]).replace("cos", "np.cos").replace("sin", "np.sin"))\nprint(f"a2 =", str(expressions[0]).replace("cos", "np.cos").replace("sin", "np.sin"))\nprint(f"a3 =", str(expressions[1]).replace("cos", "np.cos").replace("sin", "np.sin"))\n'