Stability and the CFL condition

In the previous lesson, we studied the numerical solution of the linear and non-linear convection equations, using the finite-difference method. We began by discretizing the one-way wave equation using a classic forward-time/backward space scheme, and computing the solution using an initial condition consisting of a square pulse.

Computing with finer spatial grids while keeping the time step constant, you encountered four cases:

1. a coarse run with nx=41 smoothed the square pulse and attenuated its height,

2. a finer run with nx=81 improved the solution, but smoothing still ocurred,

3. a surprising case with nx=101matched the exact solution, and

4. a further refinement with nx=121 destroyed the solution.

In this lesson, we will explore why changing the discretization parameters can affect your solution in such a drastic way. The central question we want to answer is:

Why did refining the spatial grid improve the linear-convection calculation, make it an exact grid shift at one setting, and then destroy it?

Let’s begin by importing our favorite Python libraries for numerical computing.

The code below corresponds to the same function we used in the previous lesson, but incorporating the vectorized update in space, leaving only the loop in time. Be sure to review the way array slicing works in this code, sketched in Figure 5 of Lesson 6, and the Python refresher that came after.

How far does the wave travel in one step?¶

Look again at the discretized equation implemented by our function:

Pay attention to the coefficient multiplying the spatial difference: c * dt / dx. The speed $c$ has units of distance per time, so $c\Delta t$ is the distance the exact wave travels during one time step. Dividing by the grid spacing expresses that distance in grid intervals. We give this dimensionless ratio a name:

This is the Courant number, also known as the CFL number (Courant–Friedrichs–Lewy) Courant et al., 19281967. For the positive speed considered here, $C=0.8$ means that the wave travels eight-tenths of a grid interval per time step; $C=1$ means one complete interval.

The ratio is already present in Equation 1 and in our code. Refining the spatial grid while keeping the time step fixed makes $\Delta x$ smaller and the Courant number larger. What values did it take in the four runs from the previous lesson?

The four Courant numbers are 0.4, 0.8, 1.0, and 1.2. The run that shifted the pulse exactly had $C=1$; the failed run had $C=1.2$. These observations suggest that this ratio matters, but they do not yet explain the failure or establish a stability condition.

To investigate, we will ask a more focused question: what happens to a small disturbance in the initial data?

Perturbation study¶

Imagine running the same calculation twice, with just one initial value changed by 10^-3 in the second run. Both runs use the same grid, time step, method, and inflow value. Does their difference remain small as we advance in time, or does the calculation amplify it?

Use a constant background, $u=1$, so that the disturbance is easy to isolate. In the second array, add 10^-3 at $x=0.5$. The unperturbed exact solution stays constant; the equation transports a disturbance without increasing its height. Our numerical update may behave differently.

We will repeat this paired experiment on nx = 81, 101, 121, keeping $c=1$ and $\Delta t=0.02$. These are the three runs with $C=0.8$, 1.0, and 1.2. Advance each pair for 25 steps, reaching the same physical time $T=0.5$ on every grid.

At each step, measure the largest absolute difference between the two arrays:

Here, $a$ denotes the unperturbed run and $b$ the perturbed run. Initially, $E^0=10^{-3}$. Taking absolute values matters because a disturbance can develop both positive and negative differences.

The perturbation occupies one grid point, so its physical width changes with the grid. We are testing how the method amplifies a grid disturbance; this is not a convergence comparison for one fixed continuous initial profile.

Each run starts from fresh arrays. The point at nx // 4 lies at $x=0.5$ on all three grids; // is integer division. The left boundary remains $u=1$ in both runs, because the function updates only u[1:].

We call advance_linear_convection() with num_steps=1 to inspect the difference after every step. This uses the same update as advancing all 25 steps at once, while letting us record its behavior along the way. The FTBS stencil can pass a disturbance at most one grid point to the right per step. Even on the coarsest grid, the perturbed point is more than 25 grid intervals from the outflow, so the disturbance cannot leave the array during this experiment.

Before running: predict which of the three cases will reduce, preserve, or amplify the maximum difference. Treat the prediction as a hypothesis to test.

Let’s plot the recorded differences. A logarithmic vertical scale lets us see shrinking and growing disturbances on the same axes. Equal vertical distances represent equal multiplicative changes, rather than equal additive changes.

A few selected steps make the amount of amplification easier to compare:

For $C=0.8$, the maximum difference decreases as the disturbance spreads over neighboring points. For $C=1$, its maximum stays at 10^-3. For $C=1.2$, a difference initially one-thousandth of the background grows to about 1 after only 25 steps. The calculation has turned a small change in the initial data into a large change in the result.

This is the question behind numerical stability: does the method control the amplification of small disturbances? More precisely, over a fixed physical time interval, we seek a bound on amplification that remains independent of the mesh as we refine it under the stated conditions. Stability need not mean that every disturbance decays; the $C=1$ result illustrates controlled propagation without decay.

These three runs provide evidence for particular choices of grid and time step. They do not prove a general bound. We still need to explain why the coefficient in Equation 1 separates these behaviors.

Where does the information come from?¶

We introduced the Courant number by asking how far the exact wave travels during one time step. Now turn that question around: where was the information needed at $x_i$ one time step earlier?

For constant positive speed $c$, the solution travels along a characteristic without changing its value. Tracing that characteristic backward from $(x_i,t_{n+1})$ gives

The exact solution therefore needs the value at $x_i-c\Delta t$. Our numerical update, Equation 1, uses only two old values: $u_{i-1}^n$ and $u_i^n$. Figure 1 puts these two descriptions on the same space–time diagram.

In the case drawn, $0<C<1$: the characteristic lands between $x_{i-1}$ and $x_i$. At $C=1$, it lands exactly on $x_{i-1}$. For $C>1$, it lands farther to the left, outside the interval spanned by the two input points.

This is the local picture of a domain of dependence: the earlier data that can influence a particular solution value. To see its meaning over several steps, trace the stencil backward again. Two steps earlier, the FTBS value can depend on points $i$, $i-1$, and $i-2$; after $m$ steps, its numerical dependence extends at most from $x_{i-m}$ to $x_i$. The exact characteristic reaches $x_i-mc\Delta t$ over the same interval. We are considering points whose backward stencils remain inside the spatial domain; at an inflow boundary, the prescribed boundary data also enter the dependence.

For the characteristic to stay within the numerical dependence interval, we need

We have assumed positive speed and positive step sizes. The limiting case $C=0$ corresponds to no travel during a step. Notice that equality is allowed: traveling exactly one grid interval does not put the required information beyond the stencil.

If $C>1$, the exact wave travels farther during a fixed time interval than information can pass through the numerical stencil. Refining with that same supercritical Courant number does not remove the mismatch. The numerical method cannot, in general, converge to a solution that depends on data it cannot reach. This information-containment requirement is the CFL condition for the stencil considered here.

The geometry explains why $C=1$ is a meaningful threshold, but it does not tell us how the update combines the available values. Containing the characteristic is a necessary condition for convergence of these explicit transport schemes; it is not, by itself, a proof that disturbances remain controlled. To explain the decay, preservation, and growth measured in our perturbation study, we will next examine the weights multiplying the two old values.

How the update controls a disturbance¶

We can establish a stability bound directly from the update. First, collect the terms multiplying each old value in Equation 1:

When $0\leq C\leq1$, the two weights, $1-C$ and $C$, are nonnegative and sum to one. The new value is therefore a weighted average of the two old values: it lies between them. Such an average is called a convex combination.

To explain the perturbation experiment, we need to apply this observation to the difference between the two runs.

Subtract the two updates¶

Let $u_{a,i}^n$ and $u_{b,i}^n$ denote the unperturbed and perturbed runs, and write their difference as

Here, $e$ is the difference between two numerical solutions, not the error relative to an exact solution. Both runs use the same $c$, $\Delta x$, and $\Delta t$, so their updates have the same weights:

Subtract the first update from the second and collect the differences:

The difference obeys the same update as the solution. This follows from the linearity of this scheme with constant speed; we cannot assume it for a nonlinear update.

Bound the largest difference¶

The quantity measured in our experiment was

Every old difference has magnitude at most $E^n$. To bound the new difference, use the triangle inequality: the magnitude of a sum is no larger than the sum of the magnitudes. For $0\leq C\leq1$, the weights are nonnegative, giving

This bounds every updated point, including the rightmost point. At the left boundary, both runs keep the same prescribed value, so their difference is zero. Taking the maximum over the whole array therefore gives $E^{n+1}\leq E^n$. Repeating the argument over successive steps yields

The largest initial disturbance cannot be amplified. This bound holds for any initial difference, not just the single-point disturbance we tested. Its amplification bound is one, independent of the grid spacing and number of steps: it establishes stability for the stated scheme and boundary treatment.

Read the experiment again¶

The decrease observed at $C=0.8$ is consistent with Equation 12. The bound permits this decrease, but does not require every disturbance to decay.

At $C=1$, Equation 9 becomes $e_i^{n+1}=e_{i-1}^n$: the disturbance shifts one grid point to the right at each step. Its maximum stays unchanged while it remains in the domain, exactly as we measured. The same shift explains the square-pulse result from Lesson 6. Equality belongs in the stability condition.

At $C=1.2$, the weight $1-C$ is negative, so the nonnegative-weight argument no longer applies. That alone does not prove growth. Our experiment has exhibited a growing disturbance; an analytical explanation of that growth is the next question.

A disturbance that grows when $C>1$¶

To prove that the method can amplify a disturbance, we only need to find one that grows. Choose a perturbation whose sign alternates at neighboring grid points:

For even $i$, $(-1)^i=1$; for odd $i$, $(-1)^i=-1$. The initial differences therefore look like $\varepsilon,-\varepsilon,\varepsilon,-\varepsilon,\ldots$, all with magnitude $\varepsilon$. We can make $\varepsilon$ as small as we like.

This is a different disturbance from the single-point perturbation in our experiment. We choose it because adjacent differences are exact opposites, which makes its evolution easy to calculate. For now, consider an alternating pattern extending over the grid without boundaries; we will return to the finite interval below.

Follow one update¶

Suppose the difference at step $n$ has the form

where $A_n$ is its signed amplitude and $A_0=\varepsilon$. Its left neighbor has the opposite sign:

Substitute both expressions into Equation 9:

The alternating pattern is preserved. Only its amplitude changes, according to

We call $G$ the amplification factor for this pattern. Its sign tells us whether the signs flip at each step; its magnitude tells us whether the disturbance shrinks, stays the same size, or grows.

Repeat the update¶

Starting from $A_0=\varepsilon$, successive steps give

Continuing in the same way,

Now let $C>1$. Then $1-2C<-1$, so

The magnitude increases at every step, and the signs reverse. For example, at $C=1.2$, $G=-1.4$: each step multiplies the magnitude by 1.4. After 20 steps, it is about 837 times its initial value. Even an arbitrarily small initial disturbance can be amplified by a large factor.

At $C=1$, $G=-1$: the signs flip but the magnitude stays unchanged, consistent with shifting an alternating pattern by one grid point. At $C=0.8$, $G=-0.6$, so its magnitude decreases by a factor of 0.6 per step. These factors describe this alternating pattern; they do not give the maximum difference for the earlier single-point disturbance.

What about the boundary and a fixed final time?¶

Our numerical runs have a prescribed left boundary, with zero difference between the two inflow values. We cannot impose the alternating pattern at that boundary as well. Instead, set $e_0^0=0$ and alternate the initial differences at the interior points.

After $n$ steps, the backward stencil of point $i$ reaches only as far as $i-n$. Thus, at points more than $n$ grid intervals from the left boundary ($i>n$), all the initial differences contributing to the update belong to the alternating pattern. Equation 19 holds there exactly. A numerical check of this prediction must use such points, rather than a region already influenced by the inflow.

To connect growth to our definition of stability, hold $C>1$ and a physical final time $T$ fixed. With $c>0$,

using grids for which $T$ is an integer number of steps. As $\Delta x$ decreases, $n$ increases. The amplification $(2C-1)^n$ therefore grows without bound. On our finite domain, choose a sufficiently short fixed $T$ so that an interior point satisfies $x_i>n\Delta x=cT/C$ throughout refinement. The growing pattern at that point then remains unaffected by the inflow. Since the maximum difference is at least as large as the difference at that point, there can be no mesh-independent amplification bound.

Together, the two arguments establish the sharp condition $0\leq C\leq1$ for this FTBS scheme with the boundary treatment considered: within that interval, no initial difference is amplified in maximum magnitude; for $C>1$, an alternating initial difference provides a counterexample. This does not mean every initial state grows when $C>1$—a constant state, for example, remains constant. Stability must control all admissible small disturbances.

Make a robust solver with an agent¶

Our short function makes FTBS easy to inspect, but it leaves its assumptions to the user (a.k.a. caller). We now know enough to build a solver that enforces those assumptions, chooses a stable time step, and reaches a requested final time. An agent can draft the validation and tests while the numerical reasoning is fresh in our minds.

The activity below supplies the design. Your job is to turn it into a clear specification, inspect the agent’s implementation, and establish whether its tests provide useful evidence. The result will solve constant-speed transport on a uniform grid with constant inflow; it will not decide whether that model or grid is appropriate for an engineering application.

1. Diagnose four weaknesses¶

The following calls all return arrays without raising an exception. Use a small initial profile so that you can inspect the results directly. Its grid is $x=[0,1,2,3]$, so the spacing is 1.

With the intended positive-speed update at $\Delta t=0.2$, the result would be [1.0, 1.8, 1.2, 1.0]. Compare that reference with the outputs and use the table to explain each weakness.

A visible failure can still justify a check before computation. For the last case, another conditional in the old function cannot recover missing information: we need a better interface.

2. Derive the numerical requirements¶

For $c<0$, information arrives from the right. Replace the backward difference with a forward difference:

Here $C=c\Delta t/\Delta x$ is signed. The weights $1+C$ and $-C$ are nonnegative and sum to one when $-1\leq C\leq0$. Together with the positive-speed result, this gives $|C|\leq1$, provided the spatial difference uses the upstream neighbor. Hold the left endpoint fixed for $c>0$ and the right endpoint fixed for $c<0$. Because the speed is constant, choose the direction once per solver call. For $c=0$, the state does not change.

Next, let $q$ be a target Courant magnitude, with $0<q\leq1$. For nonzero speed and requested duration $T>0$, choose

The ceiling rounds up to an integer number of steps. Consequently $\Delta t\leq\Delta t_{\mathrm{limit}}$, while $N\Delta t=T$ in exact arithmetic. Report the realized signed $C=c\Delta t/\Delta x$, whose magnitude may be smaller than $q$. Handle $c=0$ or $T=0$ separately, before any division that assumes nonzero values.

3. Use this design brief¶

The new interface is

solve_transport(u0, x_limits, c, t_final, courant=0.8, max_steps=100_000)

u0 represents values at equally spaced points including both domain endpoints. The solver constructs x from x_limits and len(u0), and derives dx. It cannot check whether the caller actually sampled the intended profile at those locations.

For this exercise, assume ordinary real numeric lists, NumPy arrays, and scalar parameters, with magnitudes for which grid and time-step calculations are representable in floating point. We do not require general-purpose type validation or special handling of extreme floating-point ranges. courant is the target magnitude $q$, not the signed realized $C$.

Return a dictionary with these keys:

The step budget is a simple limit on requested work, not a runtime or memory guarantee. Nonuniform grids, variable speeds, and changing inflow values are outside this solver’s scope. Keep the implementation focused on this brief.

4. Put the brief into a specification¶

Copy the template below into your notebook and fill it from the brief. You are organizing supplied requirements, not inventing a new solver design. For each validation requirement, name the defect or assumption it addresses and one way to check the promised behavior.

Outcome: What calculation should this solver deliver?
Context and assumptions: What equation, grid, and boundary model does it support?
Interface and outputs: Copy the signature and define every returned field.
Required behavior: State the update, direction selection, time selection,
    copying rules, and zero-evolution cases.
Rejected requests: List invalid inputs and the required exception behavior.
Acceptance evidence: Map each requirement to a hand result, invariant,
    comparison, or rejected-call test.
Agent permissions: State what the agent may add and what it must preserve.

Use the independent expectations in subsection 6 to complete the evidence section. If a requirement is unclear to you, resolve it against the brief before asking the agent to implement it.

5. Ask an agent to implement and test it¶

Attach your working notebook, including the original solver, derivations, and completed specification. Use this prompt:

Implement solve_transport according to the specification in my notebook.
Use the displayed FTBS update for positive speed, the derived forward-space
update for negative speed, and the derived final-time step selection.

Add only new, unexecuted cells at the end: one implementation cell, one
cell containing tests, and one short markdown cell mapping the tests to
requirements and identifying any limitations. Use only NumPy and the Python
standard library. Keep the numerical updates visible in the function and
use triple-single-quoted docstrings. Do not add a package or testing framework.

Use straightforward checks for the stated contract. Do not broaden the
accepted input types or add general-purpose validation machinery. Assume
ordinary real numeric inputs in representable ranges, as stated in the brief.
Do not add special handling for mixed Boolean data, extreme floating-point
values, or roundoff-triggered changes to the step count. Identify concerns
outside the contract in a short limitations note rather than implementing them.

Preserve every existing cell. Do not execute code, install anything, access
external resources, or write my verdict. If the specification has a
consequential ambiguity, ask rather than invent behavior.

Write callable tests using assertions and numpy.testing. Give them a
run_solver_checks(solver) entry point so I can test a deliberately defective
replacement later. Use the independent expected values in my specification;
do not compute expected answers by duplicating the solver's algorithm.
Use the supplied acceptance cases to test the checks, reported quantities,
and numerical results, including both speed directions. Keep the suite small. For refused calls, check the exception type and that
the message explains the failed requirement, without requiring exact wording.

The agent supplies implementation and test-writing labor. You remain responsible for checking both: a passing test written from the same mistaken assumption as the solver is weak evidence.

6. Audit, run, and defend the result¶

Before executing, inspect the function against your specification. Check that conversion to floating point precedes evolution, the two slices use old values, negative speed changes both stencil and inflow endpoint, and the time-selection and budget checks occur before the loop. Inspect the tests too: do they use independent answers and actually exercise the intended branch?

Use these acceptance cases. Unless stated otherwise, use floating-point inputs and the default step budget.

Use np.testing.assert_allclose for computed floating-point values; small roundoff differences are not implementation defects. Use exact checks for step counts, required keys, and exception types. In the final-time case, allow a small floating-point tolerance when checking that the realized Courant magnitude does not exceed the target; this is a comparison tolerance, not permission to choose a materially larger step.

After reviewing the cells, run the tests yourself. Record each failure as an implementation defect, a faulty test, or a specification ambiguity, with a reason. Send the agent a focused revision request and rerun affected checks after inspecting the change.

Finally, check that the tests can detect a known mistake. Add this temporary replacement in your working notebook; it deliberately corrupts the right inflow for negative speed when endpoint values differ:

def solve_transport_wrong_inflow(u0, x_limits, c, t_final,
                                 courant=0.8, max_steps=100_000):
    '''Deliberately impose the wrong inflow value for negative speed.'''
    result = solve_transport(u0, x_limits, c, t_final,
                             courant=courant, max_steps=max_steps)
    if c < 0 and result['num_steps'] > 0:
        result['u'][-1] = float(u0[0])
    return result

run_solver_checks(solve_transport_wrong_inflow)

An assertion should fail on a negative-speed case with unequal endpoint values. Identify that test and explain the mismatch. If the suite passes, it has missed the injected defect: improve the test, not the defective replacement. Keep this deliberate failure separate from the cell that runs the checks on the accepted solver.

Debrief: what became more robust?¶

The caller no longer has to coordinate grid spacing, time step, spatial bias, inflow endpoint, and step count by hand. The solver handles those dependent choices, reports what it did, and refuses the invalid requests named in its contract. The tests provide evidence that the promised behavior is implemented, and the injected defect checks that an important test can fail for the right reason.

These are real responsibilities in scientific software. We collected them in one small function so that you could inspect them. In a larger application, input processing, mesh and field objects, solver routines, and a separate test suite can share these responsibilities. A numerical update need not repeat every input check at every step. Where to put a check, and what it costs, are design decisions.

This solver has a deliberately limited contract. We have not investigated every input type, overflow or underflow, or extreme grid and time scales. Nor do passing tests and a stable update establish convergence, sufficient accuracy on the chosen grid, or a valid physical model. Further use requires evidence appropriate to that use; high-consequence work would demand much more than this notebook exercise.

The agent reduces the labor of drafting implementation and tests. It does not remove the work of understanding the contract, inspecting the code, and obtaining independent evidence. Judge the improvement by the failures prevented and the behavior verified, rather than by the number of checks the agent adds.