A collision is a short interaction in which bodies exert large forces on one another. The word "short" is relative to the surrounding motion: during the contact interval, external impulse is often small enough that the colliding bodies' total momentum can be treated as constant. Kinetic energy needs separate accounting.
For two bodies with constant masses, the exact system impulse-momentum balance is
\[m_1\mathbf v_{1f}+m_2\mathbf v_{2f}=m_1\mathbf v_{1i}+m_2\mathbf v_{2i}+\mathbf J_{\mathrm{ext}}\]Subscripts \(i\) and \(f\) refer to states immediately before and after contact. The external impulse is
\[\mathbf J_{\mathrm{ext}}=\int_{t_i}^{t_f}\mathbf F_{\mathrm{ext}}\,dt\]When this impulse is negligible, total momentum is conserved. That approximation can hold even when gravity acts, because a millisecond of weight impulse may be tiny compared with the contact impulse. It must be estimated, not assumed from the word "collision."
The contact impulse on body 1 and body 2 satisfies
\[\mathbf J_{12}=-\mathbf J_{21}\]for an isolated action-reaction pair. Each body's momentum changes, while their internal impulses cancel in the total.
Total translational kinetic energy before impact is
\[K_i=\frac12m_1v_{1i}^2+\frac12m_2v_{2i}^2\]An elastic collision conserves total kinetic energy of the modeled bodies as well as momentum:
\[K_f=K_i\]In an inelastic collision, translational kinetic energy decreases while total energy is still conserved when internal energy, permanent deformation, sound, heat, and other stores are included. A perfectly inelastic collision is the special case in which the bodies leave with the same velocity. "Perfectly" means maximally inelastic for the given one-dimensional momentum constraint; it does not mean all energy disappears.
A superelastic collision can have \(K_f>K_i\) if stored internal, chemical, or deformation energy becomes translational kinetic energy. Momentum can still be conserved.
If two bodies stick and external impulse is negligible, their shared final velocity is
\[\mathbf v_f=\frac{m_1\mathbf v_{1i}+m_2\mathbf v_{2i}}{m_1+m_2}\]This is the initial center-of-mass velocity. The kinetic-energy decrease in one dimension or for collinear relative motion is
\[K_i-K_f=\frac12\mu|\mathbf v_{1i}-\mathbf v_{2i}|^2,\qquad \mu=\frac{m_1m_2}{m_1+m_2}\]The reduced mass \(\mu\) is positive, so the decrease is nonnegative. It is zero only when the bodies initially have the same velocity, meaning there is no relative impact motion.
import numpy as np
m1, m2, u1, u2 = 2.0, 3.0, 4.0, 0.0
shared_velocity = (m1 * u1 + m2 * u2) / (m1 + m2)
initial_energy = 0.5 * m1 * u1**2 + 0.5 * m2 * u2**2
final_energy = 0.5 * (m1 + m2) * shared_velocity**2
print(f"Shared velocity: {shared_velocity:.3f} m/s")
print(f"Energy transferred to internal modes: {initial_energy-final_energy:.3f} J")
For a direct one-dimensional elastic collision with no external impulse, momentum and kinetic-energy conservation determine
\[v_{1f}=\frac{m_1-m_2}{m_1+m_2}v_{1i}+\frac{2m_2}{m_1+m_2}v_{2i}\]and
\[v_{2f}=\frac{2m_1}{m_1+m_2}v_{1i}+\frac{m_2-m_1}{m_1+m_2}v_{2i}\]These formulas apply to velocities along one common line. They are not component-wise rules for arbitrary two-dimensional impacts. Useful limits expose mistakes: equal masses exchange velocities; a light object bouncing from a much heavier stationary body nearly reverses its velocity; and if both initial velocities are equal, no closing collision occurs.
The relative-velocity form is simpler:
\[v_{2f}-v_{1f}=-(v_{2i}-v_{1i})\]The relative speed of separation equals the relative speed of approach.
For a frictionless impact along contact normal \(\hat{\mathbf n}\), a simplified normal restitution law is
\[e=-\frac{(\mathbf v_{2f}-\mathbf v_{1f})\cdot\hat{\mathbf n}}{(\mathbf v_{2i}-\mathbf v_{1i})\cdot\hat{\mathbf n}}\]when the denominator represents approaching motion and is nonzero. Ideal elastic impact has \(e=1\); perfectly inelastic normal response has \(e=0\). Many passive macroscopic impacts use values between zero and one, but the coefficient can depend on speed, temperature, shape, and material history. Active or energy-releasing impacts can produce an effective value above one.
Restitution constrains only the normal relative velocity. It says nothing by itself about tangential friction, spin, rolling, deformation duration, or energy divided among those modes.
Let the contact impulse on body 2 be
\(j\hat{\mathbf n}\)and the impulse on body 1 be opposite. For particles with no rotation,
\[j=-\frac{(1+e)(\mathbf v_{2i}-\mathbf v_{1i})\cdot\hat{\mathbf n}}{1/m_1+1/m_2}\]under a consistent normal convention. Apply an impact impulse only when bodies are approaching. Applying it to separating bodies produces a nonphysical second bounce.
def collide(m1, m2, u1, u2, restitution):
relative_speed = u2 - u1 # Normal points from body 1 to body 2.
if relative_speed >= 0:
return u1, u2
impulse = -(1 + restitution) * relative_speed / (1 / m1 + 1 / m2)
return u1 - impulse / m1, u2 + impulse / m2
v1, v2 = collide(m1, m2, u1, u2, 0.8)
print("Final velocities (m/s):", v1, v2)
Reversing the chosen normal reverses the signed scalar convention but must leave the physical velocity update unchanged when every related sign is transformed consistently.
Momentum supplies one conservation equation per spatial component. Kinetic energy adds one scalar equation for elastic impact, but these equations alone may still leave many possible final states. Contact geometry, impact parameter, smoothness, friction, and rotation provide additional information.
For ideal smooth spheres, impulse acts along the line joining centers at contact, so the tangential relative velocity is unchanged. Rough contact can apply tangential impulse and exchange translational and rotational energy. If spin is modeled, angular momentum and rotational kinetic energy must join the ledger.
Simultaneous contacts cannot always be resolved correctly as an arbitrary sequence of pairwise collisions. The chosen order can change the result. Dense granular systems and stacked rigid bodies need a coupled contact method or a clearly documented approximation.
In the center-of-momentum frame, total momentum is zero before and after an isolated collision. For a two-body elastic collision, each body's momentum reverses appropriately while its magnitude in that frame remains fixed. This often makes collision geometry easier to analyze.
Kinetic-energy change between the same initial and final states is independent of a Galilean shift between inertial frames when total mass and total momentum are conserved. The individual kinetic energies are frame-dependent, but the collision classification is not changed by selecting another inertial frame.
Evaluate momentum and energy from the exact same snapshots. A solver that checks one quantity before penetration correction and another afterward creates a meaningless comparison. Include external impulse explicitly and keep internal-energy transfer separate from numerical loss.
Finite time steps can miss fast collisions or allow deep overlap. Position correction can inject momentum or energy if it is mixed with the velocity impulse. Reduce the step, change the contact tolerance, and confirm that outgoing states converge. Event-driven methods must locate contact time accurately; time-stepping methods must document how overlap and repeated impulses are prevented.
momentum_residual = m1 * v1 + m2 * v2 - (m1 * u1 + m2 * u2)
energy_loss = initial_energy - (0.5 * m1 * v1**2 + 0.5 * m2 * v2**2)
restitution_residual = (v2 - v1) + 0.8 * (u2 - u1)
assert np.isclose(momentum_residual, 0.0)
assert np.isclose(restitution_residual, 0.0)
assert energy_loss >= 0
print("Momentum residual and energy loss:", momentum_residual, energy_loss)
Collisions are constrained transfers, not merely abrupt changes in velocity. Momentum tracks the system-wide impulse balance, while kinetic energy and restitution describe what the contact does with relative motion.
Two bodies of mass 2 and 3 kilograms approach at relative speed 4 metres per second. The question for this calculation is precise: how does kinetic energy lost change when we change coefficient of restitution, while keeping the other choices fixed? The earlier sections explain the physical relationship. Here we turn one part of that relationship into a calculation that you can run, inspect, and change in Google Colab.
The collision is isolated, one-dimensional, and frictionless; restitution lies between zero and one. Keep that sentence beside the results. A formula can return a number outside its useful range, and a plotting library will usually draw that number without asking whether the physical assumptions still hold. The selected interval is from 0 to 1 dimensionless. It gives us a specific domain to discuss instead of leaving the graph's range to a software default.
Use the notebook link above the lesson, run the setup cell, and then run the cells in order. The code below is also included there. NumPy handles the arrays of numerical values, SymPy keeps the equation symbolic, and Matplotlib draws the results. If these libraries or notebook cells are unfamiliar, start with Your first physics notebook.
In this small example, q stands for coefficient of restitution, expressed in dimensionless. The output is
kinetic energy lost, expressed in J. The equation below uses the numerical values in those
units. Any fixed coefficients therefore carry the units needed to make the relation
consistent. A value supplied in another unit must be converted before it enters this
expression.
This is a function of one changing input. It is not a statement that the physical system has only one variable. We deliberately hold the other inputs fixed so that one cause of variation is visible. Changing several inputs at once would make it harder to tell which change produced the output difference. Write down the held-fixed choices before trying a different experiment.
At 0.8 dimensionless, the expression gives 3.456 J. This is our reference calculation. It is a prediction of the stated example, not a new measurement. Its displayed digits help us compare implementations. They should not be read as a claim that physical inputs are known to six significant figures.
import numpy as np
import sympy as sp
import matplotlib.pyplot as plt
q = sp.Symbol("q", real=True)
expression = sp.sympify("9.6*(1-q**2)", locals={"q": q})
reference_input = 0.8
reference_output = float(expression.subs(q, reference_input))
print("Input:", reference_input, "dimensionless")
print("Predicted kinetic energy lost:", reference_output, "J")
sp.Symbol declares a mathematical symbol. It does not assign a numerical measurement
to that symbol. subs then replaces the symbol with our chosen input, and float
converts the resulting symbolic number into a standard Python numerical value.
Keeping these steps separate lets us inspect the equation before numerical substitution.
It also prevents a later plot from becoming the only surviving record of the model.
Before drawing hundreds of points, calculate three. The lower endpoint gives 9.6 J. The midpoint, at 0.5 dimensionless, gives 7.2 J. The upper endpoint gives 0 J. Read these together with the assumptions: they describe different choices of coefficient of restitution within the same example, with the other inputs held fixed.
These three values are a useful first check, but they do not establish the shape between them. Two endpoints can be equal even when the curve bends or oscillates in between. Likewise, a large change between endpoints does not tell us whether the relationship is linear. We need the equation and intermediate samples to answer those questions.
evaluate = sp.lambdify(q, expression, modules="numpy")
inputs = np.array([0.0, 0.5, 1.0])
outputs = np.asarray(evaluate(inputs), dtype=float)
outputs = np.broadcast_to(outputs, inputs.shape)
for input_value, output_value in zip(inputs, outputs):
print(f"{input_value:.6g} dimensionless -> {output_value:.6g} J")
assert np.all(np.isfinite(outputs))
lambdify turns the symbolic expression into a numerical function. We explicitly select
NumPy as its numerical backend so it can accept an array. The same expression now has
two useful forms: one for algebra and one for calculating many cases. The assertion
checks that the chosen inputs produce finite outputs. It would catch an infinity or an
undefined numerical result, but a finite value can still be physically inappropriate.
The horizontal axis shows coefficient of restitution; the vertical axis shows kinetic energy lost. Every plotted point is a prediction of the same expression. We mark the reference case separately so that a reader can connect the earlier printed value to its place on the curve. The axis labels include units because a shape without a scale cannot tell us the size of a physical effect.
input_grid = np.linspace(0.0, 1.0, 301)
output_grid = np.asarray(evaluate(input_grid), dtype=float)
output_grid = np.broadcast_to(output_grid, input_grid.shape)
fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(input_grid, output_grid, label="Stated model")
ax.scatter([reference_input], [reference_output], label="Reference case", zorder=3)
ax.set_xlabel("Coefficient of restitution (dimensionless)")
ax.set_ylabel("Kinetic energy lost (J)")
ax.legend()
ax.grid(alpha=0.25)
fig.tight_layout()
plt.show()
Look for a sign change, a turning point, or a region where the curve becomes steep. Connect each feature to the expression rather than to the colour of the line. If the curve has a minimum or maximum, changing the input in opposite directions may have similar effects on the output near that point. If it is nearly straight over a short interval, a local linear approximation may be useful even when the complete expression is nonlinear.
Increasing the number of samples makes the line smoother and may reveal features that a coarse grid missed. It does not improve the physical inputs or enlarge the model's validity range. For quantities that can take only discrete values, use the line as a guide to the relationship and interpret only the allowed input values. A plotting function does not enforce such restrictions for you.
A derivative asks a more specific question than the whole graph: how quickly does kinetic energy lost change near the reference input? Its units are the output unit divided by the input unit. For this calculation, the derivative at 0.8 dimensionless is -15.36 in those ratio units. Locally, the output decreases as the input increases. This statement is local; the sign or size of the slope can be different elsewhere on the curve.
\[f'(q) = - 19.2 q\]For a small input change, multiply the derivative by that change to estimate the output change. Then compare that estimate with a fresh evaluation of the full expression. The distinction between these two calculations matters: the derivative gives a tangent approximation, while substitution evaluates the chosen expression at the new point.
slope_expression = sp.diff(expression, q)
slope_at_reference = float(slope_expression.subs(q, reference_input))
changed_input = 0.805
linear_prediction = reference_output + slope_at_reference * (
changed_input - reference_input
)
full_prediction = float(evaluate(changed_input))
print("Symbolic derivative:", slope_expression)
print("Tangent prediction:", linear_prediction, "J")
print("Full expression:", full_prediction, "J")
print("Approximation difference:", full_prediction - linear_prediction, "J")
Here the changed input is 0.805 dimensionless. The tangent estimate is 3.3792 J, while the full expression gives 3.37896 J. If the two values are close, it means the chosen input change is small enough for this local approximation at this point. It does not mean the entire model is linear. If the slope is zero, a quadratic or higher-order term may supply the leading change instead.
This is also the beginning of uncertainty propagation. If the input has a small stated standard uncertainty, the magnitude of the slope times that uncertainty estimates one contribution to output uncertainty. Do not automatically use the same rule for a wide interval, a discontinuity, a turning point with zero first derivative, or several correlated uncertain inputs. Those situations need more of the model than one slope.
We can estimate the same slope from values on either side of the reference point. Subtract the lower output from the upper output and divide by twice the input step. Parentheses are essential: the entire output difference must be divided by the entire input interval. Without them, Python can execute a different mathematical expression while still producing a plausible number.
The calculation below uses three decreasing step sizes. Each estimate is compared with the symbolic derivative, so the reference is independent of the numerical difference formula. Near an endpoint, it is acceptable here to evaluate a small extension only because the displayed mathematical expression is smooth there; such a check does not assert that extrapolated values describe additional physical cases.
input_span = 1.0
for fraction in [1e-2, 1e-3, 1e-4]:
step = input_span * fraction
upper = float(evaluate(reference_input + step))
lower = float(evaluate(reference_input - step))
numerical_slope = (upper - lower) / (2 * step)
difference = numerical_slope - slope_at_reference
print(f"step={step:.6g}: slope={numerical_slope:.6g}, difference={difference:.6g}")
For a smooth function, reducing the step often reduces the approximation error at first. Very small steps can eventually make subtraction of nearly equal floating-point values less reliable. Therefore “smaller” is not an unlimited guarantee of “better.” Compare the observed behaviour over a useful range of steps instead of choosing the smallest number the computer can represent.
An absolute slope difference is useful when the expected slope is near zero. A relative difference would divide by a very small reference and could look enormous even when the absolute discrepancy is tiny. Away from zero, a relative difference can help compare quantities with very different scales. The tolerance should reflect the question and the numerical method, rather than a universal decimal threshold copied into every test.
Make one copy of the reference calculation before editing it. Choose a second value of coefficient of restitution inside the stated interval, predict whether kinetic energy lost will rise or fall, and then evaluate the expression. Place the prediction and the result in a text cell. If they disagree, identify whether the mistake was in the sign, the unit conversion, the algebra, or the assumption about how the curve behaves.
Next choose a shorter plotting interval around the reference point. Compare the local curve with its tangent approximation. Then widen the interval again. The difference shows why a useful local estimate need not describe the whole range. Keep the axis labels and numerical scale visible when comparing the plots; automatic axis limits can make two different effects appear equally large.
Finally, return to the physical limit stated at the start: The collision is isolated, one-dimensional, and frictionless; restitution lies between zero and one. Describe one change to the setup that would violate that limit. State what new input, equation, or observation would be needed before trusting an answer in the changed situation. This last step connects the Python experiment to the physics instead of treating successful execution as the only criterion for a useful result.
Before saving, restart the runtime and run all cells in order. A result that depends on an earlier, now-deleted input cell is not reproducible from the saved notebook. Keep unrounded values for later calculations, round only the reported output, and label any new data as measured, assumed, or generated. The graph and printed values should remain traceable to the same expression and the same choices of inputs.
There is another way to check that the expression and its derivative describe the same relationship. Add up the small output changes predicted by the slope across the chosen input interval. In the continuous limit, this integral equals the output at the upper endpoint minus the output at the lower endpoint. The result has units of J: the input unit in the integration cancels the input unit in the denominator of the slope.
This statement is about the mathematics of this particular relationship. It does not assign a new physical meaning to an arbitrary area under a graph. For example, an area under force versus distance represents work, while an area under a different pair of quantities may represent something else or have no standard name. Always carry the axis units through the calculation before naming an integral.
slope_function = sp.lambdify(q, slope_expression, modules="numpy")
endpoint_change = float(evaluate(1.0) - evaluate(0.0))
for points in [101, 201, 401]:
grid = np.linspace(0.0, 1.0, points)
slopes = np.asarray(slope_function(grid), dtype=float)
slopes = np.broadcast_to(slopes, grid.shape)
integrated_change = np.trapezoid(slopes, grid)
print(
points,
"integrated change:",
integrated_change,
"endpoint change:",
endpoint_change,
)
The trapezoidal rule joins neighbouring slope samples with straight segments and adds their signed areas. A negative contribution remains negative; replacing every slope by its magnitude would calculate a different quantity. Increasing the sample count provides a practical refinement check. If the values do not settle toward the endpoint difference, first inspect the interval, singularities, and array ordering.
Do not conclude that two agreeing calculations prove the model is physically complete. Both use the same selected relationship, so an omitted effect would be omitted from both. Their agreement is evidence that the numerical tools implement that relationship consistently. Comparing with an independent experiment or a more appropriate physical model is a separate test.
Keep this comparison in the notebook when you change the expression. It is particularly useful after adding a coefficient, changing a sign convention, or converting an input unit. A code cell that produces a smooth graph can conceal these mistakes; a second calculation of the same output change gives you another place to notice them.