How to Check Your Classical Mechanics Homework the Hard Way

In my third year of undergrad, in 2023, I was taking a classical mechanics class alongside an applied numerical methods class. For me, as for a lot of physics majors, classical mechanics was a first introduction to Lagrangian and Hamiltonian methods for finding a system’s equations of motion. If you know those equations and the initial conditions, you can work out what happens in the system at any given time. At least in principle.

The “old-school” way is to use Newtonian methods: define a coordinate system, draw free-body diagrams, write out the forces, and break them into components along your axes. Then use Newton’s second law, $\mathbf{F}=m\mathbf{a}$, to find the accelerations and equations of motion.

That’s fine for simple problems, but it can get annoying pretty quickly. You have to keep track of tension forces, normal forces, and other constraint forces, even when you don’t actually care what those forces are.

The Lagrangian approach trades difficult setup for a bit more algebra. Choose independent generalized coordinates that describe the motion, then write the kinetic and potential energies in terms of those coordinates, their velocities, any constants, and possibly time. For the ideal constraints in the examples here, holonomic (to use the fancy word from Goldstein’s pg. 12, 3rd ed.) and workless, choosing the coordinates takes care of the constraints. You don’t have to solve for every tension and normal force along the way.

With $L=T-V$ and the Euler–Lagrange equations, you can work out the accelerations. Hamilton’s equations describe the same motion using generalized coordinates and their canonical momenta. There’s still care to be had in choosing the coordinates and writing the energies, but a lot of what follows is algebra.

And, in the same way I tell the students I TA, “arithmetic should only be done in the comfort of one’s own home…” Algebra belongs there, too.

While I was learning all of this, I was also trying to come up with a final project for my applied numerical methods class. Finding equations of motion this way is pretty algorithmic, which seemed perfect for a symbolic homework checker. We were already using MATLAB in the numerical methods class, so that’s what I used.

I wrote a MATLAB program that takes symbolic expressions for a system’s kinetic and potential energies and works through Lagrange’s and Hamilton’s equations. It handles systems with one or two generalized coordinates: enough for a simple Atwood machine, a double pendulum, a double Atwood machine, or the block and wedge above.

I was reinventing the wheel, but if that wheel is rolling down a frictionful wedge sliding on a frictionless surface then I guess it’s worth it…

After presenting this project in my applied numerical methods class, I ended up with an A−.

Apparently there isn’t much numerical about doing everything symbolically. Who would have thought?

The block and wedge

The wedge can slide along the floor while the block slides down its hypotenuse. Both surfaces are frictionless. Once we choose generalized coordinates for the motion, we can write the total kinetic and potential energies in those coordinates and get started with Lagrange’s and Hamilton’s equations.

Letting MATLAB do the algebra

Once the energies are entered, MATLAB handles the derivatives and solves the resulting equations together.

% calculateEulerLagrange.m — first coordinate
% declare the coordinates and their velocities as symbolic functions of time.
syms q1(t) q1dot(t) q2(t) q2dot(t) t;

% take the partial derivatives of L with respect to q1 and its velocity.
partial_lagrangian_partial_q1 = diff(lagrangian, q1(t));
partial_lagrangian_partial_q1dot = diff(lagrangian, q1dot(t));
% the velocity derivative is the canonical momentum p1.
p1 = partial_lagrangian_partial_q1dot;

% differentiate that momentum with respect to time.
time_derivative_of_partial_lagrangian_partial_q1dot = ...
    diff(partial_lagrangian_partial_q1dot, t);

% replace derivatives of the velocity functions with acceleration symbols.
compact_time_derivative_1 = subs( ...
    time_derivative_of_partial_lagrangian_partial_q1dot, ...
    [diff(q1dot(t), t), diff(q2dot(t), t)], [q1ddot, q2ddot]);

% replace derivatives of the coordinates with velocity symbols.
compact_time_derivative_1_more = subs( ...
    compact_time_derivative_1, ...
    [diff(q1(t), t), diff(q2(t), t)], [q1dot, q2dot]);

% form the first euler–lagrange equation (the other coordinate is similar).
euler_lagrange_eq_1 = simplify( ...
    partial_lagrangian_partial_q1 - compact_time_derivative_1_more == 0);

% solveEulerLagrangeEquations.m — solve the coupled pair
% isolate the second acceleration, then substitute into the first equation.
sol2_for_q2ddot = solve(euler_lagrange_equation_2, q2ddot);
q1ddot_solution = simplify(solve(subs( ...
    euler_lagrange_equation_1, q2ddot, sol2_for_q2ddot), q1ddot));
% substitute back to obtain the second acceleration.
q2ddot_solution = simplify(subs(sol2_for_q2ddot, q1ddot, q1ddot_solution));

Run it in the browser!

For this article, I ported the MATLAB codebase into JavaScript so it can run in your browser. It follows the same workflow: enter the energies and let the solver do the algebra. Start with an example, inspect its energies, then change them or enter your own. A simple Atwood machine needs only one coordinate; the other examples use two. The playbacks at the top use the original preset energies.