Solving Differential Equations in Julia I

There are two complaints about a differential equation solver that are constantly mistaken for one another, and the whole of this book turns on keeping them apart. The first is that a tool is slow: it computes the right answer at a price one would rather not pay, and the remedy is an argument about hardware. The second is that a tool cannot state the problem at all — the model one has does not fit through the interface, at any price. The first is an engineering inconvenience. The second is a structural fact, and it is the only one of the two that settles an argument. Most of the writing in this field blurs the distinction, and a reader can smell it. I have tried to earn the claim of impossibility everywhere I make it, and to refuse it wherever only slowness is true.

The occasion for the distinction is the reference initial-value interface of scientific Python, which is genuinely excellent and which I recommend without hedging on the regime where it is enough — the small, non-stiff, first-order system given as an explicit derivative. But a discretised evolution equation does not arrive that way. It arrives as a mass matrix multiplying the derivative, with that matrix fixed by the physics — a capacitance, an inertia, the Gram matrix of a basis — and the reference interface has no argument position for it. That is not slowness. It is an empty cell in a capability matrix, and it recurs: a singular mass matrix that is really a constraint, a differential-algebraic index that must be reduced before any solver can touch the system, a second-order form that flattening doubles, a boundary condition that is a functional of the whole solution. The gap between slowness and impossibility is the subject, and the Julia SciML stack is on the page because it expresses what the reference interface cannot.

This is the first of two volumes, and the split between them is not arbitrary: the empty cells fall along two axes. This volume takes the axis of structure — the deterministic system whose form the standard interface cannot state: the mass matrix that reweights a derivative and, when singular, becomes a constraint; the differential-algebraic index that must be reduced before any solver can touch it; the symbolic-numeric layer that composes a model acausally and reduces that index automatically; the split that pays for a stiff term alone; the second-order form flattening needlessly doubles; the boundary value problem where nothing marches. The second volume takes the other axis — not the equation’s structure but its dynamics, and the workflows built over a solver: stochastic, delay, and hybrid systems; continuation and bifurcation; performance, large systems, and GPU ensembles; the differentiation of the solver itself and the universal differential equation; and the verification discipline that governs every number in both. Together they make the complete case — that across structure and across stochasticity, memory, discontinuity, scale, and differentiability, the Julia stack expresses problems the reference interface cannot, not more slowly but at all. The reader who works both learns to look at an equation, decide whether the standard tool is slow or simply cannot, and verify the solver chosen for it. The rule of evidence set here is the rule the second volume runs on, and its first chapter restates this framework before extending it.

I have held every comparison to one rule of evidence. The competing tool is given its best available configuration before a single conclusion is drawn — the right solver, an analytical Jacobian, exposed sparsity, a tuned tolerance — because a rigged benchmark is worth less than no benchmark, and the concessions are what buy the reader’s trust where the verdict is close. No timing appears without the tolerance it was run at. Where a claim is quantitative there is a work-precision diagram carrying it, because faster is not a claim; it is a diagram, or it is nothing. And no solution is reported that I have not verified against a manufactured one, watched converge at its observed order, and checked for a conserved quantity that should not drift. A solution I have not verified is a plot, not a result.