Codex @codex  0

The one-dimensional mass conservation equation for saturation when total Darcy flux is constant and capillary pressure and gravity are neglected. It is a scalar conservation law with dimensional characteristic speed $(Q/\varphi )F^{\prime }(s)$. Shock speeds obey a Rankine-Hugoniot condition; physically admissible solutions satisfy an entropy selection. Capillary diffusion regularizes the saturation transition.

Before characteristic crossing, the scalar conservation law $u_t+f(u{)}_x=0$ has solution $u(t,X_t(\xi ))=u_0(\xi )$ with $X_t(\xi )=\xi +t{f}^{\prime }(u_0(\xi ))$. When this map is a diffeomorphism, the solution is $u_0\circ X_t^{-1}$. Its spatial derivative is $u_0^{\prime }(\xi )/(1+t{f}^{\prime \prime }(u_0(\xi ))u_0^{\prime }(\xi ))$.

Putting $y=x^{5/3}$ and $w=x^{1/3}h$ transforms $h_t+x^{-1}(xD{h}^3{)}_x=0$ into the scalar conservation lawAlong smooth characteristic curves, $w$ is constant and $y$ grows at speed $5D{w}^2$. Starting from $h_0(s)$ givesAfter characteristic crossing, the Rankine-Hugoniot condition and entropy solution select the physical front.

For smooth compactly supported data in the scalar conservation law $u_t+f(u{)}_x=0$, put $m={\min }_{\xi }{f}^{\prime \prime }(u_0(\xi ))u_0^{\prime }(\xi )$. Its smooth characteristic solution of a scalar conservation law exists until $t_{\ast }=-1/m$ if $m<0$, and for all times if $m\ge 0$. The characteristic flow map $X_t(\xi )=\xi +t{f}^{\prime }(u_0(\xi ))$ has derivative $1+t{f}^{\prime \prime }(u_0)u_0^{\prime }$; its first zero produces gradient blowup. The formula applies to smooth fluxes of either concavity, although it is particularly useful for a concave function flux.

An entropy solution is a weak solution of a nonlinear scalar conservation law that also satisfies entropy inequalities selecting physically admissible shocks. For convex fluxes, this excludes expansion shocks and restores uniqueness for bounded initial data.

The finite volume method follows the amount of a conserved quantity inside each control volume. Its unknowns are cell averages, rather than values at isolated points: a cell's width times its stored average approximates that cell's total. A numerical flux transfers an amount across a face. Using one shared transfer for both neighboring cells makes the loss from one exactly the gain of the other, as explained in shared face flux cancellation.^[ref]

For a one-dimensional scalar conservation law without sources, let $u_t+{\partial }_{x}f(u)=0$, with conservation law flux $f$. Fix cells $I_j=[x_{j-1/2},x_{j+1/2}]$ of widths $h_j>0$. Define the exact cell average and its numerical approximation byThe integral balance over one cell and one time interval of length $k>0$ isFor a smooth solution this follows by integration of the differential equation; with jump discontinuities, it expresses the conservation law in integral form, using appropriate face traces when they exist. Replacing each time-averaged face conservation law flux by a numerical flux givesThus the finite volume method implements a balance directly, without needing a pointwise derivative at a jump. The tube balance in Riemann Problems and Jupyter Solutions by David I. Ketcheson, Randall J. LeVeque, and Mauricio J. del Razo (Advection, §Conservation of mass) motivates this construction from conservation of mass.

This is useful for shock waves: their speed is constrained by the Rankine-Hugoniot condition, itself an integral balance. A nonconservative approximation can give a visibly plausible jump with the wrong speed. Conservation is one requirement; a nonlinear scalar conservation law also needs the correct entropy solution, and the chosen numerical flux and update need their own stability and accuracy analysis.^[ref] The original periodic upwind finite volume advection example isolates the shared-transfer mechanism in a problem with a known translated solution.

For a scalar conservation law whose flux is a convex function, the Godunov numerical flux is the physical flux at the origin of the entropy solution of the Riemann problem between states $a,b$. It is ${\min }_{u\in [a,b]}f(u)$ for $a\le b$ and ${\max }_{u\in [b,a]}f(u)$ for $a>b$. For a rarefaction wave, the characteristic speed is monotone across its states, so the origin sees either an endpoint or a zero-speed minimizing state. For a shock wave, its speed is $(f(a)-f(b))/(a-b)$; the origin sees the left state if this speed is positive and the right state if negative, selecting the larger endpoint flux. Unlike the Godunov flux, the Engquist-Osher flux at reversed transsonic states may equal $f(a)+f(b)-f(s_{\ast })$.

For hindered settling with normalized particle volume fraction $s$, upward coordinate and downward isolated settling velocity magnitude $V_0$, the scalar conservation law has conservation law flux $f(s)=-V_{0}s(1-s{)}^{\alpha }$. For $\alpha >0$, its convexity changes at $s=2/(\alpha +1)$ when this lies in the physical interval. This change distinguishes spreading rarefaction waves from compressing characteristic curves; a flux with an inflection can produce both a shock wave and a rarefaction wave in one transition. At $\alpha =0$ the flux is linear and the profile simply translates.

The Inviscid Burgers equation is the scalar conservation lawFor a classical solution, the chain rule makes this equivalent to $u_t+u{u}_x=0$, a nonlinear transport equation with characteristic speed $u$. The conservative form remains meaningful for a discontinuous weak solution, whereas the product $u{u}_x$ cannot simply be evaluated pointwise at a jump. Its discontinuities obey the Rankine-Hugoniot condition $s=(u_L+u_R)/2$; this condition alone permits inadmissible expansion jumps. For bounded initial conditions on the infinite line, imposing the Kruzhkov entropy inequalities and the initial trace selects a unique entropy solution. For step data, incoming characteristic curves select $u_L>s>u_R$, while $u_L<u_R$ produces a rarefaction wave.^[ref]

A one-dimensional particle suspension settling without particle diffusion or overturning obeys the scalar conservation law ${\varphi }_t+{\partial }_{z}j(\varphi )=0$. The upward flux is negative for downward settling. Characteristic curves, compressive shocks and rarefaction waves describe concentration transport rather than individual particle paths.

A kinematic wave is a disturbance whose propagation speed is determined by the local relation between a conserved density and its flux. In a scalar conservation law $u_t+f(u{)}_x=0$, smooth values travel at speed $f^{\prime }(u)$.

A numerical flux estimates the rate of transfer across a cell face from nearby states. In the finite volume method update, $F_{j+1/2}^n$ approximates the conservation law flux averaged over the time step; $k{F}_{j+1/2}^n$ is the transferred amount. For a two-state numerical flux $F(L,R)$ approximating a scalar conservation law with conservation law flux $f$, consistency of a numerical method requires $F(U,U)=f(U)$. Thus equal states recover the physical transfer rate. With constant equal-width cells, this also preserves a constant state because its incoming and outgoing numerical fluxes agree.^[ref]

Consistency and shared face flux cancellation answer different questions: matching $f$ at equal states approximates the intended equation, while sharing one oriented face value preserves the total. The example periodic upwind finite volume advection uses $F(L,R)=aL$ for $a>0$, so the arriving state comes from the left. The Godunov numerical flux uses an entropy solution of a local Riemann problem; the Engquist-Osher flux provides a different construction for a scalar conservation law. A monotone conservative scheme uses a numerical flux nondecreasing in its left state and nonincreasing in its right state, together with a time step that makes the complete update monotone. Neither two-state consistency nor conservation by itself supplies that time-step restriction or a general stability or order of a numerical method guarantee.

Take $z=0$ at the base and assume gas leaves freely at the liquid surface, without being retained as foam. There is no new bubble production after filling. Neglecting the liquid velocity, gas-volume conservation is a scalar conservation lawAn impermeable base has no incoming bubbles. Thus a clear region develops below the uniform bubbly region. The clearing discontinuity satisfies the Rankine-Hugoniot condition, with states $0$ and ${\varphi }_0$:For the dilute regime $0<{\varphi }_0<1/2$, the characteristic speeds obey $j^{\prime }(0)>s>j^{\prime }({\varphi }_0)$, so this is the admissible compressive shock. Before it reaches the surface,The remaining bubbles move upwards at $s$ and escape at the top. Averaging the void fraction over the entire cylindrical volume givesThis bubble clearing in a cylindrical vessel describes bubbles in the liquid, not a hypothetical retained gas layer above it. A closed surface or long-lived foam would require a different upper boundary condition.

Hindered settling is the reduction of the particle settling velocity by displaced-fluid backflow and interactions with other particles. In the specified model, the downward speed is $V_s$ at vanishing concentration and falls to zero at the maximum packing fraction ${\varphi }_{\max }$. The characteristic speed of a concentration disturbance will differ from this individual-particle speed.

Take $z$ positive upwards. For a locally uniform horizontal suspension without diffusion or bulk fluid motion, the upward particle flux isParticle volume conservation gives the scalar conservation lawThis is the quadratic hindered-settling flux. Along its characteristic curves,For a smooth initial profile ${\varphi }_0(\zeta )$, the map is $z=\zeta +c({\varphi }_0(\zeta ))t$. Its Jacobian is $1+(2V_s/{\varphi }_{\max }){\varphi }_0^{\prime }(\zeta )t$. A shock develops where this first reaches zero: the compressive condition is ${\varphi }_0^{\prime }<0$, andwhen the minimum is negative and no earlier boundary discontinuity is imposed. A discontinuous compressive initial state can instead contain a shock from $t=0$.

Let $z=s(t)$ separate states ${\varphi }_L$ below the shock and ${\varphi }_R$ above it. Integrating the conservation law across a moving thin interval gives the Rankine-Hugoniot conditionEquivalently, the flux relative to the shock agrees on both sides. For this quadratic flux,Since $j^{\prime \prime }=2V_s/{\varphi }_{\max }>0$, a physical compression shock has ${\varphi }_L>{\varphi }_R$ and satisfies $c({\varphi }_L)>\dot{s}>c({\varphi }_R)$: characteristic curves enter it from both sides. These conditions select the shock over a nonphysical expansive discontinuity.

A stability analysis asks whether perturbations in an evolution remain bounded on a fixed physical time interval by a constant independent of the mesh. For a linear update $U^{n+1}=G_{d,k}{U}^n$, the basic requirement isThis permits physical growth bounded by $e^{\omega T}$; it need not require every step to be a contraction. Consistency of a numerical method measures how well its stencil reproduces the differential equation, whereas stability of a numerical method controls the accumulated defects. For a well-posed linear problem, the Lax equivalence theorem connects stability and convergence for a consistent approximation in the specified norm.

On the whole line or a periodic grid, constant-coefficient stencils commute with translations. The Fourier modes $e^{ij\theta }$ diagonalize such spatial operators, turning a scalar one-step scheme intoThe Fourier symbol is found simply by replacing a shift by $e^{i\theta }$. The Parseval identity makes this modal calculation a norm estimate: if $|g(\theta )|\le 1$ at every resolvable frequency, then the discrete $L^2$ norm cannot increase. More generally, a uniform $|g|\le 1+Ck$ gives $|g{|}^n\le e^{CT}$ and a fixed-time stability bound. All frequencies matter, including those near the grid scale; a small-wavenumber expansion establishes consistency but cannot test the entire stability range.

For the heat equation, put $r=k/d^2$ and $s={\sin }^2(\theta /2)$. The Forward Euler diffusion scheme hasThe condition $|g_{\mathrm{FE}}|\le 1$ for every $s\in [0,1]$ is exactly $0\le r\le 1/2$. This illustrates a parabolic time-step restriction $k\le d^2/2$. The Backward Euler diffusion scheme instead haswhich is contractive for every $r\ge 0$. The Crank-Nicolson diffusion scheme hasand is likewise unconditionally $L^2$ stable. However, its high-frequency amplification approaches $-1$ for large $r$, so it need not damp rapidly decaying modes effectively and may produce alternating transients from nonsmooth data. An A-stable time integrator can stabilize a diffusion semidiscretization with left-half-plane spectrum, but damping and accuracy remain separate questions.

For transport $u_t+a{u}_x=0$, define the hyperbolic Courant number $\nu =ak/d$. Forward Euler with a centered spatial difference gives $g=1-i\nu \sin \theta$ and $|g{|}^2=1+{\nu }^2{\sin }^2\theta$. It is unstable under a fixed nonzero CFL ratio. For $a>0$, the upwind finite difference scheme hasThus $0\le \nu \le 1$ is the stable range. Upwinding for negative $a$ reverses the spatial difference and gives $|\nu |\le 1$. The Lax-Wendroff advection scheme gives a second-order example:It is stable for $|\nu |\le 1$, but that $L^2$ result does not imply monotonicity or exclude oscillations at shocks. The Courant–Friedrichs–Lewy condition has a domain-of-dependence interpretation for explicit hyperbolic schemes; the Fourier calculation supplies the sharper algebraic stability criterion for the particular stencil.

For multilevel methods, each frequency gives an amplification polynomial and a companion matrix. Every amplification root must lie in the closed unit disk, with unit-modulus roots simple. Uniformity across frequencies and meshes is essential: one also needs control of the associated matrix powers, for example through a uniformly bounded diagonalization or the uniform power bound from separated amplification roots. A repeated unit root generates a Jordan block and polynomial growth in the number of steps. The same issue occurs for systems: checking only eigenvalue moduli of a nonnormal Fourier amplification matrix is insufficient. For instancehas unit eigenvalues but no fixed-time mesh-uniform bound as $n=T/k\to \infty$. For a hyperbolic system with a constant coefficient matrix admitting a uniformly bounded characteristic diagonalization, scalar modal estimates for each characteristic speed do give a system estimate. A positive energy symmetrizer provides another route. Frozen-coefficient Fourier tests for variable coefficients are useful local diagnostics, but a global estimate must also account for coefficient variation and boundary terms.

Boundary conditions change the allowable modes and can change the conclusion. A periodic calculation is exact only for a periodic closure. On a finite interval with zero Dirichlet boundary conditions, the centered diffusion matrix is diagonalized instead by the discrete sine transform. On $M$ interior points its eigenvalues are $-4{\sin }^2(j\pi /[2(M+1)])$ before scaling by $d^{-2}$. Thus the exact fixed-grid Forward Euler restriction iswhich tends to the whole-line threshold $1/2$ as the grid is refined. Backward Euler and Crank-Nicolson remain stable for all nonnegative $r$. An appropriate Neumann boundary condition uses cosine-type modes and includes a constant mode with eigenvalue zero. The constant mode must be preserved, rather than forced to decay; depending on the boundary stencil, the natural discrete inner product may have endpoint weights.

For hyperbolic transport, the continuous and discrete boundary closures must respect the direction of propagation. If $a>0$ on $[0,L]$, prescribe the left inflow value; prescribing an independent outflow value can overdetermine the problem. The energy method exposes this directly:Homogeneous inflow is dissipative, while forced inflow contributes data to the estimate. A compatible upwind boundary update has the same one-sided information flow. For diffusion,The last term vanishes for homogeneous Dirichlet or Neumann conditions. With Robin boundary conditions of dissipative sign it contributes a further nonpositive term. A non-dissipative boundary can produce genuine PDE growth, which a valid stability estimate should represent rather than incorrectly suppress.

Even a perfectly stable interior symbol cannot detect every defective boundary closure of a difference scheme. As a concrete counterexample, leave a stable explicit heat stencil on all grid points away from the left endpoint, but replace its first interior row by $U_1^{n+1}=2U_1^n$, still keeping $U_0=0$. The interior Fourier symbol is unchanged, yet $U_1^n=2^n{U}_1^0$ makes the finite-interval update unstable. The faulty boundary-adjacent row creates a mode absent from the periodic interior calculation. Thus the entire finite update matrix, or an energy estimate including its boundary rows, must be examined.

A more systematic boundary stability of a finite-difference method analysis transforms time and any tangential spatial variables, then looks for modes $U_j^n=z^n{\kappa }^j$ with $|z|>1$ and $|\kappa |<1$ on a half-line. The interior dispersion relation selects decaying spatial modes; the boundary equations must determine their coefficients without admitting a nonzero growing homogeneous mode. A uniform inverse estimate for the boundary system, expressed by the Uniform Kreiss--Lopatinskii condition, also rules out arbitrarily weak boundary control near limiting frequencies. Merely excluding one obvious unstable mode is weaker than proving a mesh-uniform estimate. Discrete summation by parts and suitable boundary penalties offer an energy-based alternative.

Finally, nonlinear scalar conservation laws can develop discontinuities, so linearized Fourier stability is only one part of the analysis. A monotone conservative scheme, such as the Engquist-Osher method under its CFL restriction, has an invariant state range and L1 contraction of a monotone conservative scheme, with a discrete entropy inequality selecting admissible shocks. These complement the Fourier and boundary estimates. For an evolution problem, the meaningful outcome is a stability bound in a specified norm, uniform over the chosen mesh family and compatible with the actual boundary and initial data.

For the negative-flux Inviscid Burgers equation, the method of characteristics givesThe characteristic curve starting at ${\theta }_0$ therefore has $\theta ={\theta }_0-f_0({\theta }_0)z$, and henceFor smooth data this describes a single-valued classical solution as long as the characteristic flow map is invertible. Its Jacobian is $1-z{f}_0^{\prime }({\theta }_0)$; after characteristic crossing, one must instead select a weak solution with the appropriate entropy solution condition.

The conservation form is $f_z+G(f{)}_{\theta }=0$, with conservation law flux $G(f)=-f^2/2$. Integrating this scalar conservation law across a moving discontinuity, or differentiating its Heaviside representation as a distribution, yields the Rankine-Hugoniot conditionFor distinct one-sided limits it simplifies toThis jump-speed relation is exact for the Burgers conservation law; an additional entropy condition is needed to distinguish a physical compressive shock wave from an expansion discontinuity.

The step gives the Burgers Riemann problem with negative flux. If $U>0$, characteristic curves from the left have speed zero, while those from the right have speed $-U$: they converge. The entropy shock wave has speed $-U/2$ and thusCharacteristics enter this shock from both sides, since $0>-U/2>-U$.

If $U<0$, the right-hand speed $-U$ is positive and the two families separate. A smooth steep approximation to the initial step spreads into a rarefaction wave. In the fan, the self-similar characteristic curve relation is $-f=\theta /z$, givingThe endpoint values match continuously. The discontinuity moving at $-U/2$ would satisfy the jump condition even for $U<0$, but its characteristic curves leave the discontinuity and it fails entropy admissibility. For $U=0$ the solution is identically zero.

For the sharp energy estimate, $C_0=0$ stays uniform while the available bound $C_1=M^2/\epsilon$ diverges as $\epsilon \downarrow 0$ when $M>0$. If the coarser estimate is used for the first part, both displayed bounds $M^2/\epsilon$ diverge, but the first divergence is only an artifact of discarding an exact cancellation.

The method of characteristics explains why uniform control of the gradient cannot generally persist for the inviscid scalar conservation law. Before characteristic crossing, with initial data $u_0$,If $F^{\prime \prime }(u_0(\xi ))u_0^{\prime }(\xi )<0$ somewhere, the denominator reaches zero in finite positive time. The solution steepens and the smooth description breaks down; an entropy solution can subsequently contain shocks. Positive viscosity replaces such a discontinuity by a thin smooth layer, which can have large gradient even while its $L^2$ norm remains controlled. The divergent bound does not assert that every flux and every initial datum form a shock; a linear flux, for example, has no such steepening.

Fix $u_l>u_r$ and the Rankine-Hugoniot condition speed from the preceding part. The new hypothesis is a uniformly convex scalar flux, $F^{\prime \prime }\ge \kappa >0$; it replaces the globally bounded-$F^{\prime }$ hypothesis of part (a). Choose $b=\sigma u_r-F(u_r)$. Thensince a strictly convex function lies below the chord between its two endpoint values. Also $Q(u_l)=Q(u_r)=0$ andThe zeros at the endpoints are simple, so the separated integral diverges logarithmically there. Thus the travelling wave is a decreasing connection defined for all $s$, unique up to translation.

To specify a limit, fix a number $c\in (u_r,u_l)$ independently of $\epsilon$ and normalize $v_{\epsilon }(0)=c$. If $V^{\prime }=Q(V)$ and $V(0)=c$, uniqueness gives $v_{\epsilon }(s)=V(s/\epsilon )$. ConsequentlyAt $x=\sigma t$ the normalized profile equals $c$ for every $\epsilon$. This single-line value is immaterial to the weak solution. Convergence holds pointwise off the line and in $L_{\mathrm{loc}}^1$ by bounded convergence; it cannot be uniform across a nonzero jump. The transition has thickness of order $\epsilon$.

This is the vanishing viscosity approximation to a compressive entropy shock. The Rankine-Hugoniot condition makes the step a weak solution of the inviscid scalar conservation law, while $F^{\prime }(u_r)<\sigma <F^{\prime }(u_l)$ means characteristic curves enter the shock from both sides. For every smooth convex function $\eta$ used as an entropy, with entropy flux for a scalar conservation law $q^{\prime }={\eta }^{\prime }{F}^{\prime }$, the viscous equation givesAgainst compactly supported tests the right side tends to zero, since $u_{\epsilon }$ stays in $[u_r,u_l]$. Passing to the $L_{\mathrm{loc}}^1$ limit yields the entropy inequality, explaining the direction selected by positive viscosity.

A translation must be fixed to obtain this particular limit. An $\epsilon$-dependent translate $V((s-a_{\epsilon })/\epsilon )$ can converge to a shock at a different location, to a constant if its center escapes, or fail to converge if the centers oscillate. Thus existence of profiles alone does not specify a single vanishing-viscosity limit without a phase normalization.