Space & Time

In the first module, we studied numerical integration methods for the solution of ordinary differential equations (ODEs), using the phugoid model of glider flight as a motivation. In this module, we will study the numerical solution of partial differential equations (PDEs), where the unknown is a multi-variate function. The problem could depend on time, $t$, and one spatial dimension $x$ (or more), which means we need to build a discretization grid with each independent variable.

We will start our discussion of numerical PDEs with 1D linear and non-linear convection equations, the 1D diffusion equation, and 1D Burgers’ equation.

1D linear convection¶

The one-dimensional linear convection equation is the simplest, most basic model that can be used to learn something about numerical solution of PDEs. It’s surprising that this little equation can teach us so much! Here it is:

The equation represents a wave propagating with speed $c$ in the $x$ direction, without change of shape. For that reason, it’s sometimes called the one-way wave equation (sometimes also the advection equation).

With an initial condition $u(x,0)=u_0(x)$, the equation has an exact solution given by:

Go on: check it. Take the time and space derivative and stick them into the equation to see that it holds.

Look at the exact solution for a moment ... we know two things about it:

1. its shape does not change, being always the same as the initial wave, $u_0$, only shifted in the $x$-direction; and

2. it’s constant along so-called characteristic curves, $x-ct=$constant. This means that for any point in space and time, you can move back along the characteristic curve to $t=0$ to know the value of the solution.

Why do we call the equations linear? PDEs can be either linear or non-linear. In a linear equation, the unknown function $u$ and its derivatives appear only in linear terms, in other words, there are no products, powers, or transcendental functions applied on them.

The most important feature of linear equations is: solutions can be superposed to generate new solutions that still satisfy the original equation. This is super useful!

Finite-differences¶

In the previous lessons, we discretized time derivatives; now we have derivatives in both space and time, so we need to discretize with respect to both these variables.

Imagine a space-time plot, where the coordinates in the vertical direction represent advancing in time—for example, from $t^n$ to $t^{n+1}$—and the coordinates in the horizontal direction move in space: consecutive points are $x_{i-1}$, $x_i$, and $x_{i+1}$. This creates a grid where a point has both a temporal and spatial index. Here is a graphical representation of the space-time grid:

For the numerical solution of $u(x,t)$, we’ll use subscripts to denote the spatial position, like $u_i$, and superscripts to denote the temporal instant, like $u^n$. We would then label the solution at the top-middle point in the grid above as follows: $u^{n+1}_{i}$.

Each grid point below has an subscript index, corresponding to the spatial position and increasing to the right, and an superscript index, corresponding to the time instant and increasing upwards. A small grid segment would have the following values of the numerical solution at each point:

Another way to explain our discretization grid is to say that it is built with constant steps in time and space, $\Delta t$ and $\Delta x$, as follows:

Discretizing our model equation¶

Let’s see how to discretize the 1D linear convection equation in both space and time. By definition, the partial derivative with respect to time differentiates only with time and not with space; its discretized form changes only the $n$ indices. Similarly, the partial derivative with respect to $x$ differentiates with space not time, and only the $i$ indices are affected.

We’ll discretize the spatial coordinate $x$ into points indexed from $i=0$ to $N$, and then step in discrete time intervals of size $\Delta t$.

From the definition of a derivative (and simply removing the limit), we know that for $\Delta x$ sufficiently small:

This formula could be applied at any point $x_i$. But note that it’s not the only way that we can estimate the derivative. The geometrical interpretation of the first derivative $\partial u/ \partial x$ at any point is that it represents the slope of the tangent to the curve $u(x)$. In the sketch below, we show a slope line at $x_i$ and mark it as “exact.” If the formula written above is applied at $x_i$, it approximates the derivative using the next spatial grid point: it is then called a forward difference formula.

But as shown in the sketch below, we could also estimate the spatial derivative using the point behind $x_i$, in which case it is called a backward difference. We could even use the two points on each side of $x_i$, and obtain what’s called a central difference (but in that case the denominator would be $2\Delta x$).

We have three possible ways to represent a discrete form of $\partial u/ \partial x$:

• Forward difference: uses $x_i$ and $x_i + \Delta x$,

• Backward difference: uses $x_i$ and $x_i- \Delta x$,

• Central difference: uses two points on either side of $x_i$.

The sketch above also suggests that some finite-difference formulas might be better than others: it looks like the central difference approximation is closer to the slope of the “exact” derivative. We’ll see later how to make this observation rigorous.

The three formulas are:

Euler’s method is equivalent to using a forward-difference scheme for the time derivative. Let’s stick with that, and choose the backward-difference scheme for the space derivative. Our discrete equation is then:

where $n$ and $n+1$ are two consecutive steps in time, while $i-1$ and $i$ are two neighboring points of the discretized $x$ coordinate. With given initial conditions, the only unknown in this discretization is $u_i^{n+1}$. We solve for this unknown to get an equation that lets us step in time, as follows:

It is useful to sketch a grid segment, showing the grid points that influence our numerical solution. This is called a stencil. Below is the stencil for solving our model equation with the finite-difference formula we wrote above.

And compute!¶

Alright. Let’s get a little Python on the road. First: we need to load our array and plotting libraries, as usual.

We also set notebook-wide plotting parameters for the font family and the font size by modifying entries of the rcParams dictionary.

As a first exercise, we’ll solve the 1D linear convection equation with a square wave initial condition, defined as follows:

We also need a boundary condition on $x$: let $u=1$ at $x=0$. Our spatial domain for the numerical solution will only cover the range $x\in (0, 2)$.

Now let’s define a few variables; we want to make an evenly spaced grid of points within our spatial domain. In the code below, we define a variable called nx that will be the number of spatial grid points, and a variable dx that will be the distance between any pair of adjacent grid points. We also can define a step in time, dt, a number of steps, nt, and a value for the wave speed: we like to keep things simple and make $c=1$.

We also need to set up our initial conditions. Here, we use the NumPy function np.ones() defining an array which is nx-element long with every value equal to 1. How useful! We then change a slice of that array to the value $u=2$, to get the square wave, and we print out the initial array just to admire it. But which values should we change? The problem states that we need to change the indices of u such that the square wave begins at $x = 0.5$ and ends at $x = 1$.

We can use the np.where() function to return a list of indices where the vector $x$ meets some conditions. The function np.logical_and() computes the truth value of x >= 0.5 and x <= 1.0, element-wise.

With the list of indices, we can now update our initial conditions to get a square-wave shape.

Now let’s take a look at those initial conditions we’ve built with a handy plot.

It does look pretty close to what we expected. But it looks like the sides of the square wave are not perfectly vertical. Is that right? Think for a bit.

Now it’s time to write some code for the discrete form of the convection equation using our chosen finite-difference scheme.

For every element of our array u, we need to perform the operation:

We’ll store the result in a new (temporary) array un, which will be the solution $u$ for the next time-step. We will repeat this operation for as many time-steps as we specify and then we can see how far the wave has traveled.

We first initialize the placeholder array u to hold the values we calculate for the $n+1$ time step, beginning with a copy of the initial condition.

Then, we may think we have two iterative operations: one in space and one in time (we’ll learn differently later), so we may start by nesting a spatial loop inside the time loop, as shown below. You see that the code for the finite-difference scheme is a direct expression of the discrete equation.

Here, nt counts updates rather than stored time levels. Because range(nt) starts at zero and performs exactly nt passes, the final physical time in this example is nt * dt.

Note 1—We stressed above that our physical problem needs a boundary condition at $x=0$. Here we do not need to impose it at every iteration because our discretization does not change the value of u[0]: it remains equal to one and our boundary condition is therefore satisfied during the whole computation!

Note 2—We will learn later that the code as written above is quite inefficient, and there are better ways to write this, Python-style. But let’s carry on.

Now compare the computed profile with the exact square wave translated to the same physical time.

That’s funny. Our square wave has definitely moved to the right, but it’s no longer in the shape of a top-hat. What’s going on?

Does a finer grid always help?¶

The rounded profile suggests that a finer spatial grid might improve the solution. Test that idea while keeping dt, nt, and $c$ fixed: only the number of spatial points changes below. Before executing the cell, predict what each refinement will do. Then compare the computed profiles with the exact translated square wave and record whether the numerical values remain between 1 and 2.

Treat this as an observation rather than a convergence study, and do not change dt to repair any surprising result. We will return to what happens—and why—in the next lesson.

Spatial truncation error¶

Recall the backward-difference approximation we are using for the spatial derivative:

We obtain it by using the definition of the derivative at a point, and simply removing the limit, in the assumption that $\Delta x$ is very small. But we already learned with Euler’s method that this introduces an error, called the truncation error.

We can determine the spatial order of this error by expanding $u(x_i-\Delta x)$ in a Taylor series about $x_i$. Complete that step on paper before reading the result.

The dominant term that is neglected in the finite-difference approximation is of $\mathcal{O}(\Delta x)$. We also see that the approximation converges to the exact derivative as $\Delta x \rightarrow 0$. That’s good news!

In summary, the chosen “forward-time/backward space” difference scheme is first-order in both space and time: the truncation errors are $\mathcal{O}(\Delta t, \Delta x)$. We’ll come back to this!

A controlled spatial refinement study¶

The square wave makes numerical diffusion easy to see, but its jumps do not satisfy the smoothness assumed in the Taylor expansion above. To measure the formal spatial behavior of the scheme, use a smooth Gaussian pulse on the same unit background. Its exact solution is the same profile translated by $c t$.

For each grid, compare the numerical and exact solutions at the same final time with the discrete $L_1$ error in Equation 14:

When the grid spacing is reduced, estimate the observed order from two consecutive errors with Equation 15:

Use a very small time step for every spatial grid so that time-discretization error is much smaller than the spatial error over this range. Repeating the study with half that time step provides a sensitivity check: the spatial conclusions should barely change.

The error decreases by nearly a factor of two whenever $\Delta x$ is halved, and the observed order approaches 1. Halving the already small time step barely changes either the errors or the orders. Together, those observations support the first-order spatial behavior predicted by the truncation-error analysis over this tested grid range.

This controlled result does not explain the failed refinement in the earlier square-wave experiment. Keep that discrepancy in mind for the next lesson.

Non-linear convection¶

Let’s move on to the non-linear convection equation, using the same methods as before. The 1D convection equation is:

The only difference with the linear case is that we’ve replaced the constant wave speed $c$ by the variable speed $u$. The equation is non-linear because now we have a product of the solution and one of its derivatives: the product $u\,\partial u/\partial x$. This changes everything!

We’re going to use the same discretization as for linear convection: forward difference in time and backward difference in space. Here is the discretized equation:

Solving for the only unknown term, $u_i^{n+1}$, gives an equation that can be used to advance in time:

There is very little that needs to change from the code written so far. In fact, we’ll even use the same square-wave initial condition. But let’s re-initialize the variable u with the initial values, and re-enter the numerical parameters here, for convenience (we no longer need $c$, though).

How does it look?

Changing just one line of code in the solution of linear convection, we are able to now get the non-linear solution: the line that corresponds to the discrete equation now has un[i] in the place where before we just had c. So you could write something like:

for n in range(nt):  
  un = u.copy() 
  for i in range(1, nx): 
    u[i] = un[i] - un[i]*dt/dx*(un[i]-un[i-1]) 

We’re going to be more clever than that and use NumPy to update all values of the spatial grid in one fell swoop. We don’t really need to write a line of code that gets executed for each value of $u$ on the spatial grid. Python can update them all at once! Study the code below, and compare it with the one above. Here is a helpful sketch, to illustrate the array operation—also called a “vectorized” operation—for $u_i-u_{i-1}$.

Hmm. That’s quite interesting: like in the linear case, we see that we have lost the sharp sides of our initial square wave, but there’s more. Now, the wave has also lost symmetry! It seems to be lagging on the rear side, while the front of the wave is steepening. Is this another form of numerical error, do you ask? No! It’s physics!

Same equation, different algorithms¶

The front of the wave steepens because larger values of $u$ move faster than smaller values. As that steepening produces a sharp jump, a choice that looked like ordinary algebra begins to matter numerically.

For a differentiable function, the chain rule gives

We can therefore write the nonlinear convection equation in conservative form:

As long as $u$ is smooth, Equation 16 and Equation 20 describe the same continuous equation. They do not, however, lead to the same discrete algorithm.

Discretize the two forms¶

The update used in the preceding computation applies a backward difference directly to $u\,\partial u/\partial x$; we will call it the pointwise update. It is Equation 18. Applying the same forward-time, backward-space pattern to the flux derivative instead gives the conservative update:

The difference becomes visible by factoring the flux difference:

The pointwise update multiplies the backward difference by $u_i$, while the conservative update effectively uses the average $(u_i+u_{i-1})/2$. When neighboring values are close, these multipliers become close under refinement. Across a jump, their difference is not small.

A balance built into the algorithm¶

The word conservative describes an algebraic property of the discrete update. Sum Equation 21 over the updated nodes $i=1,\ldots,N$:

Every interior flux appears once with a plus sign and once with a minus sign, so the sum telescopes: only the two boundary fluxes remain.

If we move the boundary-flux term in Equation 23 to the left and substitute a computed solution, the number left over is the discrete balance residual. More generally, a residual measures the amount by which computed values fail to satisfy an equation. For the conservative update, this residual is zero in exact arithmetic and should be close to zero in floating-point arithmetic. Over several time steps, we add the boundary contribution from each step and compare it with the total change from the initial state to the final state; the agent activity below calls this quantity the cumulative balance residual, $R$. A small residual verifies this particular balance identity; it does not, by itself, prove that the flux or the rest of the implementation is correct.

If the endpoint fluxes are equal, the discrete total $\Delta x\sum_i u_i$ is unchanged apart from round-off. The pointwise update does not have this flux-difference structure and does not satisfy the same identity. This cancellation principle is central to conservative methods for nonlinear hyperbolic problems LeVeque, 2002.

When the solution remains differentiable, the chain rule supports the equivalence of the two continuous forms, and consistent discretizations should approach the same smooth solution. At a jump, the ordinary derivative used in the pointwise argument is no longer defined there, so that chain-rule argument is not enough to determine how the jump propagates. The conservative form retains a meaningful balance across the jump.

This short bridge does not yet develop the integral conservation law, weak solutions, shock speeds, or general numerical fluxes. Those ideas belong in Module 3. For now, the important warning is that algebraically equivalent smooth equations can produce algorithms with materially different behavior once discontinuities enter the computation.

Check the functions against your hand calculation before asking an agent to build on them. The constant-state checks exercise a different property, so keep both kinds of evidence.

Compare the algorithms with an agent¶

The two update functions are short enough to derive, implement, and check yourself. Comparing them over several grids and tracking a flux balance is repetitive work that can usefully be delegated. The agent will first draft the comparison driver and later help you search for a blind spot in one diagnostic. In both roles, the agent proposes work for you to inspect; it will not choose the equations, change the algorithms, execute the code, or decide what the evidence means.

Record expectations before delegation¶

In your notebook, record your answers to these questions before attaching it to an agent:

1. What maximum change do you expect after either function advances a constant state by one step?

2. What magnitude do you expect for the conservative scheme’s balance residual, obtained by moving all terms in Equation 23 to one side?

3. Does the pointwise update have an algebraic reason to produce the same residual?

4. As the grid is refined, should the two algorithms approach one another more clearly for the smooth hump or for the discontinuous square pulse? Why?

These are hypotheses to test, not outputs that the agent should be instructed to manufacture.

Use this bounded task brief¶

Add the following brief as a Markdown cell in your notebook. Read every requirement and resolve any mismatch with your own function names before delegating.

Agent algorithm-comparison task brief¶

• Goal and scope: Add a comparison driver that compares my existing pointwise_step() and conservative_step() functions. Do not implement, rewrite, or repair either method. Prefer named functions and straightforward, readable loops. The comparison concerns discrete behavior only; deriving shock speeds, introducing another scheme, and making the final numerical judgment are out of scope.

• Problem data: Use $x\in[0,2]$, nx_values = [81, 161, 321, 641, 1281], and t_compare = 0.15. For each grid, start with a target dt = 0.2 * dx, round the step count to reach the common final time, and then reset dt = t_compare / num_steps so every run ends at exactly the same time. Compare (a) the smooth hump $1+\exp[-((x-0.7)/0.2)^2]$ and (b) the square pulse with background 1 and value 2 on $0.5\le x\le1$. Each run must begin from a fresh copy of the same initial array.

• Required diagnostics: First report the maximum one-step change produced by each method from a constant array with value 1.5. Write a reusable runner named run_balance_comparison(step_function, flux_function, initial_state, dt, dx, num_steps). It must return a dictionary containing u_final, initial_total, final_total, total_change, and balance_residual.

  For every run, calculate these three diagnostics:

  $$\begin{aligned} M(u) &= \Delta x\sum_i u_i,\\[0.75em] R &= M_{\text{final}}-M_{\text{initial}} +\sum_n\Delta t\left[f(u_N^n)-f(u_0^n)\right],\\[0.75em] D &= \Delta x\sum_i \left|u_i^{\text{pointwise}}-u_i^{\text{conservative}}\right|. \end{aligned}$$

  (24)

  Here, $M$ is the discrete total, $R$ is the cumulative balance residual, and $D$ is the inter-algorithm difference at the final time. Evaluate $M$ for the initial and final states. In $R$, use the boundary values from the state before each update. Evaluate $D$ using the two final solutions on the same grid.

• Required output: Produce one readable table for each profile with nx, dx, adjusted dt, step count, $D$, both balance residuals, and both initial-to-final total changes. Produce a figure showing both final solutions on the finest grid for each profile and a log-log plot of $D$ against $\Delta x$. Label methods, profiles, axes, and final time. Expose numerical values; do not replace them with pass/fail messages or a conclusion. Use only the NumPy and Matplotlib imports already present.

• Permissions and stopping rule: Inspect this notebook and add exactly three Python cells followed by one Markdown cell that records the added functions, experiment, plots, and any assumptions. Do not execute any cell, alter existing cells, write another file, install anything, or use the network. Stop after inserting those four cells and report that the comparison driver is ready for human audit.

• Done when: The unexecuted cells contain the constant-state comparison, reusable balance runner, two-profile grid study, required tables and plots, plus the change record. An expected trend is not a completion condition; preserve and expose surprising results.

The brief separates the model and evidence you fixed from the code organization the agent may choose. Compare it with the minimum sufficient specification before granting access.

Invoke the agent¶

Save your notebook so the attachment contains your checked step functions, independent expectations, and task brief. Attach that notebook to your agent interface and send this short request:

If you do not have access to an agent that can edit a notebook, ask it to return the same four proposed cells in chat, add them yourself only after inspection, and follow the audit below.

Audit before execution¶

Inspect the notebook diff or the four inserted cells before running them. Confirm that the code:

• calls your existing step functions without redefining or modifying them;

• starts every method and grid from an independent copy of identical initial data;

• uses the requested grids, profiles, and common physical time, and reports the adjusted dt and step count;

• accumulates the boundary-flux term from the old state before each call to the step function;

• uses the plus sign in the balance residual derived from Equation 23, with right-boundary flux minus left-boundary flux;

• computes $D$ with the two solutions on the same grid and includes the factor $\Delta x$;

• exposes totals and residuals for both methods rather than checking conservation only for the method named conservative_step; and

• changes no earlier cell or out-of-scope artifact.

If any item fails inspection, ask the agent to revise only its added cells, then inspect them again. Once the code passes, execute one inserted cell at a time. Recalculate one row’s $D$ and residual independently from the stored arrays. A clean plot or an agent statement that the method conserves is not a substitute for that arithmetic.

Test the diagnostics with known defects¶

A check that passes working code is not necessarily capable of detecting broken code. In this final exercise, you will challenge the diagnostics with two known defects and record which checks detect each one. Complete the steps in order, in new cells below the agent’s comparison driver.

First, save the accepted baseline results. Then add the two deliberately defective functions below under new names. Do not overwrite the correct functions. The first omits the factor $1/2$ in the intended flux; the second reads a newly updated left neighbor while marching across the array.

def conservative_step_missing_half(u, dt, dx):
    '''Deliberately use the wrong quadratic flux.'''
    u_next = u.copy()
    u_next[1:] = (
        u[1:] - dt / dx * (u[1:]**2 - u[:-1]**2)
    )
    return u_next


def pointwise_step_mixed_levels(u, dt, dx):
    '''Deliberately mix old and new time levels.'''
    u_next = u.copy()
    for i in range(1, u.size):
        u_next[i] = (
            u_next[i]
            - dt / dx * u_next[i]
            * (u_next[i] - u_next[i - 1])
        )
    return u_next

If you need to, review the Python refresher above that explained why the vectorized implementation of the spatial update is safe. The function pointwise_step_mixed_levels() is doing precisely what that admonition warned against: a “left-to-right” loop for the spatial update.

Step 1 — Check one update against known values¶

The earlier hand calculation gave the expected results for one update from $u=[1,2,1]$ with $\Delta t/\Delta x=0.1$. Reuse those values as an independent reference. Add and run this code in a new notebook cell:

u_defect_test = np.array([1.0, 2.0, 1.0])
pointwise_expected = np.array([1.0, 1.8, 1.1])
conservative_expected = np.array([1.0, 1.85, 1.15])

missing_half_actual = conservative_step_missing_half(
    u_defect_test, dt=0.1, dx=1.0
)
mixed_levels_actual = pointwise_step_mixed_levels(
    u_defect_test, dt=0.1, dx=1.0
)

print('Missing-half result: ', missing_half_actual)
print('Expected conservative:', conservative_expected)
print('Mixed-level result:  ', mixed_levels_actual)
print('Expected pointwise:  ', pointwise_expected)

missing_half_detected = not np.allclose(
    missing_half_actual, conservative_expected
)
mixed_levels_detected = not np.allclose(
    mixed_levels_actual, pointwise_expected
)
assert missing_half_detected
assert mixed_levels_detected

The cell should run without an assertion error. Here, a passing assert means that the comparison successfully detected a mismatch. You should obtain [1.0, 1.7, 1.3] from the missing-half function and [1.0, 1.8, 1.08] from the mixed-level function.

Record: the known-value check detects both defects. For the mixed-level function, notice that the first updated value is correct; the error appears at the next point, after the loop reads the newly updated left neighbor: i=2 uses the updated value of u[1]. Follow the arithmetic on paper using the update of $u=[1,2,1]$ with $\Delta t/\Delta x=0.1$.

Step 2 — Show what the constant-state check misses¶

Now apply both defective functions to a constant array. Add and run:

constant_defect_test = np.full(5, 1.5)

for method_name, step_function in (
    ('missing half', conservative_step_missing_half),
    ('mixed levels', pointwise_step_mixed_levels),
):
    actual = step_function(
        constant_defect_test, dt=0.1, dx=1.0
    )
    maximum_change = np.max(
        np.abs(actual - constant_defect_test)
    )
    print(f'{method_name:12s}: {maximum_change:.3e}')
    np.testing.assert_allclose(actual, constant_defect_test)

Both maximum changes should be zero, and both assertions should pass. This does not mean the functions are correct: every spatial difference is zero for a constant state, so neither defect is exercised.

Record: the constant-state check detects neither defect.

Step 3 — Use unequal boundary fluxes to detect the wrong flux¶

A conservative method may have a nonzero change in its discrete total when flux enters or leaves through the boundaries. The requirement is not total_change == 0; it is balance_residual == 0.

Use one step from an array whose endpoint values differ. The runner must receive the intended nonlinear_flux even when it calls the defective step function. Add and run:

u_balance_test = np.array([1.0, 2.0, 1.5])
balance_test_results = {}
for method_name, step_function in (
    ('conservative_step', conservative_step),
    (
        'conservative_step_missing_half',
        conservative_step_missing_half,
    ),
):
    balance_test_results[method_name] = run_balance_comparison(
        step_function, nonlinear_flux, u_balance_test.copy(),
        dt=0.1, dx=1.0, num_steps=1,
    )

print(f"{'method':>32} {'total change':>16} "
      f"{'balance residual':>20}")
for method_name, result in balance_test_results.items():
    print(
        f"{method_name:>32} "
        f"{result['total_change']:16.8e} "
        f"{result['balance_residual']:20.8e}"
    )

Before checking the output, calculate the intended boundary contribution:

The correct method should report total_change $=-0.0625$ and a residual near zero. The missing-half method should report total_change $=-0.125$ and balance_residual $=-0.0625$. Confirm those expectations explicitly:

correct_balance = balance_test_results['conservative_step']
wrong_balance = balance_test_results[
    'conservative_step_missing_half'
]

np.testing.assert_allclose(
    correct_balance['total_change'], -0.0625, atol=1e-14
)
np.testing.assert_allclose(
    correct_balance['balance_residual'], 0.0, atol=1e-14
)
np.testing.assert_allclose(
    wrong_balance['total_change'], -0.125, atol=1e-14
)
np.testing.assert_allclose(
    wrong_balance['balance_residual'], -0.0625, atol=1e-14
)

Record: with unequal endpoint fluxes, the balance residual detects the missing factor $1/2$. The correct total is not constant, but its change is fully explained by the boundary flux.

Step 4 — Expose the balance check’s blind spot¶

Test the red-team hypothesis with a square pulse on a background of 1. This profile is nonconstant but has equal endpoint values. Use the coarsest grid from the agent study so the check runs quickly:

nx_square = 81
x_square = np.linspace(0.0, 2.0, num=nx_square)
dx_square = 2.0 / (nx_square - 1)
num_steps_square = int(round(0.15 / (0.2 * dx_square)))
dt_square = 0.15 / num_steps_square
u_square = np.ones(nx_square)
u_square[np.logical_and(x_square >= 0.5, x_square <= 1.0)] = 2.0

correct_square = run_balance_comparison(
    conservative_step, nonlinear_flux, u_square.copy(),
    dt_square, dx_square, num_steps_square,
)
wrong_square = run_balance_comparison(
    conservative_step_missing_half, nonlinear_flux,
    u_square.copy(), dt_square, dx_square, num_steps_square,
)
square_difference = dx_square * np.sum(
    np.abs(correct_square['u_final'] - wrong_square['u_final'])
)

print('Correct residual:   ',
      correct_square['balance_residual'])
print('Wrong-flux residual:',
      wrong_square['balance_residual'])
print('Difference between final states:', square_difference)

np.testing.assert_allclose(
    correct_square['balance_residual'], 0.0, atol=1e-12
)
np.testing.assert_allclose(
    wrong_square['balance_residual'], 0.0, atol=1e-12
)
assert not np.allclose(
    correct_square['u_final'], wrong_square['u_final']
)

Both residuals should be near zero even though the final arrays disagree; with these settings, square_difference is about $4.48\times10^{-1}$. During this short run, the pulse remains separated from the boundaries and both endpoint values remain 1. The intended and incorrectly scaled fluxes therefore both have zero net boundary difference at every step. The global balance check cannot distinguish them in this case.

Record: the equal-endpoint square-pulse balance test misses the wrong flux. The nonzero solution difference shows disagreement, but without an independent reference it does not, by itself, identify which solution is correct.

Step 5 — Record what each diagnostic can detect¶

Add this completed detection table to your notebook beside the numerical output:

Your conclusion should say that a passing diagnostic supports only the property and test case it actually exercises. The known-value check exposes both injected errors; the constant-state check exposes neither; and flux balance detects the wrong flux only when the boundary data make the discrepancy visible. No single passing check proves that an implementation is correct.

See Test the checks for the general role of this adversarial step.

Completion rubric¶

This is a Complete/Revise checkpoint. The activity is complete only when the specification, code audit, executed evidence, defect challenge, and learner judgment are all visible in your notebook.

What’s next?¶

The next lesson explains the stability restriction behind the failed fixed-dt refinement experiment. Module 3 will return to the conservative form from a control-volume balance, define how conservation laws represent discontinuous solutions, determine shock propagation, and apply those ideas to traffic flow.