Method of Characteristics

• Method of Characteristics is a technique that reduces first-order PDEs to ODEs along characteristic curves, enabling clear analysis of transport dynamics.
• It is applied across scalar equations, nonlinear systems, and multidimensional problems, including fluid dynamics, neutron transport, and population balances.
• Combined with numerical schemes like finite volume and spline interpolation, it delivers scalable and stable solutions even in shock and entropy-dissipative scenarios.

The method of characteristics is a classical method for gaining understanding in the solution of a partial differential equation, and it was used to solve the initial value problem for partial differential equations of first order. In its standard form, the method replaces a PDE by ordinary differential equations posed along characteristic curves, so that the solution propagates along paths determined by the transport structure of the equation. Across the literature represented here, that basic reduction is used in scalar first-order equations, nonlinear hyperbolic systems, population balance equations, neutron transport, radiation transport, blood flow, magnetohydrodynamics, and several generalized settings including fractional and partially integrable systems (Wu, 2010, Ponomarev et al., 2019, Li et al., 2012).

1. Classical reduction along characteristic curves

For a first-order PDE

$$a(x,t)\frac{\partial u}{\partial x} + b(x,t)\frac{\partial u}{\partial t} = c(x,t),$$

the classical characteristic construction seeks curves $(x(s),t(s))$ such that

$$\frac{dx}{ds} = a(x,t), \qquad \frac{dt}{ds} = b(x,t), \qquad \frac{du}{ds} = c(x,t).$$

In this form, the PDE is reduced to an ODE system along the characteristic parameter $s$ (Wu, 2010).

A quasilinear formulation used for online simulation is

$$w_t + v\big( t, x, w(t, \cdot)_{[0, \ell]} \big) w_x = f\big( t, x, w(t, \cdot)_{[0, \ell]} \big),$$

with boundary data $w(t,0)=u(t)$. The corresponding characteristic equations are

$$\dot\xi(t) = v\big( t, \xi(t), w(t, \cdot)_{[0, \ell]} \big), \qquad \dot\omega(t) = f\big( t, \xi(t), w(t, \cdot)_{[0, \ell]} \big),$$

and the PDE solution at $(t,x)$ is equal to $\omega(t)$ where $x=\xi(t)$ (Ponomarev et al., 2019). This formulation makes explicit that characteristic curves are not merely geometric objects: they are the trajectories along which the state is updated.

The same principle appears in population balance equations. For

$(x(s),t(s))$0

the method is applied to the hyperbolic part

$(x(s),t(s))$1

with characteristic relation

$(x(s),t(s))$2

A first-order discretization along characteristics yields

$(x(s),t(s))$3

where $(x(s),t(s))$4 and $(x(s),t(s))$5 is obtained by linear interpolation (Li et al., 2012). In that setting, the higher dimensional population balance equation is transformed into a series of lower dimensional convection-diffusion-reaction equations in the physical space.

2. Systems, second-order equations, and generalized characteristic structures

For the hyperbolic Monge-Ampère equation

$(x(s),t(s))$6

three alternative derivations of the method of characteristics lead to the same characteristic condition. Writing $(x(s),t(s))$7, $(x(s),t(s))$8, $(x(s),t(s))$9, the characteristic equation is

$$\frac{dx}{ds} = a(x,t), \qquad \frac{dt}{ds} = b(x,t), \qquad \frac{du}{ds} = c(x,t).$$0

with discriminant

$$\frac{dx}{ds} = a(x,t), \qquad \frac{dt}{ds} = b(x,t), \qquad \frac{du}{ds} = c(x,t).$$1

so there are two real, distinct roots,

$$\frac{dx}{ds} = a(x,t), \qquad \frac{dt}{ds} = b(x,t), \qquad \frac{du}{ds} = c(x,t).$$2

The method gives rise to two mutually coupled systems of ordinary differential equations, and the paper solves them with explicit one-step methods (Euler, Runge-Kutta) and spline interpolation (Bertens et al., 2021).

For two-dimensional compressible Euler flows, the method is used for both the direct and the adjoint problem. The characteristic directions follow from the singularity of the block linear system associated with the Cauchy problem, with streamlines $$\frac{dx}{ds} = a(x,t), \qquad \frac{dt}{ds} = b(x,t), \qquad \frac{du}{ds} = c(x,t).$$3 defined by $$\frac{dx}{ds} = a(x,t), \qquad \frac{dt}{ds} = b(x,t), \qquad \frac{du}{ds} = c(x,t).$$4 and characteristic families $$\frac{dx}{ds} = a(x,t), \qquad \frac{dt}{ds} = b(x,t), \qquad \frac{du}{ds} = c(x,t).$$5 defined by $$\frac{dx}{ds} = a(x,t), \qquad \frac{dt}{ds} = b(x,t), \qquad \frac{du}{ds} = c(x,t).$$6. A linear algebra analysis shows that along $$\frac{dx}{ds} = a(x,t), \qquad \frac{dt}{ds} = b(x,t), \qquad \frac{du}{ds} = c(x,t).$$7 there is only one independent characteristic equation for each family, while along streamlines $$\frac{dx}{ds} = a(x,t), \qquad \frac{dt}{ds} = b(x,t), \qquad \frac{du}{ds} = c(x,t).$$8 there are two linearly independent characteristic equations (Ancourt et al., 2023).

A constructive generalization to systems of PDE considers

$$\frac{dx}{ds} = a(x,t), \qquad \frac{dt}{ds} = b(x,t), \qquad \frac{du}{ds} = c(x,t).$$9

with Lipschitz or $s$0 coefficients and initial data on a hyperplane. The solution is constructed by discretizing the domain of the solution along the direction perpendicular to the initial condition hyperplane, defining upper and lower bounds on successive hyperplanes, and letting the number of discretization hyperplanes tend to infinity. In that limit, the upper and lower bounds approach each other, yielding a unique function for the solution, denoted $s$1, and producing a local unique Lipschitz or $s$2 solution (Shalchian, 2018).

The fractional method of characteristics replaces integer-order derivatives and differentials by modified Riemann-Liouville derivatives and fractional differentials. For

$s$3

the fractional characteristic equations are

$s$4

This is presented as a fractional generalization of the classical method (Wu, 2010).

A separate line of generalization reduces $s$5-dimensional nonlinear PDEs to systems of $s$6-dimensional first order PDEs, integrating the original system with respect to only one independent variable and thereby producing partially integrable systems (Zenchuk, 2010). A related construction combines the characteristics method and Hopf-Cole transformation to partially integrate a large class of multidimensional nonlinear PDEs, including a specific generalization of the equation describing the dynamics of two-dimensional viscous fluid and a generalization of the Korteweg-de Vries equation (Zenchuk, 2012).

3. Numerical realizations and algorithmic patterns

In numerical practice, the method of characteristics is often coupled to interpolation, finite volume or finite element discretization, and explicit or semi-implicit ODE solvers. In population balance equations, the characteristic step in the $s$7 variables is followed by finite element discretization in the physical space, producing a family of convection-diffusion-reaction equations indexed by the internal coordinate (Li et al., 2012).

For first-order quasilinear PDEs in an online simulation environment, a dynamic pool of characteristics is tracked by the vectors $s$8 and $s$9. Simultaneous integration of all current characteristics is carried out with Simulink’s built-in solver, while Simulink’s zero-crossing detection algorithm is used to create new characteristics dynamically and discard the old ones. New characteristics are inserted at the inlet when the spatial gap to existing characteristics is at least $$w_t + v\big( t, x, w(t, \cdot)_{[0, \ell]} \big) w_x = f\big( t, x, w(t, \cdot)_{[0, \ell]} \big),$$0, the value gap to input is at least $$w_t + v\big( t, x, w(t, \cdot)_{[0, \ell]} \big) w_x = f\big( t, x, w(t, \cdot)_{[0, \ell]} \big),$$1, or after a fixed time interval at least $$w_t + v\big( t, x, w(t, \cdot)_{[0, \ell]} \big) w_x = f\big( t, x, w(t, \cdot)_{[0, \ell]} \big),$$2; characteristics are removed once they have traveled beyond $$w_t + v\big( t, x, w(t, \cdot)_{[0, \ell]} \big) w_x = f\big( t, x, w(t, \cdot)_{[0, \ell]} \big),$$3. The implementation is bundled in the FOQL_PDE block (Ponomarev et al., 2019).

For one-dimensional blood flow, the numerical method of characteristics is formulated in Riemann invariant variables

$$w_t + v\big( t, x, w(t, \cdot)_{[0, \ell]} \big) w_x = f\big( t, x, w(t, \cdot)_{[0, \ell]} \big),$$4

with characteristic speeds

$$w_t + v\big( t, x, w(t, \cdot)_{[0, \ell]} \big) w_x = f\big( t, x, w(t, \cdot)_{[0, \ell]} \big),$$5

At each time step, the method estimates the foot of the characteristic explicitly,

$$w_t + v\big( t, x, w(t, \cdot)_{[0, \ell]} \big) w_x = f\big( t, x, w(t, \cdot)_{[0, \ell]} \big),$$6

interpolates the previous solution at that foot using piecewise linear interpolation, and updates the source term by explicit Euler. Boundary conditions are formulated using characteristics, and network branching conditions are enforced through nonlinear algebraic equations solved by Newton’s method (Acosta et al., 2014).

In unstructured-mesh shallow water computation, the finite volume characteristics scheme is a predictor-corrector procedure in which numerical fluxes are reconstructed using the method of characteristics while an Eulerian finite volume method discretizes the conservation law. For each edge midpoint, the method backtracks along the characteristic to locate the departure point, interpolates the state there, updates along the projected characteristic with source terms, transforms projected velocities back to Cartesian components, and then computes the numerical flux. The method is formulated on unstructured triangular meshes and avoids the need to solve a Riemann problem (Ziggaf et al., 2021).

For the hyperbolic Monge-Ampère equation, explicit one-step methods are paired with local B-spline interpolation because numerical characteristics do not, in general, land on the grid. Forward Euler is paired with linear splines, Modified Euler with quadratic splines, and RK4 with quartic splines (Bertens et al., 2021).

4. Shocks, entropy, and numerical stability

The method of characteristics continues to encode the solution structure after shock formation, but classical single-valued reconstruction fails once characteristics overlap. A convex hull formulation addresses this issue for nonlinear hyperbolic conservation laws and Hamilton-Jacobi equations. From the multivalued solutions determined by the method of characteristic, the convex hull algorithm extracts the entropy dissipative solutions, even after the formation of shocks, and it also computes entropy conservative, or multivalued, solutions. The construction applies to both convex or non-convex flux/Hamiltonians (LeFloch et al., 2014).

Stability depends strongly on how the ODEs along characteristics are integrated. For non-dissipative hyperbolic PDEs with periodic boundary conditions, the method of characteristics coupled with simple Euler, modified Euler, and Leap-frog was analyzed in detail. All three methods are unstable for periodic boundary conditions; simple Euler has mild instability at the lowest and highest wavenumbers, modified Euler has strongest instability at the middle of the spectrum, and Leap-frog has by far the strongest instability among all three solvers. The same study comments that the fourth-order Runge–Kutta solver within the MoC framework is also unstable in the tested construction (Lakoba et al., 2016).

A closely related analysis shows that changing periodic boundary conditions to nonreflecting boundary conditions can alter stability in a decisive way. Non-periodic boundary conditions may render some of the MoC schemes stable for most practical computations, even though they are unstable for periodic boundary conditions. Specifically, MoC-SE and MoC-ME become practically stable for most modes, while MoC-LF remains unstable, with nearly the same growth rates as for periodic boundary conditions. This contradicts a statement found in some literature that an instability detected by the von Neumann analysis for a given numerical scheme implies an instability of that scheme with arbitrary boundary conditions (Lakoba et al., 2016).

The mechanism identified in that analysis is structural rather than incidental. With nonreflecting boundary conditions, error modes become monotonically spatially growing or decaying exponentials; there are no bulk modes resembling Fourier harmonics. By contrast, for MoC-LF, parasitic root multipliers coincide with true root multipliers at mid-spectrum, so unstable Fourier-like modes continue to exist (Lakoba et al., 2016). This is a precise correction to a common misconception about the transfer of periodic-boundary stability conclusions to non-periodic settings.

5. Representative application domains

The modern literature uses the method of characteristics across a wide range of transport-dominated models. The following cases are explicitly documented in the cited works.

Domain	Characteristic role	Representative paper
Population balance equations	transforms the higher dimensional equation into a series of lower dimensional convection-diffusion-reaction equations	(Li et al., 2012)
Neutron transport	integrates along straight particle tracks; CAD-based geometry and ray tracing via MCAM; time-dependent spherical transport avoids negative flux problem	(Chen et al., 2013, Barroso, 2017)
Blood flow	diagonal nonlinear hyperbolic system solved in characteristic variables with unconditional stability	(Acosta et al., 2014)
MHD and resistive RMHD	Alfvén waves handled by the method of characteristics; electromagnetic field updated with appropriate characteristic velocities	(Iwasaki et al., 2011, Takamoto et al., 2011)
Radiation magnetohydrodynamics	invariant specific intensity transported along null geodesics; moments closed from a swarm of samples	(Ryan et al., 2019)
Hybrid plasma simulation	Vlasov equation solved by following collisionless phase point trajectories in phase space	(Javaheri et al., 2015)

In neutron transport, the method is widely used for neutron transport calculation, but highly heterogeneous geometry makes geometry modeling and ray tracing a central difficulty. A CAD-based automatic modeling tool, MCAM, was introduced for geometry modeling and ray tracing, and a new MOC code was developed and integrated into the SuperMC system. The diamond-difference scheme was applied to MOC to reduce the spatial discretization errors of the flat flux approximation (Chen et al., 2013). A separate time-dependent neutron transport formulation in one-dimensional spherical symmetry uses the characteristic method to overcome the negative flux problem that arises when the system is very opaque and the angular neutron flux can become negative when it is extrapolated in spatial meshes; an additional term is required to maintain neutron conservation in a moving system in Lagrangian coordinates (Barroso, 2017).

In smoothed particle magnetohydrodynamics, a nonlinear Riemann problem with magnetic pressure is solved for compressive waves, while the method of characteristics is applied for Alfvén waves. The relevant perpendicular variables satisfy

$$w_t + v\big( t, x, w(t, \cdot)_{[0, \ell]} \big) w_x = f\big( t, x, w(t, \cdot)_{[0, \ell]} \big),$$7

and the interface state is recovered from the characteristic invariants (Iwasaki et al., 2011). In resistive relativistic MHD, the numerical flux of the fluid is calculated by using an approximate Riemann solver, and the electromagnetic field is updated by using the method of characteristics with appropriate characteristic velocities rather than the speed of light in all regimes (Takamoto et al., 2011).

In covariant radiation magnetohydrodynamics, the radiative transfer equation is solved in invariant form deterministically for each sample, each sample transports a multigroup representation of the invariant specific intensity along a null geodesic, and the samples in each zone are used to construct a Delaunay triangulation of the unit sphere in the comoving frame to calculate the Eddington tensor and average source terms (Ryan et al., 2019). In the hybrid kinetic-fluid plasma scheme, electrons are treated kinetically by following characteristic phase-point trajectories in phase space, while ions obey the fluid equations; the characteristic equations and the fluid equations are both solved using the Leapfrog-Trapezoidal method, initialized by the Euler-Trapezoidal scheme (Javaheri et al., 2015).

6. Accuracy, scalability, and computational consequences

The method of characteristics is frequently adopted because it reorganizes transport-dominated dynamics into update rules aligned with propagation directions. In population balance equations, this leads directly to parallelism: after characteristic reduction, the equations for different internal grid points are independent at each time step except for a minor sequential dependency for constructing $$w_t + v\big( t, x, w(t, \cdot)_{[0, \ell]} \big) w_x = f\big( t, x, w(t, \cdot)_{[0, \ell]} \big),$$8, so the bulk of the workload per time step is embarrassingly parallel. Strong scaling results report computation time decreasing from $$w_t + v\big( t, x, w(t, \cdot)_{[0, \ell]} \big) w_x = f\big( t, x, w(t, \cdot)_{[0, \ell]} \big),$$9 s on 1 processor to $w(t,0)=u(t)$0 s on 16 processors, with speedup increasing from $w(t,0)=u(t)$1 to $w(t,0)=u(t)$2; weak scaling tables show computation times remaining nearly constant (Li et al., 2012).

The blood-flow numerical method of characteristics is explicitly described as efficient and unconditionally stable. Because no CFL time-step restriction is needed, it can use $w(t,0)=u(t)$3 up to 12x larger than allowed by explicit discontinuous Galerkin for stability in a network of 64 arteries, and the reported speedups are $w(t,0)=u(t)$4 compared with dG, with comparable accuracy in physiologically relevant scenarios (Acosta et al., 2014). A plausible implication is that characteristic-based marching becomes especially attractive when wave speeds are much larger than advective velocities and explicit Eulerian schemes are dominated by stability rather than accuracy constraints.

In CAD-based neutron transport, benchmark accuracy is reported numerically: the relative error of $w(t,0)=u(t)$5 between SuperMC and reference is only $w(t,0)=u(t)$6 for the ISSA 1D problem; agreement for the multi-cell full assembly is within $w(t,0)=u(t)$7; and C5G7 MOX benchmark results are within $w(t,0)=u(t)$8 (Chen et al., 2013). In the time-dependent spherical transport formulation, the consistency of the numerical solution is checked by making at each time-step the balance of neutrons in the system, and the reported neutron conservation error is less than $w(t,0)=u(t)$9 in one heterogeneous pellet simulation and less than $$\dot\xi(t) = v\big( t, \xi(t), w(t, \cdot)_{[0, \ell]} \big), \qquad \dot\omega(t) = f\big( t, \xi(t), w(t, \cdot)_{[0, \ell]} \big),$$0 in the multi-layered hybrid fission-fusion-fission system (Barroso, 2017).

For the hyperbolic Monge-Ampère equation, the documented convergence behavior follows the time integrator: Forward Euler with linear spline gives first-order convergence, Modified Euler with quadratic spline gives second-order convergence, and RK4 with quartic spline gives fourth-order convergence (Bertens et al., 2021). In the population balance setting, the method of characteristics has first-order accuracy in the internal coordinate $$\dot\xi(t) = v\big( t, \xi(t), w(t, \cdot)_{[0, \ell]} \big), \qquad \dot\omega(t) = f\big( t, \xi(t), w(t, \cdot)_{[0, \ell]} \big),$$1, while the finite element method in space with linear elements ($$\dot\xi(t) = v\big( t, \xi(t), w(t, \cdot)_{[0, \ell]} \big), \qquad \dot\omega(t) = f\big( t, \xi(t), w(t, \cdot)_{[0, \ell]} \big),$$2) shows rates close to optimal and quadratic elements ($$\dot\xi(t) = v\big( t, \xi(t), w(t, \cdot)_{[0, \ell]} \big), \qquad \dot\omega(t) = f\big( t, \xi(t), w(t, \cdot)_{[0, \ell]} \big),$$3) show higher rates as expected (Li et al., 2012).

Characteristic-based finite volume variants target a similar balance between fidelity and cost. The finite volume characteristics scheme for the shallow water equations is accurate, conservative, non-oscillatory, and avoids the need to solve a Riemann problem; it is reported to be stable for high CFL numbers, for example $$\dot\xi(t) = v\big( t, \xi(t), w(t, \cdot)_{[0, \ell]} \big), \qquad \dot\omega(t) = f\big( t, \xi(t), w(t, \cdot)_{[0, \ell]} \big),$$4, and to outperform the Roe scheme in relative $$\dot\xi(t) = v\big( t, \xi(t), w(t, \cdot)_{[0, \ell]} \big), \qquad \dot\omega(t) = f\big( t, \xi(t), w(t, \cdot)_{[0, \ell]} \big),$$5-error for dam-break problems on unstructured meshes (Ziggaf et al., 2021). A related finite volume method constructs intermediate states at interfaces using the method of characteristics, generates entropic solutions, avoids solving Riemann problems, and introduces a diffusion control parameter to increase the accuracy of the scheme; the numerical examples demonstrate the ability to capture contact discontinuity and shock wave profiles with high accuracy and low cost as well as robustness (Aboussi et al., 2023).

In radiation transport, the MOCMC method closes the radiation moment equations with a swarm of characteristic samples and reports convergence at least as $$\dot\xi(t) = v\big( t, \xi(t), w(t, \cdot)_{[0, \ell]} \big), \qquad \dot\omega(t) = f\big( t, \xi(t), w(t, \cdot)_{[0, \ell]} \big),$$6, rather than the canonical Monte Carlo $$\dot\xi(t) = v\big( t, \xi(t), w(t, \cdot)_{[0, \ell]} \big), \qquad \dot\omega(t) = f\big( t, \xi(t), w(t, \cdot)_{[0, \ell]} \big),$$7, where $$\dot\xi(t) = v\big( t, \xi(t), w(t, \cdot)_{[0, \ell]} \big), \qquad \dot\omega(t) = f\big( t, \xi(t), w(t, \cdot)_{[0, \ell]} \big),$$8 is the number of samples per zone (Ryan et al., 2019). This suggests that, in settings where the transport operator can be evaluated deterministically along characteristics while the moments are evolved conservatively, the method can combine directional resolution with a reduced noise floor relative to standard Monte Carlo sampling.