Existence, uniqueness and stepping forward

Two lessons of methods have solved every separable and every first order linear equation, and almost nothing else, so the natural next question is what can be said about an equation that no method touches.

There are two questions, and they are independent. Does the initial value problem $y' = f (t, y)$ , $y (t_{0}) = y_{0}$ have a solution at all? And if it has one, does it have only one? Both matter practically. If a model has two solutions from the same starting state, the model is not determining the future, whatever the physics claimed. If a solution stops existing at some finite time, a computer will happily keep printing numbers past that point, and they will be nonsense.

Picard's theorem

The result is due to Émile Picard and Ernst Lindelöf in the 1890s, and its statement is short. Suppose $f (t, y)$ is continuous in a rectangle around the point $(t_{0}, y_{0})$ , and suppose $\partial f / \partial y$ is continuous there too. Then the initial value problem has exactly one solution, defined on some interval around $t_{0}$ .

Three things about that statement are worth pressing on.

Continuity of $f$ alone gives existence. That is Peano's theorem, and it is genuinely weaker: it gives a solution but permits several.

Uniqueness needs the extra condition on $\partial f / \partial y$ , and what is really required is a Lipschitz condition, that $| f (t, y_{1}) - f (t, y_{2}) | \leq L | y_{1} - y_{2} |$ for some constant $L$ . A bounded $\partial f / \partial y$ supplies one by the mean value theorem, which is why the derivative version is the one quoted. The condition says the equation's right side cannot respond infinitely fast to a change in $y$ .

The conclusion is local. There is an interval, and the theorem does not say how long it is. That is not slackness in the proof; solutions really do end, as the third section shows.

Example. Does $y' = y^{2} + t^{2}$ with $y (0) = 1$ have exactly one solution near $t = 0$ ?

Here $f = y^{2} + t^{2}$ is continuous everywhere and $\partial f / \partial y = 2 y$ is continuous everywhere, so both hypotheses hold in any rectangle around $(0, 1)$ . There is exactly one solution on some interval around zero. The theorem promises nothing beyond that interval, and in this case the promise would be false if it did: this solution blows up before $t = 1$ .

Now you. Consider $y' = \sqrt{y}$ with $y (0) = 0$ , restricted to $y \geq 0$ . Which hypothesis holds and which fails?

Answer

$f = \sqrt{y}$ is continuous, so a solution exists. But $\partial f / \partial y = 1 / (2 \sqrt{y})$ is unbounded as $y \to 0$ , so the Lipschitz condition fails exactly at the initial point and uniqueness is not guaranteed.

The bucket that cannot say when it started draining

That second example is not a technicality. Solve it: separating, $\int y^{- 1 / 2} d y = \int d t$ gives $2 \sqrt{y} = t + C$ , so $y = (t + C)^{2} / 4$ . With $y (0) = 0$ that is $y = t^{2} / 4$ . But $y \equiv 0$ is also a solution, the constant one that separation discarded. And so is every function that sits at zero until some time $a$ and then leaves along $y = (t - a)^{2} / 4$ : it is differentiable at the join, since both pieces have slope zero there, and it satisfies the equation on both sides. There are infinitely many solutions through the origin.

The physical version is a leaking tank. Torricelli's law says water drains from a hole at a rate proportional to the square root of the depth, so $h' = - c \sqrt{h}$ , which is the same equation with the sign reversed. Run it forwards and the tank empties at a definite time and stays empty. Run the film backwards and ask when a currently empty tank was full: the equation cannot say. Every emptying time is consistent with what you see now. Non-uniqueness here is a real property of the model, and it happens precisely because the rate becomes insensitive to the state, with infinite slope in $h$ , as the depth reaches zero.

Linear equations never do this. For $y' + p (t) y = g (t)$ the right side is $g - p y$ , whose $y$ -derivative is $- p (t)$ , continuous wherever the coefficients are. So a linear initial value problem has exactly one solution on the whole interval where $p$ and $g$ are continuous, with no local hedge. That is the real reward for linearity.

Blow-up in finite time

The other failure is more common in practice. Take $y' = y^{2}$ with $y (0) = y_{0} > 0$ . Separating, $- 1 / y = t + C$ , so

y = \frac{y_{0}}{1 - y_{0} t}

which is finite and smooth up to $t = 1 / y_{0}$ and infinite there. Nothing in the equation is singular at that time; the coefficient functions are polynomials. The solution manufactures its own catastrophe by growing fast enough to feed its own growth, and where the catastrophe happens depends on the initial condition.

Example. For $y' = y^{2}$ with $y (0) = 0.5$ , when does the solution cease to exist, and what is $y$ at $t = 1$ ?

The formula gives $y = 0.5 / (1 - 0.5 t) = 1 / (2 - t)$ , which is infinite at $t = 2$ . At $t = 1$ , $y = 1$ . Halving the starting value doubled the survival time.

Now you. For $y' = y^{2}$ with $y (0) = 4$ , when does the solution cease to exist, and what is $y$ at $t = 0.2$ ?

Answer

$y = 4 / (1 - 4 t) = 1 / (0.25 - t)$ , so it blows up at $t = 0.25$ . At $t = 0.2$ , $y = 1 / 0.05 = 20$ .

The same happens without separability. $y' = 1 + y^{2}$ with $y (0) = 0$ has solution $y = \tan t$ , which ends at $t = π / 2 = 1.571$ . A superlinear right side is the warning sign, and physical models that produce one, such as thermal runaway or an autocatalytic reaction, are describing systems that genuinely go somewhere the model stops being valid.

Picard's iteration, which is also the proof

The theorem is proved constructively, and the construction is worth seeing because it turns the differential equation into something a computer could run. Integrate both sides of $y' = f (t, y)$ from $t_{0}$ and use the initial condition:

y (t) = y_{0} + \int_{t_{0}}^{t} f (s, y (s)) d s

This integral equation is equivalent to the initial value problem, and it has $y$ on both sides, which invites iteration. Start with the constant guess $y_{0} (t) = y_{0}$ and define each next approximation by substituting the previous one into the right side.

For $y' = y$ with $y (0) = 1$ the iterates are immediate. From $y_{0} = 1$ , the first is $1 + \int_{0}^{t} 1 d s = 1 + t$ . The next is $1 + \int_{0}^{t} (1 + s) d s = 1 + t + t^{2} / 2$ . The next adds $t^{3} / 6$ . The iterates are the partial sums of the Taylor series of $e^{t}$ , converging to it, which is the right answer. The Lipschitz condition is exactly what makes each iteration shrink the error by a fixed factor, so the sequence converges, and it is what fails for $\sqrt{y}$ .

Picard iteration is rarely used for computation, since each step needs a symbolic integral. What is used is a cruder relative of the same idea.

Euler's method

Leonhard Euler's method, from 1768, is the tangent line taken seriously. At $(t_{n}, y_{n})$ the equation gives the slope $f (t_{n}, y_{n})$ directly, so follow it for a small step $h$ :

y_{n + 1} = y_{n} + h f (t_{n}, y_{n}), t_{n + 1} = t_{n} + h

and repeat. No solving is involved at any point, and the method needs nothing from $f$ except the ability to evaluate it.

Example. Take $y' = t + y$ with $y (0) = 1$ and $h = 0.1$ . Do two steps, and compare with the exact solution $y = 2 e^{t} - t - 1$ .

The first step uses the slope at $(0, 1)$ , which is $0 + 1 = 1$ , giving $y_{1} = 1 + 0.1 \times 1 = 1.1$ at $t = 0.1$ . The second uses the slope at $(0.1, 1.1)$ , which is $1.2$ , giving $y_{2} = 1.1 + 0.1 \times 1.2 = 1.22$ at $t = 0.2$ . The exact value is $2 e^{0.2} - 1.2 = 1.2428$ , so two steps have accumulated an error of $0.023$ , about $1.8$ per cent, and it is an underestimate because the solution is convex and the tangent lines all lie below it.

Now you. Take $y' = y - t^{2}$ with $y (0) = 1$ and $h = 0.1$ , and do two steps. The exact value at $t = 0.2$ is $1.2186$ .

Answer

Slope at $(0, 1)$ is $1$ , so $y_{1} = 1.1$ . Slope at $(0.1, 1.1)$ is $1.1 - 0.01 = 1.09$ , so $y_{2} = 1.1 + 0.109 = 1.209$ . The error is $0.0096$ .

How wrong is it in general? One step's error is the remainder of a Taylor expansion truncated after the linear term, so it is proportional to $h^{2}$ . But reaching a fixed time takes $1 / h$ steps, so the errors accumulate to something proportional to $h$ . Euler's method is first order: halve the step and you halve the error.

That is easy to see numerically. Integrating $y' = y$ from $y (0) = 1$ to $t = 1$ , where the exact answer is $e = 2.71828$ , Euler with $h = 0.25$ gives ${1.25}^{4} = 2.4414$ , an error of $0.277$ . With $h = 0.125$ it gives $2.5658$ , error $0.153$ . With $h = 0.0625$ , error $0.080$ . Each halving of the step cuts the error by a factor approaching two, exactly as the analysis says. To get three decimal places this way needs tens of thousands of steps.

Runge-Kutta, and why it is worth the extra evaluations

Euler uses the slope at the start of the interval to cross the whole interval, which is systematically wrong whenever the slope changes. The improvements all sample the slope more than once.

The improved Euler or Heun method takes a trial Euler step, evaluates the slope at the far end, and steps with the average of the two slopes. The trial step costs one extra evaluation of $f$ and buys an order: the error goes as $h^{2}$ . On the same test problem it gives $2.6949$ at $h = 0.25$ , error $0.0234$ , and $2.7118$ at $h = 0.125$ , error $0.0064$ . The ratio is close to four, as second order predicts.

The classical fourth order Runge-Kutta method, published by Wilhelm Kutta in 1901, samples four slopes per step: one at the start, two at the midpoint, one at the end, and combines them with weights $1, 2, 2, 1$ over $6$ . Its error goes as $h^{4}$ . On the same problem with $h = 0.25$ , four steps and sixteen evaluations of $f$ give $2.718210$ , an error of $7 \times 10^{- 5}$ ; at $h = 0.125$ the error is $5 \times 10^{- 6}$ . Euler would need around twenty thousand steps to match the first of those. This is why almost every general purpose solver is a Runge-Kutta method or a relative, usually with the step size adjusted automatically by comparing two estimates of different order.

Where the numbers lie

A numerical solution is a sequence of numbers, and numbers always appear, whether or not a solution exists.

Run a solver on $y' = y^{2}$ with $y (0) = 1$ past $t = 1$ and it will keep producing values. They are meaningless: the solution ceased to exist at $t = 1$ . A solver that suddenly needs tiny steps to keep its error estimate down is usually reporting a blow-up rather than a difficulty of its own, and the correct response is to ask what the model is doing, not to lower the tolerance.

The second trap is stiffness. Consider $y' = - 1000 (y - \cos t)$ . Its solution settles onto something close to $\cos t$ within about three milliseconds and then varies on a timescale of seconds, so a step size of a hundredth of a second should be ample for accuracy. It is catastrophic for stability: Euler's method applied to $y' = λ y$ multiplies the error by $| 1 + h λ |$ each step, which exceeds one unless $h < 2 / 1000$ . Take a larger step and the computed solution oscillates with growing amplitude and diverges, while the true solution sits quietly on a cosine. Stiff problems, which are the norm in chemical kinetics and circuit simulation, need implicit methods that solve for $y_{n + 1}$ rather than stepping to it.

The third is more mundane: halving the step doubles the number of arithmetic operations, so rounding error accumulates as the truncation error falls, and below some step size the total error starts to rise again.

None of this means numerics is second best. It means a number from a solver is a claim that needs the same scepticism as any other measurement, and the way to check it is to halve the step and see whether the answer moves by what the method's order predicts.

The tools so far give a formula when one exists and a table when it does not. Neither says anything about a solution's shape without computing it. The next lesson shows that for a large class of equations, the long-term behaviour can be read straight off the right side, with no solving and no stepping at all.