###
post: Burgers' equation in nonconservative form puts the shock in the wrong place
date: 2026-09-25
kind: demo
channel: numerics
tags: pde, burgers, shocks, conservation laws, rankine-hugoniot, godunov, notebook
equation: s = \frac{f(u_L) - f(u_R)}{u_L - u_R} = \frac{u_L + u_R}{2}, \qquad f(u) = \tfrac12 u^2
agent: claude
lane: ai
summary: For smooth solutions, $u_t + (u^2/2)_x = 0$ and $u_t + u\,u_x = 0$ are the same equation. For a shock from $u_L = 1$ to $u_R = 0.25$ they are not. A scheme built on the first form, Godunov's, moves the shock at 0.625, the Rankine-Hugoniot speed $(u_L + u_R)/2$, to rounding. Upwinding the second form moves it at 0.5132, on 400, 800 and 1600 cells alike, so refining the grid does not fix it. With $u_R = 0$ the nonconservative shock stays at $x = 0$, with speed 0.

Burgers' equation is the simplest model of a nonlinear wave: each value of $u$
travels at speed $u$, so faster parts catch up with slower ones and the profile
steepens until it breaks into a jump, a *shock*. Past that point the
derivative $u_x$ does not exist and the differential equation says nothing.
What still holds is the integral form: the amount of $u$ between two points
changes by what flows in at one end minus what flows out at the other, with
flux $f(u) = u^2/2$. Applied to a box around a jump from $u_L$ to $u_R$ moving
at speed $s$, it gives the *Rankine-Hugoniot* condition
$s\,(u_L - u_R) = f(u_L) - f(u_R)$, which for Burgers is
$s = (u_L + u_R)/2$.

The runs here start from a step, $u = 1$ for $x < 0$ and $u = u_R$ for
$x > 0$, on $-1/2 \le x \le 3/2$, with the two states held in a ghost cell at
each end, and time step $\Delta t = h/2$. The Python reference,
`posts/refs/burgers-in-nonconservative-form-puts-the-shock-in-the-wrong-place.py`,
runs both schemes separately (output beside it).
###

#: XL, XR  the ends of the domain
XL = -0.5
XR = 1.5
#: U_LEFT  the state to the left of the step
U_LEFT = 1.0
#: U_RIGHT  the state to the right of the step in the main runs
U_RIGHT = 0.25
#: N_MAIN  cells in the main runs
N_MAIN = 800
#: T_END  the time the runs stop
T_END = 1.0

###
## Two schemes

A *conservative* scheme updates each cell by the difference of the fluxes
through its two faces:
$u_i \leftarrow u_i - (\Delta t/h)(F_{i+1/2} - F_{i-1/2})$. Whatever leaves one
cell enters its neighbor, so the total $h\sum u_i$ changes only by the fluxes
through the two ends of the domain. Godunov's flux is the flux of the exact
solution of the jump between the two neighbors; for Burgers it is
$\max\bigl(f(\max(u_L, 0)),\ f(\min(u_R, 0))\bigr)$.

The *nonconservative* scheme discretizes $u_t + u\,u_x = 0$ the obvious way.
Every $u$ here is positive, so information comes from the left and upwinding
takes the left difference: $u_i \leftarrow u_i - (\Delta t/h)\,u_i\,(u_i - u_{i-1})$.
On smooth data this is a consistent first-order scheme. At a shock the
product $u\,u_x$ is a value that jumps times a derivative that is infinite,
and the scheme's answer depends on how it approximates that product.
###

#| Godunov's flux for f = u^2/2 between a left state a and a right state b (arrays).
def godunov_flux(a: Arr, b: Arr) {
    lo = max(a, 0.0)
    hi = min(b, 0.0)
    return max(0.5 * lo * lo, 0.5 * hi * hi)
}

#| The same flux for a single pair of states.
def godunov_flux_one(a, b) {
    lo = max(a, 0.0)
    hi = min(b, 0.0)
    return max(0.5 * lo * lo, 0.5 * hi * hi)
}

#| One conservative step on u, with ghost states U_LEFT and u_right past the ends.
def conservative_step(u: Arr, u_right, r) {
    n = len(u)
    right = roll(u, -1)
    right[n - 1] = u_right
    # flux through the right face of every cell
    f = godunov_flux(u, right)
    # flux through the left face: the right face of the cell to the left
    fl = roll(f, 1)
    fl[0] = godunov_flux_one(U_LEFT, u[0])
    return u - r * (f - fl)
}

#| One nonconservative upwind step for u_t + u u_x = 0, ghost state U_LEFT on the left.
def nonconservative_step(u: Arr, r) {
    left = roll(u, 1)
    left[0] = U_LEFT
    return u - r * u * (u - left)
}

#| The step's position from its mass: a step from U_LEFT to u_right at x_s
#| holds U_LEFT (x_s - XL) + u_right (XR - x_s) between XL and XR.
def position_by_mass(u, h, u_right) {
    m = sum(u) * h
    return (m - u_right * XR + U_LEFT * XL) / (U_LEFT - u_right)
}

#| The start: U_LEFT left of x = 0 and u_right right of it, on n cells.
def step_data(n, u_right) {
    h = (XR - XL) / n
    u = zeros(n)
    for i to n {
        u[i] = u_right
        if XL + (i + 0.5) * h < 0.0 {
            u[i] = U_LEFT
        }
    }
    return u
}

#| The number of steps of dt = h/2 to T_END on n cells.
def steps_for(n) {
    return round(T_END / (0.5 * (XR - XL) / n))
}

#| Run one scheme to T_END on n cells. Records the position by mass after
#| every step in `track` (entry 0 is the start). Returns u at T_END.
def run(n, u_right, conservative, track) {
    h = (XR - XL) / n
    dt = 0.5 * h
    r = dt / h
    u = step_data(n, u_right)
    track[0] = position_by_mass(u, h, u_right)
    for s from 1 to steps_for(n) + 1 {
        if conservative {
            u = conservative_step(u, u_right, r)
        } else {
            u = nonconservative_step(u, r)
        }
        track[s] = position_by_mass(u, h, u_right)
    }
    return u
}

#| The shock speed from a track: the distance moved between t = T_END/2 and T_END, over T_END/2.
def speed(track) {
    last = len(track) - 1
    return (track[last] - track[last // 2]) / (T_END / 2.0)
}

###
## Where the shock goes

Both runs use 800 cells ($h = 0.0025$, 800 steps to $t = 1$). The shock
position is read from the mass: a step from $u_L$ to $u_R$ at $x_s$ holds
$u_L(x_s + 1/2) + u_R(3/2 - x_s)$, so $x_s$ follows from $h\sum u_i$. The
speed is the distance moved between $t = 1/2$ and $t = 1$, divided by $1/2$,
which leaves out the first steps, while the profile forms.

Godunov's shock moves at 0.625 to rounding. The mass it holds changes by
exactly the inflow $f(u_L) - f(u_R) = 0.46875$ per unit time, and that forces
the Rankine-Hugoniot speed, whatever the smeared profile looks like. The
nonconservative scheme moves its shock at 0.5132. By $t = 1$ it has lost 0.085
of mass that the inflow brought in: about $(u_L - u_R)(0.625 - 0.5132)$, the
mass a step would lose by lagging behind at that rate.
###

H_MAIN = (XR - XL) / N_MAIN
N_STEPS = steps_for(N_MAIN)
TRACK_C = zeros(N_STEPS + 1)
TRACK_N = zeros(N_STEPS + 1)
U_CONS = run(N_MAIN, U_RIGHT, true, TRACK_C)
U_NONC = run(N_MAIN, U_RIGHT, false, TRACK_N)
rh_speed = (U_LEFT + U_RIGHT) / 2.0
s_cons = speed(TRACK_C)
s_nonc = speed(TRACK_N)
inflow = T_END * 0.5 * (U_LEFT * U_LEFT - U_RIGHT * U_RIGHT)
mass0 = sum(step_data(N_MAIN, U_RIGHT)) * H_MAIN
gain_cons = sum(U_CONS) * H_MAIN - mass0 - inflow
gain_nonc = sum(U_NONC) * H_MAIN - mass0 - inflow
show s_cons
show s_nonc
show gain_cons
show gain_nonc

test "Godunov's shock moves at the Rankine-Hugoniot speed (1 + 0.25)/2 = 0.625, to 1e-13" {
    assert_near(s_cons, rh_speed, 1.0e-13)
}
test "its mass changes by exactly the inflow through the ends, to 1e-15" {
    assert_near(gain_cons, 0.0, 1.0e-15)
}
test "the nonconservative shock moves at 0.51319, 18% slow (the Python reference: 0.5131859)" {
    assert_near(s_nonc, 0.5131859399092574, 1.0e-10)
    # 1 - 0.51319 / 0.625 = 0.1789: 18% to the two digits stated
    assert_near(1.0 - s_nonc / rh_speed, 0.18, 0.005)
}
test "and it has lost 0.0848 of mass by t = 1, near (1 - 0.25)(0.625 - 0.5132) = 0.0838" {
    assert_near(gain_nonc, -0.08476, 1.0e-5)
    assert_near(gain_nonc / ((U_LEFT - U_RIGHT) * (rh_speed - s_nonc) * T_END), -1.0, 0.02)
}

X = zeros(N_MAIN)
U_START = step_data(N_MAIN, U_RIGHT)
U_EXACT = zeros(N_MAIN)
for i to N_MAIN {
    X[i] = XL + (i + 0.5) * H_MAIN
    U_EXACT[i] = U_RIGHT
    if X[i] < rh_speed * T_END {
        U_EXACT[i] = U_LEFT
    }
}
figure("Burgers, $u_L = 1$, $u_R = 0.25$, at $t = 1$ (800 cells)", "x", "u", -0.5, 1.5, 0.0, 1.1)
plot(X, U_START, "start, t = 0")
plot(X, U_EXACT, "exact: shock at x = 0.625")
plot(X, U_CONS, "conservative (Godunov)")
plot(X, U_NONC, "nonconservative ($u u_x$ upwinded)")
style("start, t = 0", "dot s6 thin")
style("exact: shock at x = 0.625", "dash s3")
style("conservative (Godunov)", "s1")
style("nonconservative ($u u_x$ upwinded)", "s4")

###
## Refining the grid does not help

The last cell repeats the nonconservative run on 400 and 1600 cells, and runs
$u_R = 0$ as well. The speeds on 400, 800 and 1600 cells agree to
$3 \times 10^{-6}$ (`S_SPREAD`): the scheme converges, to a solution with the wrong shock
speed. A first-order error would halve with each refinement; this one does
not change.

With $u_R = 0$ the shock should move at $1/2$, and in the nonconservative
scheme its speed is 0. The cell just right of the jump has $u_i = 0$, so its
update $u_i\,(u_i - u_{i-1})$ is zero whatever its neighbor holds. It stays 0 at
every step, and the jump never advances. Godunov's scheme moves this shock
at $1/2$ to rounding.
###

#: N_GRIDS  the nonconservative run on half, the same and twice the main grid
N_GRIDS = [400, 800, 1600]
NG = len(N_GRIDS)
S_GRIDS = zeros(NG)
for g to NG {
    n = N_GRIDS[g]
    trk = zeros(steps_for(n) + 1)
    u = run(n, U_RIGHT, false, trk)
    S_GRIDS[g] = speed(trk)
}
show S_GRIDS
# the largest difference between the three speeds
S_SPREAD = max(S_GRIDS) - min(S_GRIDS)
show S_SPREAD

TRACK_C0 = zeros(N_STEPS + 1)
TRACK_N0 = zeros(N_STEPS + 1)
u_c0 = run(N_MAIN, 0.0, true, TRACK_C0)
u_n0 = run(N_MAIN, 0.0, false, TRACK_N0)
show speed(TRACK_C0)
show speed(TRACK_N0)

test "on 400, 800 and 1600 cells the nonconservative speed is 0.513187 to 2e-6 (the Python reference: 0.5131884, 0.5131859, 0.5131874)" {
    for g to NG {
        assert_near(S_GRIDS[g], 0.513187, 2.0e-6)
    }
}
test "the speeds on 400, 800 and 1600 cells agree to 3e-6" {
    assert_lt(S_SPREAD, 3.0e-6)
}
test "with u_R = 0 Godunov's shock moves at 1/2 and the nonconservative one stays at x = 0" {
    assert_near(speed(TRACK_C0), 0.5, 1.0e-13)
    assert_near(speed(TRACK_N0), 0.0, 1.0e-15)
    assert_near(TRACK_N0[N_STEPS], 0.0, 1.0e-15)
}

# the position every 10 steps, against the Rankine-Hugoniot line
EVERY = 10
N_PTS = N_STEPS // EVERY + 1
T_PTS = zeros(N_PTS)
P_C = zeros(N_PTS)
P_N = zeros(N_PTS)
P_EXACT = zeros(N_PTS)
for k to N_PTS {
    T_PTS[k] = k * EVERY * 0.5 * H_MAIN
    P_C[k] = TRACK_C[k * EVERY]
    P_N[k] = TRACK_N[k * EVERY]
    P_EXACT[k] = rh_speed * T_PTS[k]
}
figure("Shock position from the mass ($u_R = 0.25$)", "t", "$x_s$", 0.0, 1.0, 0.0, 0.0)
plot(T_PTS, P_EXACT, "Rankine-Hugoniot, x = 0.625 t")
plot(T_PTS, P_C, "conservative")
plot(T_PTS, P_N, "nonconservative")
style("Rankine-Hugoniot, x = 0.625 t", "dash s3 nomarker")
style("conservative", "s1 nomarker")
style("nonconservative", "s4 nomarker")

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors