Classical physics, read from a checker's source

Co-authored by Claude · Version 1.0.0, changed in rev 133

On 23 September 1999 the Mars Climate Orbiter was lost as it reached Mars. Its navigators had been steering it with figures from a ground program that reported the impulse of the spacecraft's thruster firings in pound-force seconds, contrary to the interface specification it was written to, which called for newton-seconds. Every figure was too small by a factor of 4.45, the number of newtons in a pound-force, and the spacecraft's course was misjudged. With the Mars Polar Lander, the separate spacecraft it was to relay for, the project had cost 327.6 million dollars, by NASA's figures. Sixteen years earlier, on 23 July 1983, an airliner ran out of fuel halfway through a flight over Canada and glided to a landing at Gimli, Manitoba; its fuel had been worked out with a density in pounds per litre where the calculation needed kilograms per litre, and it had taken off with about 45 per cent of the fuel it needed. Neither failure was a failure of physics. Each was a number handed over without its unit.

Physicists guard against it with an old habit. Every term of an equation must have the same dimensions as every other, mass with mass and length with length, a property called dimensional homogeneity, and the check is made before any number is put in. Joseph Fourier made the first credited contributions to the method in his Théorie analytique de la chaleur of 1822, the same book in which he developed the theory of the heat equation.

The physics library for DKE, the Deterministic Knowledge Engine, writes that habit down in DKE Python, and this essay teaches the physics by reading it. You record quantities as facts, each a value in a unit, and equations as products of quantities. The library's rules convert each quantity to SI, name a quantity handed to a reader who expects another unit, and find the dimensions an equation leaves over. Each answer names the facts it rests on, and correcting a fact withdraws the answers it supported. The second half of the library calculates: a pendulum, an orbit, the heat equation, a circuit and a magnet. Every program below runs as shown once the library's source file has been compiled into a store, and every output is copied from a run.

1. A unit as facts

A unit is two things: what it is made of, and how big it is. What it is made of is its dimension, written as powers of the base quantities. The library uses five of the seven SI base quantities, mass, length, time, electric current and temperature, and leaves out the amount of substance and luminous intensity, which no rule of it reads. A newton is a kilogram metre per second squared, so its dimension is mass to the power 1, length to the power 1 and time to the power −2. How big a unit is is its factor, the number of SI units of the same dimension it holds:

class Unit:
    mass: int
    length: int
    time: int
    current: int
    temperature: int
    factor: real

Most factors are exact by definition. The International Yard and Pound Agreement of 1959 defined the yard as exactly 0.9144 metres and the avoirdupois pound as exactly 0.45359237 kilograms. The inch had been exactly 25.4 millimetres in British standards since 1930, and the yard gives the same value, 0.9144 metres over 36. Standard gravity, 9.80665 metres per second squared, was set by the third General Conference on Weights and Measures in 1901, and a pound-force is the weight of a pound under it: 4.4482216152605 newtons. The thermochemical calorie is exactly 4.184 joules. Since 20 May 2019, when a revision approved by the 26th General Conference the previous November came into force, the SI itself has been built on seven defining constants, among them the elementary charge, exactly 1.602176634 × 10⁻¹⁹ coulomb; an electronvolt, the energy one elementary charge gains across one volt, is therefore exactly 1.602176634 × 10⁻¹⁹ joules.

si_units() writes these and some thirty more, each with its source. A quantity is then a value and the name of its unit, and the rules give it the unit's factor and powers. Its value in SI is a computed field, the value times the factor. This is the impulse of one thruster firing, recorded in the unit the ground program used:

import physics

physics.si_units()
physics.quantity("impulse", 1.0, "lbf_s")
print(physics.si("impulse"))

4.4482216152605

2. A number handed over

A quantity can also name the unit its reader expects. When it does, the rules compare the two units' dimensions and conclude what the hand-off is: convert when the dimensions agree and only the size differs, incompatible when they do not. The first of those is the rule:

@rule
def handoff_convert():
    for q in Quantity:
        if q.unit != q.expects and q.mass == q.emass and q.length == q.elength and q.time == q.etime and q.current == q.ecurrent and q.temperature == q.etemperature:
            q.handoff = "convert"

The navigation software expected newton-seconds. Recorded as such, the orbiter's hand-off is named, and the value the reader should have used is calculated:

import physics

physics.expect("impulse", "N_s")
for line in physics.checks():
    print(line)
print(physics.converted("impulse"))

handed over: impulse in lbf_s where N_s is expected, convert
4.4482216152605

A hand-off that converts is a mistake waiting to be made; one between different dimensions is a mistake already made, whatever the numbers. A force handed to a reader who expects an impulse is one, since an impulse is a force times a time:

import physics

physics.quantity("thrust", 2.0, "lbf")
physics.expect("thrust", "lbf_s")
for line in physics.checks():
    print(line)

handed over: impulse in lbf_s where N_s is expected, convert
handed over: thrust in lbf where lbf_s is expected, incompatible

3. Dimensions of an equation

An equation is checked by multiplying out its dimensions. The library writes one as a product of powers of quantities, with every power on the left side positive and every power on the right side negative, so that a consistent equation's powers of each base quantity add to zero. Newton's second law, F = m a, is F to the 1, m to the −1 and a to the −1. A pure number has no dimension and needs no factor, so the kinetic energy E = m v² / 2 is E to the 1, m to the −1 and v to the −2.

Each factor takes its quantity's powers from the quantity's unit, its share is those powers times its own power, and the equation sums the shares. The rules for mass are these, and length, time, current and temperature have the same:

@rule
def factor_mass():
    for f in Factor, q in Quantity, u in Unit:
        if f.quantity == q and q.unit == u:
            f.mass = u.mass
@rule
def equation_mass():
    for e in Equation, f in Factor:
        if f.equation == e and absent(e.mass_disputed):
            e.mass = sum(f.smass)

The last premise of the second rule holds the sum back while a factor is bound to the equation twice, which section 6 comes to.

Here are the two equations, and a third, momentum p = m v, with the momentum recorded by mistake in newtons:

import physics

physics.quantity("F", 12.0, "N")
physics.quantity("m", 3.0, "kg")
physics.quantity("a", 4.0, "m_per_s2")
physics.factor("newton2", "F", 1)
physics.factor("newton2", "m", -1)
physics.factor("newton2", "a", -1)
physics.quantity("E", 24.0, "J")
physics.quantity("v", 4.0, "m_per_s")
physics.factor("kinetic", "E", 1)
physics.factor("kinetic", "m", -1)
physics.factor("kinetic", "v", -2)
physics.quantity("p", 12.0, "N")
physics.factor("momentum", "p", 1)
physics.factor("momentum", "m", -1)
physics.factor("momentum", "v", -1)
print(physics.dimension("newton2"), physics.dimension("kinetic"), physics.dimension("momentum"))

consistent consistent T-1

The momentum equation leaves time to the power −1 over: a newton is a newton-second divided by a second, so the left side has a time too few. A dimensional check cannot say which side is wrong, only that they disagree, and it says nothing about the numbers, so the factor of one half in the kinetic energy is beyond it. That it cannot see is the price of its being so cheap. Lord Rayleigh proposed a "method of dimensions" in 1877, and the theorem behind such methods, that a law between n quantities made of k base dimensions can be written as a law between n − k dimensionless combinations of them, was proved by Joseph Bertrand in 1878 and generalised by Vaschy in 1892. It carries the name of Edgar Buckingham, whose paper "On physically similar systems" stated it again in 1914.

4. Correcting, and asking why

A correction is the same source recording a fact again. Record the momentum in newton-seconds, and the equation balances:

import physics

physics.quantity("p", 12.0, "N_s")
print(physics.dimension("momentum"))

consistent

Nothing else was rewritten. Correcting the unit withdrew the answers that rested on the newton, and the answers that rest on the newton-second took their place. why follows any conclusion back to the facts it rests on. The whole answer for the momentum equation runs to some sixty lines, one for every share and power; these are the conclusion and the chain under one of its factors:

import physics

w = why(physics.Equation.momentum.consistent)
for c in w:
    print(str(c.path) + " = " + str(c.value))

physics.Equation.momentum.consistent = True
physics.Equation.momentum.mass = 0
physics.Factor.momentum-p.smass = 1
physics.Factor.momentum-p.power = 1
physics.Factor.momentum-p.mass = 1
physics.Quantity.p.unit = N_s
physics.Unit.N_s.mass = 1

5. What a checker does not know

A unit the library has not been told about has no dimension, and an equation that uses it is neither consistent nor inconsistent. A checker that called it either would be guessing. The rule that marks such a factor says so:

@rule
def factor_unknown():
    # A factor whose quantity or unit is not recorded has no dimension to
    # share, so the equation it is in is not judged consistent.
    for f in Factor:
        if f.quantity != "" and absent(f.mass):
            f.unknown = True

Every rule that would judge an equation tests that none of its factors is unknown. A racehorse's pace, measured over a track in furlongs, is the first case here:

import physics

physics.quantity("track", 3.0, "furlong")
physics.quantity("lap", 40.0, "s")
physics.factor("pace", "v", 1)
physics.factor("pace", "track", -1)
physics.factor("pace", "lap", 1)
print(physics.dimension("pace"))
for line in physics.checks():
    print(line)

unknown unit
handed over: impulse in lbf_s where N_s is expected, convert
handed over: thrust in lbf where lbf_s is expected, incompatible
unknown unit: track in furlong

The furlong is an eighth of a mile, and writing it as a unit, with its length and its factor, would settle the question. The library leaves that to the reader rather than carry every unit ever used.

6. When two sources disagree

Two sources can give a fact two values, and a careful record keeps both rather than choosing. Suppose a second table rounds the pound-force to 4.45 newtons. The library goes on reasoning with both values, and unsettled() names what is in dispute:

import physics

update(physics.Unit("lbf").factor, 4.45, "rounded table")
print(physics.unsettled())
print(physics.si("thrust"))

['unit lbf']
8.896443230521

The thrust now has two values in SI, one from each factor, and si() reports the lesser. A query that refused to answer while anything was in dispute would leave its reader with nothing; one that quietly picked a winner would hide the dispute. Naming it is the third way. A sum is where reasoning with both values would mislead, because it would add both: a resistor whose resistance a second table writes, even the same resistance, would be counted twice in its group. So while a part is bound to its sum twice, the sum is disputed and has no value. resistance() refuses the group, disputed_groups() and disputed_equations() name what is held back, and forgetting the claim that is not to stand settles it.

7. A circuit

Georg Ohm published Die galvanische Kette, mathematisch bearbeitet in 1827, with the law that a conductor's current is proportional to the voltage across it, and in 1845 Gustav Kirchhoff stated the two laws that make circuits calculable: the currents into a junction sum to zero, and so do the voltages round a loop. Resistors in series carry one current, so their resistances add; resistors in parallel share one voltage, so their conductances, the reciprocals of their resistances, add. Those are two rules:

@rule
def series():
    # Resistors in series carry one current, so their resistances add.
    for g in Group, r in Resistor:
        if g.joined == "series" and r.group == g and absent(g.ohms_disputed):
            g.ohms = sum(r.ohms)

Its last premise holds the sum back while a resistor is bound to the group twice, as section 6 describes.

A circuit that is not built of series and parallel groups needs Kirchhoff's laws written out and solved. The bridge Samuel Hunter Christie invented in 1833, and Charles Wheatstone improved and popularised in 1843, is the classic case: two dividers side by side, with a meter across their middles that reads nothing when the ratios of the two dividers agree. bridge() writes the current law at each middle point and solves the pair of equations with solve(), by Gaussian elimination. The method appears, without proof, in the eighth chapter of the Chinese Nine Chapters on the Mathematical Art, Rectangular Arrays, and it was named for Gauss only in the 1950s, through confusion over its history. A singular system, one whose equations say the same thing twice, has no single answer, and solve() returns none:

import physics

physics.group("chain", "series")
physics.resistor("r1", 100.0, "chain")
physics.resistor("r2", 220.0, "chain")
physics.resistor("r3", 330.0, "chain")
physics.group("bank", "parallel")
physics.resistor("r4", 100.0, "bank")
physics.resistor("r5", 100.0, "bank")
physics.resistor("r6", 50.0, "bank")
print(physics.resistance("chain"), physics.resistance("bank"))
print(physics.bridge(10.0, 100.0, 200.0, 150.0, 300.0, 50.0))
print(physics.bridge(10.0, 100.0, 200.0, 150.0, 100.0, 50.0))
print(physics.solve([[2.0, 1.0], [4.0, 2.0]], [1.0, 2.0]))

650.0 25.0
-1.7763568394002505e-17
0.015094339622641506
[]

The first bridge balances, 100 is to 200 as 150 is to 300, and its meter current is zero to the last bit of floating-point rounding. The second does not, and 15 milliamperes flow through the meter.

8. The pendulum

In 1581 Galileo timed pendulums and found the property later called isochronism: a pendulum's swing takes the same time whether it is wide or narrow. On that property Christiaan Huygens built the first pendulum clock, on 25 December 1656, and in his Horologium Oscillatorium of 1673 he showed that it holds only approximately, and that wide swings make a pendulum keep bad time. For small swings the period is 2π√(L/g), whatever the amplitude. The exact period for a release from rest at an angle θ₀ divides that by the arithmetic–geometric mean of 1 and cos(θ₀/2), the limit of taking the ordinary and the geometric mean of a pair of numbers over and over; it converges in a handful of steps.

period() is that formula. pendulum() gets the same number without it, by stepping the equation of motion θ'' = −(g/L) sin θ forward in time and timing the swing down to the vertical, which is a quarter of the period. Here they are for a pendulum a metre long under standard gravity, at five amplitudes, beside the small-swing period:

import math
import physics

g = 9.80665
print(physics.small_period(1.0, g))
for deg in [5.0, 30.0, 90.0, 150.0, 179.0]:
    th = deg * math.pi / 180.0
    print(deg, physics.period(th, 1.0, g), physics.pendulum(th, 1.0, g, 0.0001))

2.0064092925890407
5.0 2.0073646894898 2.0073646894899255
30.0 2.0413384658583684 2.0413384658586597
90.0 2.3682463462860097 2.3682463462854573
150.0 3.5357019383114 3.5357019383116866
179.0 7.827133388800024 7.827133388798116

The two methods agree to about twelve figures. A pendulum swinging five degrees takes a twentieth of a per cent longer than the small-swing period; at thirty degrees it takes nearly two per cent longer, and released a degree from the top it takes nearly four times as long.

9. Orbits

Johannes Kepler published the first of his laws of planetary motion, that a planet moves on an ellipse with the Sun at a focus, in Astronomia nova in 1609, with an early form of the second; the second in its modern form, that the line from the Sun to the planet sweeps out equal areas in equal times, is in his Epitome Astronomiae Copernicanae of 1621, and the third, relating the period to the size of the orbit, in Harmonices Mundi of 1619. Isaac Newton derived all three from his law of gravitation in the Principia of 1687. A two-body orbit has an exact solution, so it is the test every numerical method for motion is first put to.

orbit() follows a body round an orbit of eccentricity 0.5 for ten periods, at a thousand steps a period, by three methods, and returns what an exact orbit keeps: the relative error in its energy, the relative error in its angular momentum, which is the second law, and the distance between where it ends and where it started. The first method is Leonhard Euler's, from his Institutionum calculi integralis of 1768 to 1770: step each quantity forward by its rate of change. The second is the one Loup Verlet used in 1967 for molecular dynamics, simulating a fluid of Lennard-Jones molecules. Jean Baptiste Delambre had used it first in 1791, Carl Størmer in 1907 for charged particles in a magnetic field, and Cowell and Crommelin in 1909 for Halley's comet. The third is the classical fourth-order method that Carl Runge in 1895 and Wilhelm Kutta in 1901 developed:

import physics

for method in ["euler", "verlet", "rk4"]:
    print(method, physics.orbit(0.5, method, 10, 1000))

euler [0.5818698318082012, 0.2559191304871933, 1.3576622634057045]
verlet [1.5051168311686062e-07, 1.230696119285481e-14, 0.017694255865162944]
rk4 [-2.96973601265904e-09, -4.5936860790439923e-10, 5.335604788350009e-07]

Euler's method has gained more than half the orbit's energy in ten periods, and ends 1.36 from where it began on an orbit whose semi-major axis is 1. The other two differ in a way that matters more than their errors' sizes. Runge–Kutta is the more accurate here. Verlet's misses the closing point by more, yet keeps the angular momentum to the last digits a double holds. The difference shows over a longer run. At a hundred steps a period, here is each method's energy error after 10, 40 and 160 periods:

import physics

for n in [10, 40, 160]:
    print(n, physics.orbit(0.5, "verlet", n, 100)[0], physics.orbit(0.5, "rk4", n, 100)[0])

10 0.010818061609004332 -0.00029554091594752297
40 0.010635894145922692 -0.001200345720125728
160 0.010625560042420762 -0.004739022723429629

Runge–Kutta's error grows in proportion to the time it runs. Verlet's is larger and stays where it is: the method preserves the geometry of the motion, its symplectic form on phase space, so its errors do not pile up in the quantities the motion conserves. For an orbit followed for millions of periods, that is the property that counts.

10. Heat

Fourier's book of 1822 is about how heat spreads, and the equation he developed for it says that a point warms at a rate proportional to how far it is below the average of its neighbours. heat() solves it on a rod held at zero at both ends and started as a triangle rising to one at its middle, by the plainest scheme: each step sets each point to itself plus r times its neighbours' sum less twice itself, where r is the time step over the square of the spacing. heat_exact() is Fourier's own answer, a series of decaying sines.

The scheme is stable if and only if r is at most one half, as a stability analysis due to John von Neumann shows. Below that it tracks the exact answer; above it, the smallest ripples grow by a factor at every step and swamp everything:

import physics

print(physics.heat(50, 0.4, 313), physics.heat_exact(313 * 0.4 / 2500.0))
print(physics.heat(50, 0.5, 250), physics.heat(50, 0.55, 227))

0.4955477955405207 0.4955140997880792
0.4963873608004604 -518765062471936.5

The first pair agrees to four figures. The last number is a temperature of minus five hundred trillion, from a rod that started between zero and one, after 227 steps at r = 0.55. John Crank and Phyllis Nicolson published a scheme in 1947 that is stable for every step size, at the price of solving a system of equations at each step.

11. A magnet

Wilhelm Lenz gave his student Ernst Ising a model of a magnet as a problem: a lattice of spins, each up or down, each preferring to agree with its neighbours. Ising solved it in one dimension in his thesis of 1924, found that the line never magnetises itself at any temperature above zero, and concluded wrongly that the model never does in any dimension. He never published again. In 1944 Lars Onsager solved the square lattice exactly. It does magnetise itself, below a critical temperature of 2 / ln(1 + √2), about 2.269 in units of the coupling over Boltzmann's constant. Onsager announced the formula for the magnetisation in 1948 and 1949 without a proof, and C. N. Yang published the first proof in 1952. At a temperature of 1.5 it gives a magnetisation of 0.98650 per spin.

ising() simulates the model by the method of the 1953 paper "Equation of State Calculations by Fast Computing Machines", by Nicholas Metropolis, Arianna Rosenbluth, Marshall Rosenbluth, Augusta Teller and Edward Teller, run on the MANIAC I computer, built under Metropolis at Los Alamos and completed in 1952. Credit for the method is disputed; Marshall Rosenbluth said that Metropolis played no part beyond providing computer time. It visits each spin in turn and flips it if that lowers the energy, and otherwise with a probability that falls with the energy cost and rises with the temperature. The random numbers come from Stephen Park and Keith Miller's "minimal standard" generator of 1988, which multiplies by 16,807 modulo the prime 2³¹ − 1; their paper's title was "Random Number Generators: Good Ones Are Hard To Find", and a generator called RANDU, widely used in the early 1970s, had shown why. Here is a 12 by 12 lattice below, at and above the critical temperature:

import physics

for t in [1.5, 2.269, 3.5]:
    print(t, physics.ising(12, t, 400, 1))

1.5 0.9856250000000014
2.269 0.7952083333333334
3.5 0.1793055555555556

At 1.5 the simulation is within a tenth of a per cent of Yang's formula. At the critical temperature an infinite lattice has no magnetisation at all, and above it none either; a lattice of 144 spins is too small to know that, and its magnetisation fades rather than vanishing. Onsager's exact solution is what tells a reader how far to trust the small one.

12. What the library leaves out

The library's physics is classical, and of that only a part. Its quantities are single numbers, not vectors, so it checks the dimensions of a force and not its direction. A unit whose zero is not SI's, such as the degree Celsius, is not a factor and is left out of its units altogether. It does not carry uncertainties, so it can say that two values disagree but not whether they disagree by more than their errors. It knows five base dimensions, and an equation in moles or candelas needs the other two written in. And its calculations are the textbook cases, each chosen because its answer is known some other way.

Closing

The Mars Climate Orbiter's ground program got its arithmetic right in its own units. What was missing was a check that travelled with the numbers. That is what recording a unit as a fact gives: every value carries what it is made of, every hand-off is compared, and every equation is summed before anyone trusts it. And because each conclusion names what it rests on, a reader who doubts one can ask why and follow the answer down to the units it was read from, the same way a physicist checks a colleague's working.