Letters Before Numbers

Write down the kinetic energy of a moving mass. There are two completely different things a computer can do with it.

The numeric reading says: tell me the mass and the speed, right now, and I will hand you back a number. Mass is 2.0. Speed is 3.0. Energy is 9.0. Done. The letters are gone the instant you supply values, and what survives is a single decimal.

The symbolic reading says: I will keep the letters. I know that m is a thing, that v is a thing, and that this expression is one half times m times v squared. I do not need to know what they are. I can still differentiate it, rearrange it, solve it for v, or hand you the exact answer with a square root sign still in it.

That is the whole distinction. Everything below is a consequence of it. The reason it matters in engineering is that a physics model is a statement about relationships, not a statement about one operating point. You want to know that torque scales with current, not merely that torque was 41.2 newton metres yesterday afternoon.

A numeric system stores a block of memory: 64 bits per value, laid out end to end. That is what a NumPy array is.

A symbolic system stores a tree. The expression is broken into operations and operands, and each operation becomes a node with children. One half times m times v squared is not a formula in a string, it is a Mul node whose children are the number one half, the symbol m, and a Pow node holding v and 2.

This sounds like a technicality and it is the single most useful mental model you can have. Once you see the tree, every symbolic operation stops being magic:

• Differentiation is a recursive walk. Hit a Mul node, apply the product rule and recurse into the children. Hit a Pow node, apply the power rule. There is no search and no cleverness required. This is why symbolic differentiation is fast and reliable.
• Substitution is a find and replace on the tree. Swap the m leaf for the number 2.0 and rebuild upward.
• Simplification is rewriting the tree into a smaller equivalent tree using known rules. This one is a search, which is why it is the slow, occasionally disappointing operation.
• Code generation is a walk that prints the tree in another language. Same tree, printed as C or as NumPy.

Floating point is not the real number line. It is about sixteen significant decimal digits, and the moment you subtract two nearly equal numbers, most of those digits are noise that was already there. This is called catastrophic cancellation, and it is the most common way a correct formula produces a wrong answer.

The classic demonstration is the expression below. Mathematically it approaches one half as x goes to zero. Numerically, if you type it the obvious way, it collapses to exactly zero.

Notice what happened. The fix was not a better solver or more precision. The fix was an algebraic identity, applied before any number was ever stored. That is the symbolic layer earning its keep. Goldberg's 1991 survey remains the standard reference on why this happens and what to do about it [1].

Symbolic math has its own failure mode, and it is the mirror image. Numbers stay the same size forever: a double is eight bytes whether it holds 3 or 3 billion. Expressions do not. They grow, and they grow fast.

The standard name is expression swell. Take a product of n functions and differentiate it n times. The product rule fires at every node, every time, and the term count explodes.

Swell is the reason the honest advice is not "use symbolic math for everything". It is: do the derivation symbolically, then leave. The symbolic layer is a workshop, not a factory floor. Section 28 shows the door you leave through.

A computer algebra system, or CAS, is any program that manipulates expression trees rather than numbers. Mathematica is one. SymPy is one. Maple, Maxima and MATLAB's Symbolic Math Toolbox are others. They differ enormously in polish and price, and hardly at all in what they are for.

Strip away the marketing and there are seven things you ask a CAS to do.

Notice the pattern. Six of the seven take an expression and return an expression. Only the last one leaves the symbolic world, and it leaves on purpose.

You do not need this section to use the tools, but knowing it explains every performance surprise you will hit.

Canonical form

The first job is deciding when two expressions are the same. x + x and 2*x must become one thing. Systems do this by defining an automatic canonical form: sort the arguments of commutative operations, collect like terms, flatten nested sums. This happens on construction, before you ask for anything, which is why typing x + x in SymPy immediately prints 2*x.

The deep problem lurking here is that for sufficiently rich classes of expressions, deciding whether something equals zero is undecidable. Practical systems use heuristics plus randomized numerical testing, which is why simplify occasionally fails to prove something you can see by eye.

Rewriting and pattern matching

Most of a CAS is rules. If you see sin(a)² + cos(a)², replace it with 1. The engine matches patterns against subtrees and applies rewrites, searching for a smaller result. This is a search over a huge space, which is why simplify is slow and why targeted functions such as trigsimp, collect or cancel are usually faster and more predictable.

Polynomial algebra and Grobner bases

Anything involving systems of polynomial equations rests on Grobner bases, introduced in Bruno Buchberger's 1965 doctoral thesis [2]. A Grobner basis is a rewritten form of a polynomial system that makes questions like "does a solution exist" and "eliminate this variable" mechanical. Solving nonlinear simultaneous equations, eliminating intermediate variables from a kinematic chain, checking observability of a polynomial model: all of it runs on this machinery. Worst-case cost is doubly exponential, which is your warning label.

Integration and the Risch algorithm

Symbolic integration is genuinely different from differentiation. Robert Risch showed in 1969 that for a large class of elementary functions there is a decision procedure: it will either produce an elementary antiderivative or prove that none exists [3]. Full implementations are enormous, so real systems layer heuristics first and fall back to Risch pieces. That is why some integrals return instantly, some grind, and some come back unevaluated.

A CAS is not a wish granting machine. Some things have no closed form and no amount of software will change that.

• General polynomials of degree five and above. The Abel-Ruffini theorem says there is no formula in radicals. Ask SymPy to solve a general quintic and it hands back CRootOf objects: exact, uniquely identified roots that can be evaluated to any precision, but not written out in surds. That is the correct answer, not a failure.
• Most integrals. The indefinite integral of exp(-x²) is not elementary. A CAS returns it in terms of the error function, which is really just a name for the thing that has no simpler name. The definite integral over the whole line, on the other hand, is exactly the square root of pi, and a CAS will tell you so.
• Most nonlinear differential equations. Closed form solutions are the exception. This is why the symbolic step usually ends at "here are the equations of motion" and hands off to a numerical integrator.

This is not a defeat. Knowing that no elementary antiderivative exists is itself a result, and it tells you immediately to reach for quadrature instead of hunting for a formula that was never there.

"Physics based" modeling, sometimes called white box modeling, means the model comes from a conservation law or a constitutive relation rather than from fitting a curve to data. Newton's second law, an energy balance, Kirchhoff's laws, Maxwell's equations, Fourier's law of conduction. You start from the statement of the physics and derive your way to something you can compute.

That derivation is the part people still do by hand on paper, and it is the part that is easiest to get wrong. A missing factor of two in a Park transform, a dropped cross term in a Jacobian, a sign flipped in a coordinate rotation: none of these announce themselves. They show up three weeks later as an estimator that will not converge.

The symbolic layer is where that derivation belongs.

The smallest honest example of a symbolic derivation. A mass on a rigid massless rod. We will not write Newton's laws with free body diagrams and force components. We will write down two energies and let the algebra produce the equation of motion.

This is the Lagrangian recipe, and it is the standard way multibody and driveline models get derived. Write the kinetic energy T, write the potential energy V, form L = T - V, then apply the Euler-Lagrange equation. What comes out the other side is the equation of motion, with every cross term and coupling included automatically.

The last line is the thing you wanted, and you never touched a force diagram. For a pendulum that is a parlour trick. For a five link suspension or a planetary gearset it is the difference between a model that exists and one that does not.

What the small angle approximation actually costs

The derived equation contains sin(theta), which is why the pendulum is nonlinear. Almost every control textbook then replaces sin(theta) with theta, using the series the CAS will hand you: sin(x) = x - x³/6 + x⁵/120 - .... That substitution is what makes the pendulum a linear oscillator with a constant natural frequency.

The series is also how you quantify the error, instead of hoping for the best.

A permanent magnet synchronous machine, written in rotor coordinates. This is the derivation that R. H. Park's 1929 two-reaction theory made possible [4], and it is a good example because the useful result is a formula, not a number.

In the rotor frame the flux linkages are the magnet flux plus the current times inductance in each axis. Torque is the cross product of flux and current, scaled. Written symbolically that is three lines, and expanding them gives the standard result that every field oriented controller is built on.

Read that back and it is the two term torque equation: a magnet term proportional to the magnet flux and the q-axis current, plus a reluctance term proportional to the difference between the two inductances and the product of both currents. A surface magnet machine has the two inductances nearly equal, so the second term vanishes and torque is simply proportional to q-axis current. An interior magnet machine has a real saliency, so the second term is worth having and you must go looking for it.

Now push the symbolic step one stage further. Write the currents in polar form with a total magnitude and a current angle, differentiate the torque with respect to that angle, and solve for zero. The result is the maximum torque per ampere angle, in closed form.

The third pattern is the one you will use most often: take a nonlinear physical model, and get the linear state space matrices that an observer or a controller needs. Symbolically, that is one function call.

Consider a two node thermal model. A heat source dumps power into node one, node one conducts to node two, node two conducts to ambient. Two capacitances, two conductances. Writing the energy balance for each node gives two first order differential equations, and the Jacobian of that pair with respect to the states is the A matrix.

Those matrices arrived with the physical parameters still in them. That is the whole prize. You can read the structure directly: the coupling term appears twice with opposite roles, the input only enters the first state, and every entry is a conductance over a capacitance, which has units of one over seconds. Substituting numbers in gives poles, time constants, and a discretized model for an estimator.

The single most common reason an engineer needs symbolic math is not integration or simplification. It is derivatives, and specifically Jacobians.

• An extended Kalman filter needs the Jacobian of the state transition and the measurement model at every step.
• Nonlinear model predictive control needs gradients of cost and constraints with respect to the decision variables.
• Parameter estimation needs sensitivities, the derivative of output with respect to each parameter, which is also what tells you whether a parameter is identifiable at all.
• Uncertainty propagation in the first order sense is just the Jacobian sandwiched around a covariance.

You can always approximate a derivative with finite differences. The problem is that finite differences force you to trade two errors against each other and there is no step size that removes both.

The failure mode of symbolic work is not a crash. It is a clean, plausible, well typeset expression that is wrong, because the physics you fed in was wrong. The algebra was flawless. The Lagrangian had a sign error.

These five checks are cheap, they are mostly automatable, and between them they catch the overwhelming majority of real errors.

1. The degenerate case. Set a parameter to the value where you already know the answer. Saliency to zero, and the reluctance torque term must vanish. Damping to zero, and the response must not decay. Coupling to zero, and the two subsystems must separate. This is the single highest yield check because you always have at least one such case.
2. Equilibrium and its sign. Find where the derivative is zero, then differentiate once more. The sign of that second derivative tells you whether the equilibrium is stable, and you know from physics which it should be. A hanging pendulum is stable. If your equation says otherwise, stop.
3. A conserved quantity. If the system is conservative, integrate it and watch the energy. With a correct model and a decent integrator the drift is at machine level. With a wrong model it is not.
4. Dimensions. Every term in a sum must carry the same units. This is entirely mechanical and a computer can do it for you, which is the subject of the next section.
5. A limit or asymptote. What happens as a parameter goes to zero or infinity? A CAS will take that limit exactly. The answer usually has an obvious physical reading, and if it does not, something is wrong.

Rather than describe this, here it is done. Four realistic errors were injected into the derivations from Part III, and each of the five checks was run against each. The results below are measured, not asserted.

Two things are worth pulling out of that.

First, the sign flipped pendulum is caught instantly by the equilibrium check: the stiffness comes out as plus g over L instead of minus g over L, which says the pendulum is stable pointing upward. You do not need a simulation to see that, and you do not need to be clever. You need to have asked.

Second, the sign flipped motor torque is not caught by peak torque, which is identical at 123.75 newton metres either way. It is caught only by looking at where the peak is: the correct model puts it at plus 24.356 degrees, the broken one at minus 24.356 degrees. One asks for negative d-axis current, which is what a real interior magnet machine wants. The other asks for positive d-axis current, which magnetizes the machine further and is physically backwards. A scalar check would have missed this entirely.

Dimensional analysis is the check engineers were taught first and automate last. There is no reason for that. SymPy carries a units module, and once quantities carry their dimensions a whole class of errors becomes impossible rather than merely unlikely.

Now mistype the capacitance as joules rather than joules per kelvin, the exact slip a tired person makes at five o'clock:

And an outright category error, adding a power to a heat capacity, raises rather than returning nonsense:

Dimensions do more than catch mistakes. They tell you how many independent parameters your model actually has, which is almost always fewer than the number you wrote down.

The Buckingham Pi theorem, from Edgar Buckingham's 1914 paper [22], says it precisely: if a physical relationship involves n variables built from r independent dimensions, it can be rewritten in terms of exactly n minus r dimensionless groups. Mechanically, you write the exponents of each variable as a column, and the null space of that matrix is your set of groups.

Applied to the two node thermal model from Section 11: seven quantities (heat input, two conductances, two capacitances, time, temperature rise), a dimensional matrix of rank three, and therefore four dimensionless groups. Not seven knobs. Four.

This one contradicts an instinct, so it is worth stating flatly. simplify optimizes for a small expression tree. It does not know anything about floating point. The shortest form and the numerically safest form are different objects, and sometimes they are opposites.

The quadratic formula is the canonical case. Both roots of a monic quadratic come out of the same expression with a plus or minus. When the linear coefficient is much larger than the constant, one of those two roots is computed by subtracting two nearly equal numbers, and Section 3 already told you what happens next.

The fix is Vieta's relation: the product of the roots equals the constant term. Compute the large root by the formula that does not cancel, then get the small root by division.

The general statement is that an expression has a condition number, a measure of how much a small relative change in the inputs can amplify into the output, and algebraically equivalent forms can have wildly different ones. A CAS will happily hand you either. Choosing between them is engineering judgement that the tool cannot make for you.

The rank tests for controllability and observability are usually taught as numerical procedures on numerical matrices. Done symbolically they become far more useful, because instead of a yes or a no you get the condition under which the answer changes.

Take the A and B matrices derived in Section 11 and build the controllability matrix. Its determinant comes out in one line:

Read that back and it says everything. The system is controllable for every physically meaningful parameter value, and it loses controllability exactly when the coupling conductance goes to zero, which is the physical statement that if the core is thermally disconnected from the case you cannot influence the case by heating the core. That is not a fact you would have learned from a numerical rank test at one operating point. It came out of the algebra.

Observability behaves the same way, and the same conductance is the culprit whichever node you measure:

The map from a state space model to a transfer function is a matrix inverse, a matrix product and a simplification, all of which a CAS does exactly. What arrives is a formula in your physical parameters rather than a table of coefficients.

Three readings fall straight out of that expression, and none of them would be visible in a numerical Bode plot.

• The steady state gain is G(0) = 1/k_2a. Every capacitance and the internal coupling cancel completely. At equilibrium the temperature rise is set by the conductance to ambient and by nothing else, which is obvious in hindsight and is exactly the kind of thing that is obvious only in hindsight.
• The numerator has no zeros, only the constant k_12. A two node thermal chain measured at the far end is a pure two pole lag. Any zero you see in measured data is telling you the model is missing a path.
• The middle coefficient mixes three products, which is the first hint of the identifiability problem in the next section: several different parameter sets can produce that same number.

Here is the question that decides whether a calibration campaign is worth running at all: given perfect, noise free data from the sensors you actually have, can the parameters be recovered uniquely?

Notice this is not a question about noise, sample rates or excitation. It is a question about the structure of the model, and it can be answered before a single measurement is taken. The standard treatment casts it as a differential algebra problem [23], and the machinery underneath is the same polynomial elimination from Section 6.

For the two node model measured at the case, the transfer function has exactly three independent coefficients. There are four unknown parameters. Three equations, four unknowns, and the conclusion is immediate: the parameters are not individually identifiable. Ask SymPy to solve anyway and it tells you so by leaving one free:

This sounds abstract until you look at what it costs. Below are three parameter sets, all physically plausible, all producing case temperature histories that agree to within a ten trillionth of a degree. They disagree about the core temperature by seventeen degrees.

A continuous model has to become a difference equation before it runs in a controller. The exact answer is the matrix exponential of A times the sample period, and a CAS will compute it symbolically. It is also a good demonstration of why you should not want it to.

For a single first order state the closed form is exp(-T/tau) and you should absolutely keep it symbolic. For two states it already runs to hundreds of terms full of atan2 and complex parts. The correct move is to keep A symbolic, substitute numbers, and exponentiate numerically. Deriving symbolically does not mean staying symbolic to the last possible instant; it means staying symbolic while the symbols are still telling you something.

The practical question is which approximation to use, and the honest answer depends entirely on the sample period relative to the fast time constant.

Everything so far has assumed you can write the model as derivatives on the left and functions on the right. Physical models built from connected components usually cannot. Connect two components and you get constraints, not derivatives: currents at a node sum to zero, two shafts share an angle, a bearing holds a distance fixed. What you have is a differential algebraic equation, and handing it to an ODE integrator does not work.

Write the pendulum in Cartesian coordinates and you have the smallest honest example. Two equations of motion with an unknown rod tension, plus one constraint that the length is fixed:

The fix is mechanical: differentiate the constraint with respect to time until a derivative of the missing variable appears. The number of differentiations needed is the differentiation index, and finding it automatically is what Pantelides' algorithm does [24], which is the machinery running silently inside every Modelica tool and inside Simscape.

The answer at the bottom of that ladder is worth a moment. The tension comes out as the mass times the speed squared minus gravity resolved along the rod, all over the length squared. That is centripetal force minus the weight component, exactly as it should be, and it was never written down by anyone. It was produced by differentiating a geometric constraint twice and substituting.

An expression that is correct can still be a bad program. Two transformations do most of the work of fixing that, and both are one function call.

Horner form

A polynomial written out as a sum of powers evaluates each power separately. Nested multiplication evaluates the same polynomial with only multiplies and adds, and no exponentiation at all. The operation count barely changes. The runtime does not.

Common subexpression elimination

A Jacobian repeats the same subtrees over and over, because the same inner functions appear in every entry. cse finds them and hoists them into temporaries. On a three by two Jacobian of a fairly ordinary nonlinear model, it removed two thirds of the arithmetic.

The full path from a physical law to a function running on an automotive microcontroller has five steps, and SymPy covers three of them.

1. Derive and verify using everything in Parts III and IV. Nothing goes further until the degenerate cases pass.
2. Reduce. cse to hoist repeats, horner on polynomials, and collect to group by the variable that changes fastest. Look at count_ops before and after.
3. Generate. codegen emits a compilable C function and its header. ccode emits a single expression. Both accept a settings object that controls how constants and integer powers are printed.
4. Check the numerics on the target's arithmetic. This step is yours, not the tool's. Single precision on an embedded target has about seven significant digits, not sixteen, so every conditioning concern from Section 17 gets worse by nine orders of magnitude. Re-run the cancellation checks in float32 before you trust anything.
5. Wire the generated file into the build as an artifact, never as source. The derivation script is the source. If a parameter definition changes, you regenerate and the diff is reviewable. Patching generated code by hand is how a model and its implementation quietly diverge.

Note what SymPy did with the one half: it printed (1.0/2.0) rather than 0.5. That is a deliberate choice about how the C compiler will read the constant, and it is the kind of detail that separates a code generator from a pretty printer.

The first thing to untangle is that "Wolfram" and "Mathematica" are not two products, and people use the words interchangeably in a way that hides the structure.

• Wolfram Language is the language and the engine. Everything else is a wrapper around it.
• Mathematica is the desktop notebook application built on that language. It is what most people mean when they say "Wolfram". It is commercial, with academic and home licences considerably cheaper than the professional one.
• WolframAlpha is a natural language front end to a curated slice of the same engine plus a large curated data set. It is a consumer product, not a modeling environment.
• Wolfram Engine is the language without the notebook, free for developers, callable from Python through the wolframclient package.
• Wolfram SystemModeler is a separate Modelica based physical modeling tool that talks to the language.

The language itself is unusual and worth understanding on its own terms: everything is an expression, evaluation is term rewriting, and pattern matching is the primary control structure. A definition like f[x_] := x^2 is not a function in the C sense, it is a rewrite rule that fires whenever the pattern matches. That design is exactly why the system is so strong at symbolic work, and exactly why its performance is hard to predict.

Where it genuinely leads: symbolic integration and differential equations, special functions, exact and arbitrary precision arithmetic with real numerical rigour, and simplification. When SymPy returns an expression unsimplified and Mathematica returns three terms, this is usually why. Version 15 shipped in June 2026 with a built in AI assistant, an expanded time series framework, and additions on the algebra side including Grassmann and Clifford algebras [6].

The costs are real: it is proprietary, licensing is a recurring negotiation, and the language does not embed comfortably in an existing Python or C toolchain. In a fleet analytics or embedded workflow, that last point usually decides it.

SymPy is a computer algebra system written in pure Python, first released in 2007 and described in a 2017 paper in PeerJ Computer Science [7]. It is free, BSD licensed, and has one design decision that dominates everything else: it is a library, not an application. There is no notebook, no separate kernel, no license server. You import sympy in the same file where you import numpy.

That single property is why it wins in an engineering pipeline even though it loses to Mathematica in raw symbolic strength. The symbolic derivation happens in the same process, same repository, same version control, and same test suite as the numerical code that consumes it.

# derive, simplify, and emit runnable code, in one file
import sympy as sp
import numpy as np

x, y = sp.symbols('x y', real=True)
f = sp.Matrix([sp.sin(x)*y**2, sp.exp(-x)+y])

J = f.jacobian([x, y])            # exact 2x2 Jacobian, still symbolic
J_fast = sp.lambdify((x, y), J, 'numpy')   # leave the symbolic world

J_fast(0.3, 1.2)                      # plain numpy array, microseconds

Practical notes that will save you an afternoon:

• Assumptions matter. sp.sqrt(n**2) stays unsimplified for a bare symbol, because n could be negative. Declare n = sp.Symbol('n', positive=True) and it simplifies to n. Most "why will it not simplify" complaints are missing assumptions.
• Keep floats out. Writing 0.5 inside a symbolic expression injects a float and blocks exact cancellation. Use sp.Rational(1,2).
• Use targeted simplifiers. simplify is a search. trigsimp, cancel, collect, factor and together are direct, fast, and predictable.
• Reach for cse before code generation. Common subexpression elimination pulls repeated subtrees into temporaries and can shrink generated code dramatically.
• It emits more than Python. ccode, fcode, octave_code, rust_code and the full codegen module will write C, Fortran or MATLAB. This is the path from a symbolic derivation to something an embedded target runs.

The honest limits: SymPy is slower than a compiled CAS, its simplify gives up earlier, and hard integrals and differential equations are noticeably weaker than Mathematica's. For the derivations engineers actually do, Jacobians, coordinate transforms, series expansions, solving for a variable, none of that is the binding constraint.

NumPy does not belong in the same comparison. It has no symbols, no expression trees, and no algebra. Putting it in a list with Mathematica and SymPy is like listing a lathe alongside two CAD packages. It is the thing that runs afterward.

What NumPy provides is one data structure, the n-dimensional array, and a large set of operations that act on the whole array at once [8]. Three properties do all the work:

• Contiguous typed memory. A million float64 values occupy exactly 8 MB in one block. The same numbers in a Python list cost roughly four times that, scattered across the heap as individual objects.
• Vectorization. np.sin(x) on an array runs the loop in compiled C, not in the interpreter. On a million points that is roughly a 20 to 30 times speedup over the equivalent Python loop, measured on the machine that produced this page.
• Broadcasting. Arrays of different but compatible shapes are stretched against each other without copying, so a parameter sweep over a grid is one expression rather than two nested loops.

SciPy sits directly on top and supplies the things a physics pipeline needs next: integrators, optimizers, linear algebra beyond the basics, signal processing, and interpolation [9]. In practice "NumPy" as a shorthand usually means the NumPy and SciPy pair.

Here is the mistake that motivates this section. Someone builds a beautiful symbolic model, then evaluates it by calling .subs() in a loop over their data. It works. It is also about fifty thousand times slower than it needs to be, and they conclude that symbolic math is impractical.

subs rebuilds an expression tree per call. That is the wrong tool for evaluation. lambdify compiles the tree once into a plain Python function whose body is ordinary arithmetic, closed over NumPy's array operations. After that there is no symbolic machinery in the loop at all.

That is all it is. No tree, no dispatch, no rewriting. The names sin and exp are bound to NumPy's versions, so passing an array evaluates the whole array in compiled code.

The relationship in one line: Mathematica and SymPy are alternatives to each other, and NumPy is what either of them hands off to. Wolfram Language is the engine Mathematica is built on, and it happens to include a strong numeric side of its own, which is why it can be a whole environment rather than one stage of a pipeline.

Once you understand the symbolic and numeric split, a whole layer of tools stops looking like a random list and starts looking like answers to specific gaps.

CasADi deserves a second mention because it sits precisely on the seam this page is about. Its expressions are symbolic enough to be differentiated exactly, and restricted enough to be compiled to fast C. It gives up general simplification and integration, which it does not need, in exchange for derivatives and speed, which it does.

Everything so far started from a known law. Sparse identification of nonlinear dynamics, or SINDy, starts from measurements and produces the law [14]. The premise is simple and usually true: most physical systems have sparse dynamics in a sensible basis, meaning that of the hundreds of candidate terms you might write down, only a handful actually appear.

The recipe is three steps. Build a library of candidate terms evaluated on your data. Estimate the derivatives. Then solve for the coefficients with a regression that is pushed toward zeros, typically by least squares followed by thresholding, repeated.

Below is a real run against the Van der Pol oscillator: eight trajectories, ten candidate terms up to cubic, derivatives from finite differences on the states rather than from the true vector field.

What you get out is an equation, not a black box. That is the entire point, and it is why SINDy sits in a diagnostics workflow in a way that a neural network does not: you can read the result, check it against physics, and reject it when the terms are nonsense.

This idea is nearly a century old, from Bernard Koopman's 1931 paper, and it is one of the strangest and most useful things in dynamics [26]. It says: stop tracking the state, and track functions of the state instead. In that space, nonlinear dynamics become linear. Exactly linear, not approximately.

The catch is that the space of functions is infinite dimensional, so you have traded a hard finite problem for an easy infinite one. The practical question is whether some finite set of observables closes, meaning the operator maps that set into itself.

Sometimes it does exactly. Take a system with a slow manifold, the standard example in the field. It is genuinely nonlinear. Now watch what happens if you add the square of the first state as a third observable.

That matrix was not guessed. It comes from differentiating each observable along the flow and expressing the result back in the same three observables, which is a symbolic operation and takes four lines of SymPy. The third row is the interesting one: the derivative of x1 squared is twice x1 x1 prime, which is twice mu times x1 squared, so the observable closes on itself.

Extended dynamic mode decomposition, or EDMD, is the data driven version [27]: pick a dictionary of observables, lift your snapshot data through it, and least squares fit the matrix that advances one lifted snapshot to the next. Run against the system above, it recovers the exact operator.

Two things follow that are worth carrying away.

First, the choice of observables is a modeling decision, not a hyperparameter, and physics is what tells you what to pick. The reason the square works here is that the nonlinearity in the equation is a square. When people say EDMD is finicky, what they usually mean is that they chose a dictionary by convenience rather than by looking at the equations.

Second, once the model is linear you get the whole linear toolbox back: pole placement, LQR, Kalman filtering, all the rank tests from Part V. That is the practical prize, and it is why Koopman methods keep appearing in fault detection work. A linear observer on a lifted space can do things a linearized observer on the original space cannot, because the lift is exact rather than a first order truncation valid near one operating point.

SINDy fixes the library and picks coefficients. Symbolic regression drops that restriction and searches over the space of expression trees themselves, which is enormous, discontinuous, and cannot be gradient descended. Modern implementations use evolutionary search, and PySR is the one most people reach for [25]. The same group showed that constraining a network's architecture first and reading expressions out of its learned components recovers known force laws and generalizes better than the network it came from [15].

The output is not one equation. It is a front: for each level of complexity, the most accurate expression found at that complexity. Reading that front is the skill, because the tempting move, taking the most accurate one, is almost always wrong.

Here is a real front, on data you have already met. Figure 10 showed the pendulum period error against amplitude, and the exact answer involves a complete elliptic integral. What is the best short formula for it?

The result is a small lesson in why this method exists. The four term Taylor series, at complexity 13, achieves a root mean square error of 8.2 thousandths. A simple rational form with one fitted coefficient, at complexity 6, achieves 7.4 ten-thousandths. Less than half the complexity, and eleven times more accurate.

Taylor expansion could never have produced that form, because it is not a polynomial. The search found it because the search was not restricted to polynomials. That is the entire argument for symbolic regression in one comparison, and it is also the argument against trusting a series expansion just because it came from an exact derivation.

The fourth approach runs the other direction. Rather than extracting an equation from data, a physics informed neural network takes an equation you already have and puts it into the training objective. The network is trained to match whatever data exists and to make the residual of the governing differential equation small at points where there is no data at all, with those residual derivatives supplied by automatic differentiation [13].

The symbolic equation never becomes a number here. It stays a statement about what any admissible solution must satisfy, and it acts as an infinite supply of free training data in the regions your sensors do not reach. That is its real value: extrapolation constrained by physics rather than by the shape of the data.

The through line is the same in all four. A purely learned model that cannot be inspected is hard to certify and harder to argue with. A purely first principles model that ignores fleet data is usually wrong at the edges, which is exactly where the failures live. Each of these methods is an attempt to keep a readable symbolic object while letting data supply what derivation cannot.

A diagnostic or prognostic system is a chain, and the two sides of this page do not compete across it. They occupy different links, and putting either one in the wrong link is the usual reason these systems disappoint.

The short version: symbolic work decides what to compute, numeric work computes it. Every question of the form "can this be detected", "which faults can be told apart", "which sensor is worth adding", "what does this parameter do to remaining life" is a question about structure, and structure lives in the letters. Every question of the form "is it happening now, on this vehicle, at this mileage" is a question about numbers.

A residual is a quantity that is zero when the system is healthy and non-zero when it is not. The naive way to build one is to run a model alongside the real thing and subtract. That works, and it requires knowing every state, including the ones you do not measure.

The better way is to eliminate the unmeasured states symbolically, leaving a relation among measured signals only. In the fault diagnosis literature this is an analytical redundancy relation [31], and building one is exactly the polynomial elimination from Section 6. Isermann's survey is the standard map of the model-based family this belongs to [35].

Take the two node thermal model. Solve the second equation for the unmeasured core temperature and substitute into the first:

Look carefully at what came out. Those coefficients are the denominator of the transfer function from Section 19, applied to the measured temperature, minus the numerator applied to the known input. That is not a coincidence and it is worth stating as a rule: the analytical redundancy relation is the transfer function written as a differential operator and moved to one side. If you already derived one, you already have the other.

Run the healthy model through it and the residual is zero to numerical noise. Run four faults through it, using the healthy parameters in the relation, and every one of them lights it up.

Detection is the easy half. The residual said something is wrong and could not say what. For that you need more than one relation, which means more than one sensor, and the question of which sensor is again symbolic.

Add a core temperature sensor and the model splits into two relations, one per node, each testable on its own:

That table is the fault signature matrix, and it was produced by asking SymPy whether each partial derivative is identically zero. Read the columns and it tells you exactly what the second sensor bought and what it did not.

• The coupling conductance is uniquely isolable within this fault set. It is the only one of the four that moves both relations, and it moves them in opposite directions, because its two sensitivities are exactly plus and minus the temperature difference. That qualifier matters: add a bias fault on the core temperature sensor and it produces a residual in exactly the same direction, at which point the two are no longer separable at all. The companion page works that out.
• Lost thermal mass and extra heat generation are still indistinguishable. Both show the pattern 10. No amount of thresholding will separate them.
• Fouled heat sink and a wrong case capacitance are likewise a pair. Pattern 01 for both, and this one the steady state test does not break either, because a conductance to ambient does appear in the heat balance while a capacitance does not. That pair needs a third relation.

The way out of the first collision is not another sensor. It is a different test on the sensors you have. At steady state both derivatives vanish, so the first relation collapses to a heat balance in which a capacitance cannot appear at all. Simulating that out to equilibrium confirms it exactly.

The procedure above is general. Here it is applied to four real components, with the symbolic step written out and the consequence measured.

The one worth dwelling on

The permanent magnet case produces a result in one line that changes how the monitor is designed. Differentiate the torque equation with respect to magnet flux and you get three halves times pole pairs times the q-axis current, and nothing else. At zero q-axis current that sensitivity is exactly zero, so a demagnetization monitor built on torque is blind at cruise, at coast, and at standstill. That is most of a drive cycle.

The back-EMF route has sensitivity equal to electrical speed instead, which is non-zero whenever the machine turns. The two routes cover complementary parts of the operating envelope, and knowing that requires exactly one symbolic differentiation of a formula you already had.

That last pair of numbers is a warning about acting on the wrong signal. The optimal current angle moves measurably as the magnet weakens, which makes it a decent health feature. But failing to retune it costs only two tenths of a percent of torque, because the torque curve is flat near its peak. The angle drift is worth measuring and not worth correcting, and only the symbolic form makes that distinction visible.

Everything so far assumed you observe the plant. On a modern vehicle you almost never do. You observe a plant wrapped in a controller that has been designed, with considerable care, to make faults invisible.

Consider an actuator whose gain degrades slowly, inside a loop with integral action. Work out the closed loop symbolically and the problem states itself:

The output gain is exactly one for any positive actuator gain. That is what integral action is for. It also means the tracking error, the residual everybody reaches for first, contains no information about the fault whatsoever. Meanwhile the control effort goes as one over the gain, and has been reporting the degradation faithfully since the first hour.

Prognostics splits into two approaches that get conflated constantly. This section is the first: you know the physics of the degradation, you write the rate law, and you integrate it to a threshold. This is the physics-of-failure route, and where the integral has a closed form, the answer is a formula rather than a simulation.

Crack growth is the standard example. Paris and Erdogan's 1963 law says the crack extension per cycle is a power of the stress intensity range [33]. Substituting the stress intensity and integrating is one line:

SymPy returns a piecewise result and flags the m = 2 case separately, where the integral becomes a logarithm. That is the tool being careful on your behalf about a special case most textbook derivations skip.

With a representative steel the closed form gives 780,125 cycles, which agrees with numerical quadrature of the same integral to sixteen digits. The number is not the interesting part. This is:

The same structure appears everywhere in reliability, which is why it is worth recognizing. Every one of the standard damage laws is a power law or an exponential, from Archard's wear relation [37] to the Coffin-Manson low cycle fatigue law [38], and the log derivative is therefore a constant that you can read directly off the symbol without simulating anything.

Read that table across and it says a ten percent load reduction on a gear buys two and a half times the life, while the same reduction on a bearing buys about 1.4 times. Nobody simulated anything. The exponent is the engineering answer, and it exists only because the law was kept symbolic.

Temperature, and the trap in averaging it

For thermally activated mechanisms the knob is temperature and the law is exponential, which makes averaging dangerous in a specific way. Take a duty cycle with 55 percent of the time at 45 degrees, 30 at 70, 12 at 95 and 3 at 120. The arithmetic mean is 60.75 degrees. The damage-equivalent constant temperature, the one that produces the same integrated Arrhenius damage, is 76.89 degrees.

Using the arithmetic mean understates the damage rate by a factor of three. The three percent of time spent at 120 degrees does more harm than the fifty five percent at 45, and the only reason anyone knows that is that the rate law was integrated rather than averaged.

The running thermal example inverts just as cleanly. Let the heat sink foul exponentially, take the steady state core temperature, and solve for the time at which it reaches a limit:

The second route applies when you do not know the degradation trajectory, only that some parameter is drifting. You promote that parameter to a state, let it random walk, and estimate it jointly with the physical states from whatever you can measure. Then you project the estimate forward to the threshold.

The question to answer first is whether the augmented system is observable at all, and that is symbolic. Take Lie derivatives of the measurement along the augmented dynamics and check the rank:

This is worth reconciling with Section 20, which said the parameters were not identifiable from the case temperature. Both statements are true and the difference is what is held known. Section 20 left all four parameters free and had three transfer function coefficients to determine them. Here the two capacitances are known from mass and specific heat, and only the conductance drifts, so there is information to spare. Fixing a parameter from an independent source, which Section 20 recommended, is precisely what makes the rest recoverable, and the rank test proves it rather than assuming it.

With observability settled, the filter is ordinary. The Jacobian it needs is the symbolic one, extended by a column for the parameter.

A remaining useful life with no interval on it is not an engineering result. The symbolic form gives you two ways to attach one, and they are complementary rather than competing [34].

The first is the delta method: differentiate the closed form with respect to every uncertain input, square, weight by each variance, and add. This is one symbolic Jacobian and it is essentially free. Its real product is not the number at the bottom, it is the ranking.

The ranking is the actionable output. Half the variance in this life estimate comes from the heat load, and a further eighteen percent from the temperature limit. The degradation time constant, the parameter that a prognostics project would naturally spend its effort on, contributes three percent. That reordering of priorities came from one differentiation of a formula, and it would have taken a large sampling campaign to discover numerically.

The second method is Monte Carlo through the same closed form, which is cheap precisely because the expression is closed form. Here it gives a standard deviation of 1,550 hours against the delta method's 1,530, so the linearization was accurate to about one percent.

And a last word on the interval itself. The nominal answer here is 2,723 hours, and the fifth to ninety-fifth percentile range is roughly 140 to 5,230 hours. That is a very wide interval on a very confident-looking number, and it is the correct answer given the stated input uncertainties. A prognostic that reports 2,723 hours without it is not more accurate. It is less honest.

Symbolic computation is older than most of the numerical tools it now feeds. It was, in fact, one of the first serious applications of artificial intelligence research.

The first milestone is worth a sentence of its own. In 1961 James Slagle's SAINT program at MIT solved freshman calculus integration problems by heuristic search, and it was presented as artificial intelligence, because that is what it was. Macsyma grew out of that line of work at MIT's Project MAC from 1968 onward, and Joel Moses, who wrote its integration engine, later published a personal history of the project [16]. Nearly every CAS alive today is downstream of Macsyma, either technically or intellectually.

For a physics based model that has to end up as running code, this ordering works and the reverse does not.

1. Write the physics as symbols, with assumptions. Declare masses positive, temperatures real, angles real. The assumptions are not decoration; they are what lets simplification proceed.
2. Non-dimensionalize before anything else. Section 16. Find out how many parameters you really have before you start sweeping seven of them.
3. Derive. Apply Euler-Lagrange, take the energy balance, apply the coordinate transform. Let the tool do the algebra you would otherwise do on paper.
4. Run the five checks. Degenerate case, equilibrium sign, conservation, dimensions, limits. Section 13. Nothing proceeds until these pass, because everything downstream inherits the error.
5. Ask what the model can and cannot tell you. Controllability, observability and structural identifiability are symbolic questions with symbolic answers, and answering them now is far cheaper than discovering the answer from a calibration that will not converge. Sections 18 and 20.
6. Simplify with intent. Use the targeted function, not the general one. Stop when the expression is small enough to read.
7. Take the derivatives you will need later. Jacobians for the estimator, gradients for the optimizer, sensitivities for identifiability. Do it now, while everything is still symbolic.
8. Choose the numerically stable form, not the shortest one. Section 17. Then discretize, remembering that the exact matrix exponential wants numbers, not symbols. Section 21.
9. Reduce and generate. cse, horner, then lambdify for Python or codegen for C. Check the operation count fell.
10. Cross the border once. From here on, arrays and floats. Do not go back into the symbolic layer inside a loop.
11. Keep the derivation script in version control next to the generated code. The generated file is an artifact, not source. When a parameter definition changes, you regenerate, you do not patch.
12. If it will be monitored, do the diagnosis work while you are still in symbols. Eliminate the unmeasured states into a redundancy relation, take the signature matrix, and find out which faults are indistinguishable before anyone specifies a sensor. Sections 36 and 37.
13. If it will be forecast, invert the damage law before you simulate it. A closed form time-to-threshold gives you the sensitivity ranking for free, which is what tells you where to spend measurement effort. Sections 40 and 42.

Symbolic computation

Manipulating mathematical expressions as structured objects, exactly, with the letters intact.

CAS

Computer algebra system. Any program whose primary data type is an expression rather than a number.

Expression tree

The internal representation of an expression: operations as nodes, operands as children.

Canonical form

A standard shape every equivalent expression is rewritten into, so equality can be tested by comparison.

Expression swell

The tendency of symbolic results to grow rapidly in size during a derivation, even when the final answer is small.

Catastrophic cancellation

Loss of all significant digits when two nearly equal floating point numbers are subtracted.

Condition number

How much a small relative change in the inputs can be amplified in the output. A property of the formula you chose, not only of the problem.

Grobner basis

A rewritten form of a system of polynomial equations that makes solving and variable elimination mechanical.

Risch algorithm

A decision procedure for whether an elementary function has an elementary antiderivative, and for finding it if so.

Buckingham Pi theorem

A relationship among n variables built from r independent dimensions can be written with exactly n minus r dimensionless groups.

Structural identifiability

Whether parameters can be recovered uniquely from perfect noise free data with the sensors you have. A property of the model, answerable before any measurement.

Differentiation index

How many times the constraints of a differential algebraic system must be differentiated before it becomes solvable as an ordinary differential equation.

Acausal modeling

Declaring physical relationships without fixing what is input and what is output, and letting the tool derive a solvable form.

Jacobian

The matrix of first partial derivatives of a vector valued function. The central object in estimation, optimization and linearization.

Lambdify

SymPy's function that compiles a symbolic expression into a plain numeric function. The exit door from the symbolic world.

Common subexpression elimination (cse)

Factoring repeated subtrees out into temporary variables before generating code.

Horner form

A polynomial written as nested multiplication, so it evaluates with only multiplies and adds and no exponentiation.

Vectorization

Applying an operation to a whole array in compiled code instead of looping in the interpreter.

Broadcasting

NumPy's rule for combining arrays of different but compatible shapes without copying data.

Automatic differentiation

Computing exact derivative values by propagating derivative information through a numeric program, without forming a symbolic formula.

Koopman operator

The linear operator describing how functions of the state, rather than the state itself, evolve. Infinite dimensional in general, but sometimes closing on a small finite set.

Observable (Koopman sense)

Any scalar function of the state. The dictionary of observables is what an EDMD model is built on, and choosing it is a modeling decision.

Analytical redundancy relation

A relation among measured signals only, obtained by eliminating the unmeasured states symbolically. Zero when healthy, non-zero when not.

Residual

The evaluated value of such a relation on real data. Detection is a test on its magnitude; isolation is a test on the pattern across several residuals.

Fault signature matrix

Which residuals respond to which faults. Two faults with identical columns cannot be told apart, whatever algorithm is applied afterward.

Damage-equivalent temperature

The single constant temperature producing the same integrated Arrhenius damage as a varying duty cycle. Always higher than the arithmetic mean.

Remaining useful life (RUL)

Time or cycles from now until a damage state reaches its threshold. An RUL quoted without an interval is not an engineering result.

Delta method

First-order uncertainty propagation: square the symbolic derivatives, weight by input variances, and add. Cheap, and it ranks the drivers.

Pareto front

The set of candidate models for which nothing simpler is also more accurate. The output of a symbolic regression search.

1. Goldberg, D. (1991). What every computer scientist should know about floating-point arithmetic. ACM Computing Surveys, 23(1), 5-48. dl.acm.org/doi/10.1145/103162.103163
2. Buchberger, B. (1965/2006). An algorithm for finding the basis elements of the residue class ring of a zero dimensional polynomial ideal. Doctoral thesis, University of Innsbruck; English translation in Journal of Symbolic Computation, 41(3-4), 475-511. sciencedirect.com/science/article/pii/S0747717105001483
3. Risch, R. H. (1969). The problem of integration in finite terms. Transactions of the American Mathematical Society, 139, 167-189. See also Risch (1970), The solution of the problem of integration in finite terms, Bulletin of the AMS, 76(3), 605-608. projecteuclid.org/euclid.bams/1183531821
4. Park, R. H. (1929). Two-reaction theory of synchronous machines: generalized method of analysis, part I. Transactions of the AIEE, 48(3), 716-727. doi:10.1109/T-AIEE.1929.5055275
5. Baydin, A. G., Pearlmutter, B. A., Radul, A. A., and Siskind, J. M. (2018). Automatic differentiation in machine learning: a survey. Journal of Machine Learning Research, 18(153), 1-43. jmlr.org/papers/v18/17-468.html
6. Wolfram, S. (2026). Launching Version 15 of Wolfram Language and Mathematica. Published 16 June 2026. writings.stephenwolfram.com
7. Meurer, A., et al. (2017). SymPy: symbolic computing in Python. PeerJ Computer Science, 3, e103. peerj.com/articles/cs-103
8. Harris, C. R., et al. (2020). Array programming with NumPy. Nature, 585, 357-362. nature.com/articles/s41586-020-2649-2
9. Virtanen, P., et al. (2020). SciPy 1.0: fundamental algorithms for scientific computing in Python. Nature Methods, 17, 261-272. nature.com/articles/s41592-019-0686-2
10. Andersson, J. A. E., Gillis, J., Horn, G., Rawlings, J. B., and Diehl, M. (2019). CasADi: a software framework for nonlinear optimization and optimal control. Mathematical Programming Computation, 11, 1-36. link.springer.com/article/10.1007/s12532-018-0139-4
11. Modelica Association. Modelica Language Specification. specification.modelica.org
12. Bradbury, J., et al. (2018). JAX: composable transformations of Python and NumPy programs. github.com/jax-ml/jax
13. Raissi, M., Perdikaris, P., and Karniadakis, G. E. (2019). Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378, 686-707. sciencedirect.com/science/article/abs/pii/S0021999118307125
14. Brunton, S. L., Proctor, J. L., and Kutz, J. N. (2016). Discovering governing equations from data by sparse identification of nonlinear dynamical systems. PNAS, 113(15), 3932-3937. pnas.org/doi/10.1073/pnas.1517384113
15. Cranmer, M., et al. (2020). Discovering symbolic models from deep learning with inductive biases. Advances in Neural Information Processing Systems 33. arxiv.org/abs/2006.11287
16. Moses, J. (2012). Macsyma: a personal history. Journal of Symbolic Computation, 47(2), 123-130. sciencedirect.com/science/article/pii/S0747717110001483
17. Slagle, J. R. (1961). A heuristic program that solves symbolic integration problems in freshman calculus. Doctoral thesis, Massachusetts Institute of Technology.
18. SymPy documentation. Simplification, assumptions, code generation, and the physics.units module. docs.sympy.org
19. NumPy documentation. Broadcasting. numpy.org/doc/stable/user/basics.broadcasting.html
20. Wolfram Research. Wolfram Language and System Documentation Center. reference.wolfram.com/language
21. IEEE Standard for Floating-Point Arithmetic, IEEE 754. The double precision format underlying every numeric result on this page, with a 53 bit significand giving roughly 15.95 decimal digits.
22. Buckingham, E. (1914). On physically similar systems: illustrations of the use of dimensional equations. Physical Review, 4(4), 345-376. link.aps.org/doi/10.1103/PhysRev.4.345
23. Ljung, L. and Glad, T. (1994). On global identifiability for arbitrary model parametrizations. Automatica, 30(2), 265-276. sciencedirect.com/science/article/abs/pii/0005109894900299
24. Pantelides, C. C. (1988). The consistent initialization of differential-algebraic systems. SIAM Journal on Scientific and Statistical Computing, 9(2), 213-231. epubs.siam.org/doi/10.1137/0909014
25. Cranmer, M. (2023). Interpretable machine learning for science with PySR and SymbolicRegression.jl. arxiv.org/abs/2305.01582
26. Koopman, B. O. (1931). Hamiltonian systems and transformation in Hilbert space. PNAS, 17(5), 315-318. pnas.org/doi/10.1073/pnas.17.5.315
27. Williams, M. O., Kevrekidis, I. G., and Rowley, C. W. (2015). A data-driven approximation of the Koopman operator: extending dynamic mode decomposition. Journal of Nonlinear Science, 25, 1307-1346. link.springer.com/article/10.1007/s00332-015-9258-5
28. Mezic, I. (2005). Spectral properties of dynamical systems, model reduction and decompositions. Nonlinear Dynamics, 41, 309-325.
29. Higham, N. J. (2002). Accuracy and Stability of Numerical Algorithms, 2nd edition. SIAM. The standard reference for conditioning and for why algebraically equal formulas are not numerically equal.
30. Mattsson, S. E. and Soderlind, G. (1993). Index reduction in differential-algebraic equations using dummy derivatives. SIAM Journal on Scientific Computing, 14(3), 677-692.
31. Gertler, J. (1998). Fault Detection and Diagnosis in Engineering Systems. Marcel Dekker. The standard treatment of analytical redundancy relations and structured residual design.
32. Blanke, M., Kinnaert, M., Lunze, J., and Staroswiecki, M. (2016). Diagnosis and Fault-Tolerant Control, 3rd edition. Springer. Structural analysis, isolability, and sensor placement worked as a systematic method.
33. Paris, P. C. and Erdogan, F. (1963). A critical analysis of crack propagation laws. Journal of Basic Engineering, 85(4), 528-534. doi:10.1115/1.3656900
34. Sankararaman, S. and Goebel, K. (2015). Uncertainty in prognostics and systems health management. International Journal of Prognostics and Health Management, 6(4). papers.phmsociety.org/index.php/ijphm/article/view/2319
35. Isermann, R. (2005). Model-based fault-detection and diagnosis: status and applications. Annual Reviews in Control, 29(1), 71-85. sciencedirect.com/science/article/abs/pii/S1367578805000052
36. Krysander, M., Aslund, J., and Nyberg, M. (2008). An efficient algorithm for finding minimal overconstrained subsystems for model-based diagnosis. IEEE Transactions on Systems, Man and Cybernetics, Part A, 38(1), 197-206. The computational core of structural fault diagnosis.
37. Archard, J. F. (1953). Contact and rubbing of flat surfaces. Journal of Applied Physics, 24(8), 981-988. The wear law used in the thrust bearing case of Section 38.
38. Coffin, L. F. (1954) and Manson, S. S. (1953), independently. The low cycle fatigue relation between plastic strain range and cycles to failure, used here in its temperature-swing form for solder joints.