In Chapter 3, translational systems were modeled using Newton’s law. In Lesson 1 of this chapter, rotational systems were modeled using the rotational analogue. In Lesson 2, we introduced ideal transformers (e.g., gears and levers) to map effort/flow variables between domains.

A coupled translational–rotational mechanism adds an additional ingredient: a kinematic constraint that enforces a deterministic relationship between a translational coordinate \( x(t) \) and a rotational coordinate \( \theta(t) \). In its most general form:

\[ x = f(\theta), \qquad \dot{x} = \frac{d f}{d\theta}\dot{\theta}, \qquad \ddot{x} = \frac{d^2 f}{d\theta^2}\dot{\theta}^2 + \frac{d f}{d\theta}\ddot{\theta}. \]

For ideal, lossless couplings (no slip, no backlash, no compliance inside the mechanism), the correct force–torque mapping is obtained from power invariance: translational power equals rotational power.

\[ P = F \dot{x} = (\tau)\dot{\theta} \quad \Longrightarrow \quad \tau = F \frac{\dot{x}}{\dot{\theta}} = F \frac{d f}{d\theta}. \]

The scalar \( J(\theta) := \frac{d f}{d\theta} \) is the (configuration-dependent) kinematic Jacobian that converts angular motion to linear motion. When \( f(\theta) \) is linear (rack-and-pinion), \( J(\theta) \) is constant. When \( f(\theta) \) is nonlinear (crank–slider), \( J(\theta) \) varies with angle and changes the effective inertia seen by the actuator.

flowchart TD
  A["Choose coordinates (x, theta) and sign conventions"] --> B["Write kinematic constraint: x = f(theta)"]
  B --> C["Differentiate: xdot = (dx/dtheta)*thetadot ; \nxddot = (d2x/dtheta2)*thetadot^2 + (dx/dtheta)*thetaddot"]
  C --> D["Power consistency: tau*thetadot = F*xdot  =>  tau = F*(dx/dtheta)"]
  D --> E["Reflect elements across coupling (mass, spring, damper)"]
  E --> F["Assemble ODE in one coordinate (often theta)"]
  F --> G["Sanity checks: units, energy/power sign, limiting cases"]
  G --> H["Simulate (ODE solver) and validate against intuition"]
        

Assumptions used throughout this lesson: rigid links, ideal joints, no gear backlash, no slip, and (unless stated) viscous damping. Non-ideal gear train effects are deferred to Lesson 4.

Consider a pinion of pitch radius \( r_p \) driving a rack with displacement \( x \). With no slip at the pitch line:

\[ x = r_p \theta,\qquad \dot{x} = r_p \dot{\theta},\qquad \ddot{x} = r_p \ddot{\theta}. \]

The mapping is linear, so the Jacobian is constant: \( J = \frac{dx}{d\theta} = r_p \). Power invariance gives:

\[ (\tau)\dot{\theta} = F \dot{x} = F (r_p \dot{\theta}) \quad \Longrightarrow \quad \tau = r_p F. \]

This relation is equivalent to a constant-ratio transformer between rotational and translational domains. It also yields a systematic method to “reflect” translational elements into rotational equivalents.

2.1 Reflected mass (equivalent inertia)

Suppose the rack drives a translational mass \( m \). If we write dynamics in the rotational coordinate \( \theta \), the translational kinetic energy is:

\[ T_m = \frac{1}{2} m \dot{x}^2 = \frac{1}{2} m (r_p^2 \dot{\theta}^2) = \frac{1}{2} (m r_p^2)\dot{\theta}^2. \]

This shows the mass behaves like an additional rotational inertia \( J_{\text{ref}} = m r_p^2 \) as seen at the pinion shaft. If the pinion has physical inertia \( J_p \), then the total effective inertia is:

\[ J_{\text{eq}} = J_p + m r_p^2. \]

Proof (power/energy consistency): The generalized momentum associated with the shaft is \( \frac{\partial T}{\partial \dot{\theta}} \). With \( T = \frac{1}{2}J_p \dot{\theta}^2 + \frac{1}{2}m r_p^2 \dot{\theta}^2 \), we get \( \frac{\partial T}{\partial \dot{\theta}} = (J_p + m r_p^2)\dot{\theta} \). For an ideal lossless mechanism, the net shaft torque equals the rate of change of this momentum, hence \( \tau_{\text{net}} = (J_p + m r_p^2)\ddot{\theta} \) when no other effects are present.

2.2 Reflected spring and damper

Suppose the rack is connected to a spring \( k \) and viscous damper \( c \) to ground:

\[ F_s = k x,\qquad F_d = c \dot{x}. \]

Converting each to torque using \( \tau = r_p F \) and \( x = r_p \theta \):

\[ \tau_s = r_p (k x) = r_p k (r_p \theta) = (k r_p^2)\theta, \qquad \tau_d = r_p (c \dot{x}) = r_p c (r_p \dot{\theta}) = (c r_p^2)\dot{\theta}. \]

Hence, the translational spring/damper reflect as torsional spring/damper: \( k_t = k r_p^2 \), \( c_t = c r_p^2 \).

2.3 Example ODE (motor torque driving rack load)

Let input torque be \( \tau_{\text{in}}(t) \), and assume viscous shaft friction \( b_\theta \dot{\theta} \). The rotational-coordinate ODE becomes:

\[ (J_p + m r_p^2)\ddot{\theta} + (b_\theta + c r_p^2)\dot{\theta} + (k r_p^2)\theta = \tau_{\text{in}}(t). \]

The corresponding rack displacement is recovered by \( x(t) = r_p \theta(t) \).

Consider a planar crank–slider with crank radius \( r \), connecting rod length \( l \), crank angle \( \theta \), and slider displacement \( x \) along the line of stroke. Using basic geometry (right triangle from the rod projection), the slider position measured from the crank center along the stroke is:

\[ x(\theta) = r \cos\theta + \sqrt{l^2 - r^2 \sin^2\theta}. \]

Differentiating yields the kinematic Jacobian (mechanism “gain”):

\[ \frac{dx}{d\theta} = -r\sin\theta - \frac{r^2 \sin\theta\cos\theta}{\sqrt{l^2 - r^2 \sin^2\theta}}. \]

The velocity mapping is then \( \dot{x} = \left(\frac{dx}{d\theta}\right)\dot{\theta} \). Power invariance provides the instantaneous torque required to balance a slider force \(F\):

\[ (\tau)\dot{\theta} = F \dot{x} = F\left(\frac{dx}{d\theta}\right)\dot{\theta} \quad \Longrightarrow \quad \tau = F \frac{dx}{d\theta}. \]

Interpretation (mechanical advantage): When \( \left|\frac{dx}{d\theta}\right| \) is small (near dead-center positions), a large torque may be required to support a given slider force. This is a purely geometric effect and will appear in dynamic simulations even for simple viscous loads.

flowchart LR
  subgraph R["Rack-and-Pinion (ideal)"]
    R1["Constraint: x = r_p*theta"] --> R2["Jacobian: dx/dtheta = r_p (constant)"]
    R2 --> R3["Power: tau*thetadot = F*xdot"]
    R3 --> R4["Mapping: tau = r_p*F ; xdot = r_p*thetadot"]
    R4 --> R5["Reflections: J_ref = m*r_p^2 ; k_t = k*r_p^2 ; c_t = c*r_p^2"]
  end

  subgraph C["Crank-Slider (ideal)"]
    C1["Constraint: x = r*cos(theta) + sqrt(l^2 - r^2*sin(theta)^2)"] --> C2["Jacobian: dx/dtheta depends on theta"]
    C2 --> C3["Power: tau = F*(dx/dtheta)"]
    C3 --> C4["Effective inertia varies with theta"]
    C4 --> C5["Nonlinear ODE in theta (configuration dependent)"]
  end
        

We now derive a compact, actuator-side ODE for a crank–slider with: crank rotational inertia \( J_c \), slider mass \( m \), slider viscous load \( c \dot{x} \), and slider spring load \( k x \). The actuator applies input torque \( \tau_{\text{in}}(t) \).

4.1 Reflected inertia as a function of angle

The total kinetic energy is the sum of crank rotation and slider translation:

\[ T = \frac{1}{2}J_c \dot{\theta}^2 + \frac{1}{2}m \dot{x}^2 = \frac{1}{2}\left(J_c + m\left(\frac{dx}{d\theta}\right)^2\right)\dot{\theta}^2. \]

Define the configuration-dependent equivalent inertia: \( J_{\text{eq}}(\theta) := J_c + m\left(\frac{dx}{d\theta}\right)^2 \).

\[ T = \frac{1}{2}J_{\text{eq}}(\theta)\dot{\theta}^2. \]

4.2 Power-balance proof of the actuator-side equation

For an ideal mechanism, the rate of change of kinetic energy equals the net mechanical power injected:

\[ \frac{dT}{dt} = \tau_{\text{net}}\,\dot{\theta}. \]

Differentiate \( T = \frac{1}{2}J_{\text{eq}}(\theta)\dot{\theta}^2 \) using the chain rule:

\[ \frac{dT}{dt} = \frac{1}{2}\frac{dJ_{\text{eq}}}{d\theta}\dot{\theta}\,\dot{\theta}^2 + J_{\text{eq}}(\theta)\dot{\theta}\ddot{\theta} = \left(J_{\text{eq}}(\theta)\ddot{\theta} + \frac{1}{2}\frac{dJ_{\text{eq}}}{d\theta}\dot{\theta}^2\right)\dot{\theta}. \]

Equating to \( \tau_{\text{net}} \dot{\theta} \) and canceling \( \dot{\theta} \) (with the understanding that the identity holds by continuity at \( \dot{\theta}=0 \)) gives the fundamental actuator-side equation:

\[ J_{\text{eq}}(\theta)\ddot{\theta} + \frac{1}{2}\frac{dJ_{\text{eq}}}{d\theta}\dot{\theta}^2 = \tau_{\text{net}}. \]

Next, compute \( \tau_{\text{net}} \) by reflecting slider forces into torque:

\[ F_{\text{load}} = k x + c \dot{x} \quad \Longrightarrow \quad \tau_{\text{load}} = F_{\text{load}}\frac{dx}{d\theta}. \]

Therefore, with \( \tau_{\text{net}} = \tau_{\text{in}} - \tau_{\text{load}} - b_\theta \dot{\theta} \), we obtain a closed ODE in \( \theta \):

\[ J_{\text{eq}}(\theta)\ddot{\theta} + \frac{1}{2}\frac{dJ_{\text{eq}}}{d\theta}\dot{\theta}^2 + b_\theta \dot{\theta} + \left(k x(\theta) + c \dot{x}\right)\frac{dx}{d\theta} = \tau_{\text{in}}(t), \]

with \( x(\theta) \) and \( \dot{x} = \left(\frac{dx}{d\theta}\right)\dot{\theta} \) given by Section 3. This is a nonlinear, configuration-dependent dynamic model—derived without introducing state-space or transfer functions (reserved for later chapters).

The goal of these implementations is consistent: define the ODE in a minimal coordinate (here \( \theta \)), then simulate with a numerical integrator. We provide both library-based solvers and at least one from-scratch integrator (RK4) to expose numerical structure.

5.1 Python (NumPy/SciPy + optional SymPy)

Recommended libraries: numpy, scipy.integrate, matplotlib, and optionally sympy for symbolic differentiation of \( x(\theta) \) and \( dx/d\theta \).

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# ---------- Mechanism: rack-and-pinion ----------
def rack_pinion_rhs(t, z, Jp, m, rp, bth, c, k, tau_in_fun):
    # z = [theta, omega]
    theta, omega = z
    Jeq = Jp + m * (rp ** 2)
    tau_in = tau_in_fun(t)
    # torque from reflected spring/damper: (k rp^2)*theta + (c rp^2)*omega
    tau_load = (k * rp ** 2) * theta + (c * rp ** 2) * omega
    domega = (tau_in - bth * omega - tau_load) / Jeq
    return np.array([omega, domega])

# ---------- Mechanism: crank-slider ----------
def x_theta(theta, r, l):
    return r * np.cos(theta) + np.sqrt(l ** 2 - (r ** 2) * (np.sin(theta) ** 2))

def dx_dtheta(theta, r, l):
    denom = np.sqrt(l ** 2 - (r ** 2) * (np.sin(theta) ** 2))
    return -r * np.sin(theta) - (r ** 2) * np.sin(theta) * np.cos(theta) / denom

def dJeq_dtheta(theta, m, r, l):
    # J_eq(theta) = Jc + m*(dx/dtheta)^2  => dJ_eq/dtheta = 2*m*(dx/dtheta)*(d2x/dtheta2)
    # We compute d2x/dtheta2 numerically for simplicity (students can replace with analytic/SymPy).
    h = 1e-6
    d1 = dx_dtheta(theta + h, r, l)
    d0 = dx_dtheta(theta - h, r, l)
    d2x = (d1 - d0) / (2.0 * h)
    return 2.0 * m * dx_dtheta(theta, r, l) * d2x

def crank_slider_rhs(t, z, Jc, m, r, l, bth, c, k, tau_in_fun):
    theta, omega = z
    Jmap = dx_dtheta(theta, r, l)
    x = x_theta(theta, r, l)
    xdot = Jmap * omega

    Jeq = Jc + m * (Jmap ** 2)
    dJeq = dJeq_dtheta(theta, m, r, l)

    tau_in = tau_in_fun(t)
    F_load = k * x + c * xdot
    tau_load = F_load * Jmap

    # J_eq*thetaddot + 0.5*(dJ_eq/dtheta)*omega^2 + bth*omega + tau_load = tau_in
    domega = (tau_in - (0.5 * dJeq * (omega ** 2)) - bth * omega - tau_load) / Jeq
    return np.array([omega, domega])

# ---------- Example simulations ----------
def step_torque(t, tau0=1.0, t0=0.05):
    return tau0 if t >= t0 else 0.0

# Rack-and-pinion parameters
Jp, m, rp = 2e-4, 0.2, 0.01
bth, c, k = 1e-4, 5.0, 200.0

sol1 = solve_ivp(
    fun=lambda t, z: rack_pinion_rhs(t, z, Jp, m, rp, bth, c, k, lambda tt: step_torque(tt, 0.05, 0.02)),
    t_span=(0.0, 0.5),
    y0=np.array([0.0, 0.0]),
    max_step=1e-3
)

theta1 = sol1.y[0]
x1 = rp * theta1

plt.figure()
plt.plot(sol1.t, x1)
plt.xlabel("t [s]")
plt.ylabel("rack displacement x [m]")
plt.grid(True)

# Crank-slider parameters
Jc, m2 = 3e-4, 0.3
r, l = 0.02, 0.08
bth2, c2, k2 = 2e-4, 3.0, 50.0

sol2 = solve_ivp(
    fun=lambda t, z: crank_slider_rhs(t, z, Jc, m2, r, l, bth2, c2, k2, lambda tt: step_torque(tt, 0.02, 0.02)),
    t_span=(0.0, 0.5),
    y0=np.array([0.5, 0.0]),   # start away from singular positions
    max_step=1e-4
)

plt.figure()
plt.plot(sol2.t, sol2.y[0])
plt.xlabel("t [s]")
plt.ylabel("crank angle theta [rad]")
plt.grid(True)
plt.show()
      

5.2 C++ (from-scratch RK4; optional Boost.Odeint / Eigen)

Recommended ecosystem: a minimal RK4 is sufficient for this course stage; for larger projects, consider Boost.Odeint for ODE integration and Eigen for linear algebra (later chapters).

#include <iostream>
#include <vector>
#include <cmath>
#include <functional>

struct State {
  double theta;
  double omega;
};

static inline double x_theta(double theta, double r, double l) {
  return r * std::cos(theta) + std::sqrt(l * l - (r * r) * std::pow(std::sin(theta), 2.0));
}

static inline double dx_dtheta(double theta, double r, double l) {
  double s = std::sin(theta);
  double c = std::cos(theta);
  double denom = std::sqrt(l * l - (r * r) * s * s);
  return -r * s - (r * r) * s * c / denom;
}

// Numerical derivative for d2x/dtheta2 and then dJ_eq/dtheta
static inline double dJeq_dtheta(double theta, double m, double r, double l) {
  double h = 1e-6;
  double d1 = dx_dtheta(theta + h, r, l);
  double d0 = dx_dtheta(theta - h, r, l);
  double d2x = (d1 - d0) / (2.0 * h);
  return 2.0 * m * dx_dtheta(theta, r, l) * d2x;
}

static inline double step_torque(double t, double tau0, double t0) {
  return (t >= t0) ? tau0 : 0.0;
}

// Crank-slider RHS: returns dtheta/dt and domega/dt
State crank_slider_rhs(double t, const State& z,
                       double Jc, double m, double r, double l,
                       double bth, double c, double k,
                       const std::function<double(double)>& tau_in_fun) {
  double theta = z.theta;
  double omega = z.omega;

  double Jmap = dx_dtheta(theta, r, l);
  double x = x_theta(theta, r, l);
  double xdot = Jmap * omega;

  double Jeq = Jc + m * (Jmap * Jmap);
  double dJeq = dJeq_dtheta(theta, m, r, l);

  double tau_in = tau_in_fun(t);
  double Fload = k * x + c * xdot;
  double tau_load = Fload * Jmap;

  double domega = (tau_in - 0.5 * dJeq * omega * omega - bth * omega - tau_load) / Jeq;
  return State{omega, domega};
}

// One RK4 step
State rk4_step(double t, double h, const State& z,
               const std::function<State(double,const State&)>& f) {
  State k1 = f(t, z);

  State z2{ z.theta + 0.5*h*k1.theta, z.omega + 0.5*h*k1.omega };
  State k2 = f(t + 0.5*h, z2);

  State z3{ z.theta + 0.5*h*k2.theta, z.omega + 0.5*h*k2.omega };
  State k3 = f(t + 0.5*h, z3);

  State z4{ z.theta + h*k3.theta, z.omega + h*k3.omega };
  State k4 = f(t + h, z4);

  State out;
  out.theta = z.theta + (h/6.0)*(k1.theta + 2.0*k2.theta + 2.0*k3.theta + k4.theta);
  out.omega = z.omega + (h/6.0)*(k1.omega + 2.0*k2.omega + 2.0*k3.omega + k4.omega);
  return out;
}

int main() {
  // Parameters
  double Jc = 3e-4, m = 0.3;
  double r = 0.02, l = 0.08;
  double bth = 2e-4, c = 3.0, k = 50.0;

  auto tau_in_fun = [&](double t){ return step_torque(t, 0.02, 0.02); };

  auto rhs = [&](double t, const State& z){
    return crank_slider_rhs(t, z, Jc, m, r, l, bth, c, k, tau_in_fun);
  };

  // Simulate
  double t0 = 0.0, tf = 0.5, h = 1e-4;
  State z{0.5, 0.0}; // start away from dead centers
  int N = static_cast<int>((tf - t0)/h);

  for (int i = 0; i <= N; ++i) {
    double t = t0 + i*h;
    if ((i % 1000) == 0) {
      std::cout << t << " " << z.theta << " " << z.omega << "\n";
    }
    z = rk4_step(t, h, z, rhs);
  }
  return 0;
}
      

5.3 Java (from-scratch RK4; optional Apache Commons Math ODE)

Recommended library for production-grade ODE solving: org.apache.commons:commons-math3. Here we provide a transparent RK4 implementation aligned with the C++ version.

import java.util.function.DoubleUnaryOperator;

public class CrankSliderRK4 {

  static class State {
    double theta;
    double omega;
    State(double th, double om) { theta = th; omega = om; }
  }

  static double xTheta(double theta, double r, double l) {
    return r * Math.cos(theta) + Math.sqrt(l*l - (r*r) * Math.pow(Math.sin(theta), 2.0));
  }

  static double dxDtheta(double theta, double r, double l) {
    double s = Math.sin(theta);
    double c = Math.cos(theta);
    double denom = Math.sqrt(l*l - (r*r) * s * s);
    return -r * s - (r*r) * s * c / denom;
  }

  static double dJeqDtheta(double theta, double m, double r, double l) {
    double h = 1e-6;
    double d1 = dxDtheta(theta + h, r, l);
    double d0 = dxDtheta(theta - h, r, l);
    double d2x = (d1 - d0) / (2.0 * h);
    return 2.0 * m * dxDtheta(theta, r, l) * d2x;
  }

  static double stepTorque(double t, double tau0, double t0) {
    return (t >= t0) ? tau0 : 0.0;
  }

  static State rhs(double t, State z,
                   double Jc, double m, double r, double l,
                   double bth, double c, double k,
                   DoubleUnaryOperator tauIn) {

    double th = z.theta;
    double om = z.omega;

    double Jmap = dxDtheta(th, r, l);
    double x = xTheta(th, r, l);
    double xdot = Jmap * om;

    double Jeq = Jc + m * (Jmap * Jmap);
    double dJeq = dJeqDtheta(th, m, r, l);

    double tauInVal = tauIn.applyAsDouble(t);
    double Fload = k * x + c * xdot;
    double tauLoad = Fload * Jmap;

    double dom = (tauInVal - 0.5 * dJeq * om * om - bth * om - tauLoad) / Jeq;
    return new State(om, dom);
  }

  static State rk4Step(double t, double h, State z,
                       java.util.function.BiFunction<Double, State, State> f) {

    State k1 = f.apply(t, z);

    State z2 = new State(z.theta + 0.5*h*k1.theta, z.omega + 0.5*h*k1.omega);
    State k2 = f.apply(t + 0.5*h, z2);

    State z3 = new State(z.theta + 0.5*h*k2.theta, z.omega + 0.5*h*k2.omega);
    State k3 = f.apply(t + 0.5*h, z3);

    State z4 = new State(z.theta + h*k3.theta, z.omega + h*k3.omega);
    State k4 = f.apply(t + h, z4);

    double th = z.theta + (h/6.0)*(k1.theta + 2.0*k2.theta + 2.0*k3.theta + k4.theta);
    double om = z.omega + (h/6.0)*(k1.omega + 2.0*k2.omega + 2.0*k3.omega + k4.omega);
    return new State(th, om);
  }

  public static void main(String[] args) {

    double Jc = 3e-4, m = 0.3;
    double r = 0.02, l = 0.08;
    double bth = 2e-4, c = 3.0, k = 50.0;

    DoubleUnaryOperator tauIn = (t) -> stepTorque(t, 0.02, 0.02);

    java.util.function.BiFunction<Double, State, State> f =
      (t, z) -> rhs(t, z, Jc, m, r, l, bth, c, k, tauIn);

    double t0 = 0.0, tf = 0.5, h = 1e-4;
    int N = (int)((tf - t0) / h);
    State z = new State(0.5, 0.0);

    for (int i = 0; i <= N; i++) {
      double t = t0 + i*h;
      if (i % 1000 == 0) {
        System.out.println(t + " " + z.theta + " " + z.omega);
      }
      z = rk4Step(t, h, z, f);
    }
  }
}
      

5.4 MATLAB (ODE45) and Simulink block structure

MATLAB: use ode45 (nonstiff) to simulate the ODE in \( \theta \). Simulink: implement the same ODE using integrator chains (two integrators for \( \theta \)), plus blocks computing \( x(\theta) \), \( dx/d\theta \), \( J_{\text{eq}}(\theta) \), and the torque balance.

function lesson4_crank_slider_demo
  % Parameters
  Jc = 3e-4; m = 0.3;
  r = 0.02;  l = 0.08;
  bth = 2e-4; c = 3.0; k = 50.0;

  tau0 = 0.02; tStep = 0.02;

  rhs = @(t,z) crank_slider_rhs(t,z,Jc,m,r,l,bth,c,k,tau0,tStep);

  tspan = [0 0.5];
  z0 = [0.5; 0.0]; % [theta; omega]
  opts = odeset('MaxStep',1e-4);

  [t,z] = ode45(rhs,tspan,z0,opts);

  theta = z(:,1);
  figure; plot(t,theta); grid on;
  xlabel('t [s]'); ylabel('theta [rad]');
end

function dz = crank_slider_rhs(t,z,Jc,m,r,l,bth,c,k,tau0,tStep)
  theta = z(1); omega = z(2);

  x = r*cos(theta) + sqrt(l^2 - (r^2)*(sin(theta)^2));
  denom = sqrt(l^2 - (r^2)*(sin(theta)^2));
  Jmap = -r*sin(theta) - (r^2)*sin(theta)*cos(theta)/denom;

  xdot = Jmap*omega;

  Jeq = Jc + m*(Jmap^2);

  % numerical derivative dJeq/dtheta (simple central diff)
  h = 1e-6;
  Jp = -r*sin(theta+h) - (r^2)*sin(theta+h)*cos(theta+h)/sqrt(l^2 - (r^2)*(sin(theta+h)^2));
  Jm = -r*sin(theta-h) - (r^2)*sin(theta-h)*cos(theta-h)/sqrt(l^2 - (r^2)*(sin(theta-h)^2));
  dJmap = (Jp - Jm)/(2*h);
  dJeq = 2*m*Jmap*dJmap;

  tau_in = (t >= tStep)*tau0;

  Fload = k*x + c*xdot;
  tau_load = Fload*Jmap;

  domega = (tau_in - 0.5*dJeq*(omega^2) - bth*omega - tau_load)/Jeq;
  dz = [omega; domega];
end
      

Simulink block outline (continuous-time):

• Two cascaded Integrator blocks: \( \ddot{\theta} \rightarrow \dot{\theta} \rightarrow \theta \).
• Compute \( x(\theta) \) and \( J(\theta)=dx/d\theta \) using Trigonometric Function, Product, Sqrt blocks.
• Compute \( \dot{x} = J(\theta)\dot{\theta} \).
• Compute \( J_{\text{eq}}(\theta)=J_c + m J(\theta)^2 \).
• Compute \( \tau_{\text{load}} = (k x + c \dot{x}) J(\theta) \).
• Compute \( \ddot{\theta} = \big(\tau_{\text{in}} - \tfrac{1}{2}J_{\text{eq}}'(\theta)\dot{\theta}^2 - b_\theta \dot{\theta} - \tau_{\text{load}}\big)/J_{\text{eq}}(\theta) \).

5.5 Wolfram Mathematica (symbolic + NDSolve)

Mathematica is especially useful here for symbolic differentiation of \( x(\theta) \) and automatic generation of \( dx/d\theta \) and \( dJ_{\text{eq}}/d\theta \).

ClearAll["Global`*"];

(* Parameters *)
Jc = 3*10^-4; m = 0.3;
r = 0.02; l = 0.08;
bth = 2*10^-4; c = 3.0; k = 50.0;

x[th_] := r Cos[th] + Sqrt[l^2 - r^2 Sin[th]^2];
Jmap[th_] := D[x[th], th];

Jeq[th_] := Jc + m (Jmap[th]^2);
dJeq[th_] := D[Jeq[th], th];

tauIn[t_] := Piecewise[{ {0.02, t >= 0.02} }, 0];

(* ODE: Jeq(th)*th'' + 1/2*dJeq(th)*(th')^2 + bth*th' + (k*x + c*xdot)*Jmap = tauIn *)
eqn =
  Jeq[th[t]] th''[t] + (1/2) dJeq[th[t]] (th'[t]^2) + bth th'[t] +
   (k x[th[t]] + c (Jmap[th[t]] th'[t])) Jmap[th[t]] == tauIn[t];

sol = NDSolve[
  {eqn, th[0] == 0.5, th'[0] == 0.0},
  th, {t, 0, 0.5},
  MaxStepSize -> 10^-4
];

Plot[Evaluate[th[t] /. sol], {t, 0, 0.5}, GridLines -> Automatic,
 AxesLabel -> {"t [s]", "theta [rad]"}]
      

Problem 1 (Rack-and-pinion reflection): A pinion with inertia \( J_p \) and pitch radius \( r_p \) drives a rack with mass \( m \) connected to a spring \( k \) and damper \( c \) to ground. A motor applies torque \( \tau_{\text{in}}(t) \) to the pinion shaft. Derive the ODE in \( \theta \).

Solution: Use the constraint \( x=r_p\theta \), so \( \dot{x}=r_p\dot{\theta} \), \( \ddot{x}=r_p\ddot{\theta} \). Spring/damper forces are \( F_s = kx \), \( F_d = c\dot{x} \). Convert to torque using \( \tau=r_pF \):

\[ \tau_s = r_p(kx) = (k r_p^2)\theta,\qquad \tau_d = r_p(c\dot{x}) = (c r_p^2)\dot{\theta}. \]

The rack mass reflects as inertia \( J_{\text{ref}} = m r_p^2 \), hence \( J_{\text{eq}}=J_p+m r_p^2 \). Including shaft viscous friction \( b_\theta\dot{\theta} \), torque balance yields:

\[ (J_p + m r_p^2)\ddot{\theta} + (b_\theta + c r_p^2)\dot{\theta} + (k r_p^2)\theta = \tau_{\text{in}}(t). \]

Problem 2 (Power mapping for a general constraint): Let an ideal mechanism impose \( x=f(\theta) \). Prove that the force–torque mapping is \( \tau = F \frac{df}{d\theta} \).

Solution: For an ideal (lossless) coupling, instantaneous power is conserved: \( P_{\text{rot}} = (\tau)\dot{\theta} \) and \( P_{\text{trans}} = F\dot{x} \). Since \( x=f(\theta) \), the chain rule gives \( \dot{x} = \frac{df}{d\theta}\dot{\theta} \). Therefore:

\[ (\tau)\dot{\theta} = F\dot{x} = F\left(\frac{df}{d\theta}\dot{\theta}\right) \quad \Longrightarrow \quad \tau = F\frac{df}{d\theta}. \]

Problem 3 (Crank–slider Jacobian): For \( x(\theta) = r\cos\theta + \sqrt{l^2 - r^2\sin^2\theta} \), compute \( \frac{dx}{d\theta} \).

Solution: Differentiate term-by-term. The first term gives \( \frac{d}{d\theta}(r\cos\theta)=-r\sin\theta \). For the square-root term, apply the chain rule:

\[ \frac{d}{d\theta}\sqrt{l^2 - r^2\sin^2\theta} = \frac{1}{2}\left(l^2 - r^2\sin^2\theta\right)^{-1/2}\left(-2r^2\sin\theta\cos\theta\right) = -\frac{r^2\sin\theta\cos\theta}{\sqrt{l^2 - r^2\sin^2\theta}}. \]

Summing yields:

\[ \frac{dx}{d\theta} = -r\sin\theta - \frac{r^2\sin\theta\cos\theta}{\sqrt{l^2 - r^2\sin^2\theta}}. \]

Problem 4 (Effective inertia and actuator-side equation): Assume a crank inertia \(J_c\) and slider mass \(m\) connected by an ideal crank–slider constraint \(x=f(\theta)\). Show that the kinetic energy can be written \( T=\frac{1}{2}J_{\text{eq}}(\theta)\dot{\theta}^2 \) with \( J_{\text{eq}}(\theta)=J_c+m\left(\frac{df}{d\theta}\right)^2 \), and derive \( J_{\text{eq}}(\theta)\ddot{\theta} + \frac{1}{2}\frac{dJ_{\text{eq}}}{d\theta}\dot{\theta}^2 = \tau_{\text{net}} \).

Solution: Using \( \dot{x}=\frac{df}{d\theta}\dot{\theta} \), the kinetic energy is:

\[ T=\frac{1}{2}J_c\dot{\theta}^2 + \frac{1}{2}m\dot{x}^2 = \frac{1}{2}J_c\dot{\theta}^2 + \frac{1}{2}m\left(\frac{df}{d\theta}\right)^2\dot{\theta}^2 = \frac{1}{2}\left(J_c+m\left(\frac{df}{d\theta}\right)^2\right)\dot{\theta}^2. \]

Define \( J_{\text{eq}}(\theta)=J_c+m\left(\frac{df}{d\theta}\right)^2 \) so \( T=\frac{1}{2}J_{\text{eq}}(\theta)\dot{\theta}^2 \). For ideal couplings, power balance gives \( \frac{dT}{dt} = \tau_{\text{net}}\dot{\theta} \). Differentiate \(T\) using the chain rule:

\[ \frac{dT}{dt} = \frac{1}{2}\frac{dJ_{\text{eq}}}{d\theta}\dot{\theta}\,\dot{\theta}^2 + J_{\text{eq}}(\theta)\dot{\theta}\ddot{\theta} = \left(J_{\text{eq}}(\theta)\ddot{\theta} + \frac{1}{2}\frac{dJ_{\text{eq}}}{d\theta}\dot{\theta}^2\right)\dot{\theta}. \]

Equate to \( \tau_{\text{net}}\dot{\theta} \) and cancel \( \dot{\theta} \) to obtain:

\[ J_{\text{eq}}(\theta)\ddot{\theta} + \frac{1}{2}\frac{dJ_{\text{eq}}}{d\theta}\dot{\theta}^2 = \tau_{\text{net}}. \]

Problem 5 (Rack displacement recovery and limiting check): For the rack-and-pinion ODE in Problem 1, show that if \(k=0\), \(c=0\), \(b_\theta=0\), and \( \tau_{\text{in}}(t)=\tau_0 \) constant, then \( x(t) \) is a quadratic function of time.

Solution: With \(k=c=b_\theta=0\), the ODE reduces to:

\[ (J_p+m r_p^2)\ddot{\theta} = \tau_0 \quad \Longrightarrow \quad \ddot{\theta} = \frac{\tau_0}{J_p+m r_p^2} = \text{constant}. \]

Integrating twice gives \( \theta(t)=\theta(0)+\dot{\theta}(0)t+\frac{1}{2}\ddot{\theta}t^2 \). Since \( x=r_p\theta \), we obtain:

\[ x(t)=r_p\theta(0)+r_p\dot{\theta}(0)t+\frac{1}{2}r_p\left(\frac{\tau_0}{J_p+m r_p^2}\right)t^2, \]

which is quadratic in \(t\), consistent with constant acceleration.

We modeled coupled translational–rotational mechanisms by (i) writing kinematic constraints \(x=f(\theta)\), (ii) enforcing power invariance to map force and torque, and (iii) reflecting mass, stiffness, and damping across the coupling. Rack-and-pinion yields constant reflected parameters and a linear ODE, while crank–slider yields configuration-dependent effective inertia and a nonlinear ODE in the actuator coordinate. These time-domain models are the foundation for later chapters on transfer functions, block diagrams, and state-space representations.