Existence and Uniqueness

Learning Objectives

By the end of this section you will be able to:

1. Diagnose whether an initial value problem has a unique solution by inspecting the right-hand side $f(x, y)$.
2. State and apply the Picard–Lindelöf theorem and the weaker Peano existence theorem, and tell which one you need.
3. Recognize Lipschitz continuity geometrically, as a finite slope cone bounding $f$ in the $y$ direction.
4. Run Picard's iteration by hand on a small problem and watch it converge to the true solution.
5. Predict how a small perturbation of the initial condition spreads through time, using the Grönwall-style bound $|\Delta y(x)| \le |\Delta y_0|\,e^{L|x - x_0|}$.
6. Identify the maximal interval of existence and recognize when a solution blows up in finite time.

The Big Picture: Why We Need Two Theorems

Up to this section we have happily written equations like $\tfrac{dy}{dx} = f(x, y),\; y(x_0) = y_0$ and gone looking for a solution. A natural and slightly alarming question is:

Does a solution exist at all? And if it does, is there only one?

These two questions are not academic. If a numerical integrator can legally land on two completely different solution curves passing through the same initial condition, then every prediction it makes is suspect. If no solution exists at all, the equation is not a model — it is a nonsense statement.

The gap between "continuous" and "Lipschitz in $y$" is small but consequential — it is exactly the gap between "solutions exist" and "the world is deterministic." The rest of this section is built around making that gap visible.

The Motivating Disaster: Many Solutions From One Point

Consider the perfectly innocent-looking initial value problem:

The right-hand side $f(x, y) = 3 y^{2/3}$ is continuous everywhere — there are no poles, no logs, no divisions. We can separate variables:

Plug in $y(0) = 0$ and we get $C = 0$, so $y(x) = x^3$ is one solution. Great — done?

Not quite. Watch:

• $y(x) \equiv 0$ is also a solution. It satisfies $y'(x) = 0 = 3\cdot 0^{2/3}$ everywhere and passes through $(0, 0)$.
• Worse, for any $a > 0$ the piecewise function $y_a(x) = 0$ for $x \le a$ and $y_a(x) = (x - a)^3$ for $x > a$ is also a perfectly valid solution.
• That gives us an entire uncountable family of solutions through the single point $(0, 0)$: the trivial zero solution sits on the x-axis as long as it likes, then leaps off into a cubic at any moment.

The root cause is local: at $y = 0$ the partial derivative $\partial f / \partial y = 2\,y^{-1/3}$ blows up to $+\infty$. The right-hand side is continuous, but it is far too steep in the $y$ direction near the x-axis. The interactive visualiser below lets you feel this failure with your own mouse.

Interactive: When Does Uniqueness Fail?

Switch between equations using the buttons above the canvas. For the cube-root and square-root cases, watch the dashed purple family of "branching" solutions passing through the highlighted yellow point. Click anywhere to launch your own numerical solution — for the Lipschitz case ($dy/dx = y$) every click follows the same exponential law; for the non-Lipschitz cases nearby clicks can split apart dramatically.

Lipschitz Continuity: The Fix

The mathematical patch that pins down a unique solution is a mild but crucial smoothness assumption on $f$ in the $y$ direction.

What Lipschitz Means Geometrically

Lipschitz means that the secant slope between two points stacked vertically can never be larger than $L$. Equivalently — if $f$ is differentiable in $y$ — it just means $|\partial f/\partial y| \le L$.

Interactive: The Lipschitz Cone

For each $f$ in the picker, the visualiser below draws the steepest pair of lines that $f$ could possibly produce on a yellow disk around the clicked point. If those cyan lines stay finite as you slide the radius, $f$ is Lipschitz there. If the local $L$ shoots to infinity as the disk crosses a bad point, you have just witnessed uniqueness failing.

The Picard–Lindelöf Theorem

Three pieces of information come bundled in the theorem:

1. Existence. The solution exists.
2. Uniqueness. Two solutions of the same IVP must agree on the entire interval $[x_0 - h, x_0 + h]$.
3. Length estimate. The interval has length $2h$, where $h \le b/M$. Steep right-hand sides (large $M$) shrink the guaranteed interval.

Picard Iteration: Solving by Successive Approximation

Picard's proof is constructive — it tells you how to build the unique solution, not just that it exists. The trick is to convert the differential equation into an integral equation:

This identity holds for any candidate solution. It motivates the recursion:

On the interval where Picard–Lindelöf applies, this sequence converges uniformly to the (unique) solution. The mechanism is a contraction: the operator that maps $y_k \mapsto y_{k+1}$ shrinks the gap between any two candidates by a factor of $Lh$ per step, and $h$ is chosen small enough that $Lh < 1$.

Interactive: Watching Picard Converge

For the canonical example $dy/dx = y,\; y(0) = 1$, the recursion produces $y_n(x) = \sum_{k=0}^{n} x^k/k!$ — the partial sums of the Taylor series for $e^x$. Slide the iteration index to watch the cyan curve climb to the dashed white truth line.

Worked Example by Hand

Apply Picard's recursion to $dy/dx = y,\; y(0) = 1$ step by step. The Lipschitz constant on any rectangle is $L = 1$, so the Picard–Lindelöf hypothesis holds — uniqueness is guaranteed, and our successive iterates must all converge to the same limit.

Peano's Theorem: Existence Without Uniqueness

Peano's theorem is what saves the cube-root example from being nonsense. There, $f(x, y) = 3 y^{2/3}$ is continuous everywhere, so a solution exists through every initial point. Peano simply does not promise it is unique.

Continuous Dependence on Initial Conditions

Uniqueness is great in theory, but in practice we never know $y_0$ exactly. Measurement noise, rounding, and missing decimals all wiggle it. The same proof that gives uniqueness also tells us how much our solution wiggles back:

This is sometimes called the Grönwall estimate. The Lipschitz constant $L$ controls how fast nearby trajectories can diverge. For $dy/dx = y$ with $L = 1$, two trajectories starting $10^{-6}$ apart can be a factor of $e^{10} \approx 22\,000$ apart after only ten time units. This is the inequality behind every numerical error bound you ever encounter in scientific computing.

The Maximal Interval of Existence

Picard–Lindelöf is a local theorem. It guarantees a unique solution on a possibly tiny interval. The natural question is: how far can we extend it?

Solving by separation: $y(x) = 1/(1 - x)$. The right-hand side is smooth — but the solution blows up at x = 1. We say the solution escapes in finite time.

The maximal interval of existence $(\alpha, \omega)$ is the largest open interval on which the unique solution lives. Picard–Lindelöf can be iterated repeatedly to extend the local interval as long as $y(x)$ stays in the region where $f$ is Lipschitz. The solution can only stop existing for two reasons:

1. Blow-up: $|y(x)| \to \infty$ as $x \to \omega^{-}$.
2. Boundary escape: the trajectory reaches the boundary of the domain where $f$ is defined.

Python: Picard Iteration From Scratch

Let's implement Picard's recursion in plain Python so you can see every quadrature node and every iterate. We'll apply it to $dy/dx = y, \; y(0) = 1$ — the same example we worked by hand — and print the convergence table.

picard_iterate(f=λ x,y: y, x0=0, y0=1, x_target=1, n_iters=6)

fl = f(0.000, y_prev[0]) = 1.000   ;   fr = f(0.005, y_prev[1]) ≈ 1.005

y_0(1) = 1.000000  |y_0(1) - e| = 1.7183e+00
y_1(1) = 2.000000  |y_1(1) - e| = 7.1828e-01
y_2(1) = 2.500000  |y_2(1) - e| = 2.1828e-01
y_3(1) = 2.666667  |y_3(1) - e| = 5.1615e-02
y_4(1) = 2.708333  |y_4(1) - e| = 9.9485e-03
y_5(1) = 2.716667  |y_5(1) - e| = 1.6152e-03
y_6(1) = 2.718056  |y_6(1) - e| = 2.2607e-04
e        = 2.718282

1import math
2from typing import Callable, List
3
4def picard_iterate(
5    f: Callable[[float, float], float],
6    x0: float,
7    y0: float,
8    x_target: float,
9    n_iters: int,
10    n_quad: int = 200,
11) -> List[float]:
12    """
13    Compute Picard's successive approximations to the IVP
14        dy/dx = f(x, y),    y(x0) = y0.
15
16    We build y_{k+1}(x) = y0 + ∫_{x0}^{x} f(t, y_k(t)) dt
17    by storing each y_k as a dense table of (t, value) samples.
18    """
19    # Discretise the interval [x0, x_target] with n_quad subintervals.
20    ts = [x0 + (x_target - x0) * i / n_quad for i in range(n_quad + 1)]
21    h = (x_target - x0) / n_quad
22
23    # y_0(t) ≡ y0 — the constant initial guess.
24    y_prev = [y0 for _ in ts]
25    y_values_at_target: List[float] = [y_prev[-1]]
26
27    for k in range(n_iters):
28        y_next = [y0]                       # y_{k+1}(x0) = y0 always
29        running_integral = 0.0
30        for i in range(n_quad):
31            # Trapezoidal rule on [t_i, t_{i+1}]
32            fl = f(ts[i], y_prev[i])
33            fr = f(ts[i + 1], y_prev[i + 1])
34            running_integral += 0.5 * h * (fl + fr)
35            y_next.append(y0 + running_integral)
36
37        y_prev = y_next
38        y_values_at_target.append(y_prev[-1])
39
40    return y_values_at_target
41
42
43# Demo: dy/dx = y, y(0) = 1.  True solution is y = e^x.
44def f(_x: float, y: float) -> float:
45    return y
46
47approx = picard_iterate(f, x0=0.0, y0=1.0, x_target=1.0, n_iters=6)
48truth = math.exp(1.0)
49
50for k, val in enumerate(approx):
51    err = abs(val - truth)
52    print(f"y_{k}(1) = {val:.6f}    |y_{k}(1) - e| = {err:.4e}")
53print(f"e        = {truth:.6f}")

Run this and you reproduce the table from the worked example to six decimal places. The convergence is roughly factorial: each step kills about one more decimal because the next dropped Taylor term is $x^{n+1}/(n+1)!$.

PyTorch: Sensitivity via Autograd

The Grönwall estimate $|\Delta y(x)| \le |\Delta y_0|\, e^{L|x - x_0|}$ is an inequality. For $dy/dx = y$, equality is achieved — and we can verify it experimentally by asking PyTorch to differentiate the numerical solution with respect to the initial condition.

y0 = tensor(1.0, requires_grad=True)

y(T)              = 2.718279   (truth e^T = 2.718282)
dy(T)/dy0         = 2.718279   (Picard–Lindelöf prediction e^T)
Lipschitz bound L = 1   ⇒   |y(T)-y(T;y0+δ)| ≤ δ · e^(L·T) = δ · e ≈ 2.7183·δ

1import torch
2
3# Continuous dependence: how much does y(T) move if we wiggle y0?
4# We solve dy/dx = y from x = 0 to x = T with RK4 and let PyTorch's
5# autograd compute  dy(T)/dy0  directly.
6
7def rk4_solve(f, y0: torch.Tensor, T: float, n_steps: int) -> torch.Tensor:
8    """Forward Euler? No — that's noisy. RK4 is differentiable too."""
9    h = T / n_steps
10    x = torch.tensor(0.0)
11    y = y0
12    for _ in range(n_steps):
13        k1 = f(x, y)
14        k2 = f(x + h / 2, y + h * k1 / 2)
15        k3 = f(x + h / 2, y + h * k2 / 2)
16        k4 = f(x + h, y + h * k3)
17        y = y + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6
18        x = x + h
19    return y
20
21
22def f(_x: torch.Tensor, y: torch.Tensor) -> torch.Tensor:
23    return y                                # dy/dx = y
24
25
26y0 = torch.tensor(1.0, requires_grad=True)  # ← variable we differentiate w.r.t.
27T = 1.0
28
29yT = rk4_solve(f, y0, T, n_steps=200)       # y(T) ≈ e^T
30
31# PyTorch tracks every op above and now computes dy(T)/dy0 with one call.
32sensitivity = torch.autograd.grad(yT, y0)[0]
33
34print(f"y(T)              = {yT.item():.6f}   (truth e^T = {2.718282:.6f})")
35print(f"dy(T)/dy0         = {sensitivity.item():.6f}   (Picard–Lindelöf prediction e^T)")
36print(f"Lipschitz bound L = 1   ⇒   |y(T)-y(T;y0+δ)| ≤ δ · e^(L·T) = δ · e ≈ {2.718282:.4f}·δ")

The autograd answer $dy(T)/dy_0 \approx e^T$ reproduces the Grönwall bound exactly. This is precisely how neural-ODE training propagates gradients — the only difference is that $f$ is a neural network instead of an elementary function.

Why This Matters in Practice

Common Pitfalls

Summary

Key Takeaways

1. Existence and uniqueness are separate guarantees. Peano delivers the first; Picard–Lindelöf delivers both — and the second is what makes a differential equation a predictive model.
2. The decisive hypothesis is Lipschitz in y, which really means "the partial derivative $\partial f/\partial y$ is locally bounded." Smoothness in $x$ is not enough.
3. Picard's constructive proof is also a numerical algorithm. For $dy/dx = y, y(0) = 1$ the iterates are exactly the Taylor partial sums of $e^x$.
4. The same Lipschitz constant $L$ that guarantees uniqueness also bounds how fast nearby trajectories can spread. This is the Grönwall estimate and it is the root of every stability bound in numerical ODE solving and chaos theory.
5. Smooth right-hand sides can still produce solutions with finite maximal intervals. Uniqueness is local; always check whether your solution survives to the time you care about.

Coming Next: Chapter 21 begins our toolbox of explicit solution techniques. We'll start with Linear First-Order Equations, where the integrating factor method turns every such ODE into a single antiderivative — cleanly, and with uniqueness automatically guaranteed.