The method of characteristics allows to solve quasilinear hyperbolic equations with high accuracy by reformulating the problem into two equations that, respectively, compute the state update and track the convection along the characteristic direction. In a generic system of M quasilinear hyperbolic equations, up to M distinct characteristics directions are identified (one for each equation of the system). The reformulation of a system of quasilinear hyperbolic equations with the method of characteristics results into M systems, each one updating a state along a characteristic direction. However, for a generic quasilinear hyperbolic system, the M systems will be coupled, since the source terms of the M state update equations can depend on multiple states and, potentially, on the entire state vector. Hence, the method of characteristics cannot be applied to solve a generic quasilinear hyperbolic system. This work presents a hybrid numerical scheme that combines the method of characteristics with the method of lines to address this issue. A spatial discretization operator based on finite differences, finite elements, or spectral methods is used to reformulate a system with up to M independent characteristics to a system with only one family of characteristics. The resulting system can be solved as an initial value problem with an off-the-shelf ODE solver. We derive necessary and sufficient conditions for which the norm of the semi-discretization error of the hybrid characteristics/lines numerical methods is lower than the norm of the semi-discretization error of the method of lines that uses the same spatial discretization operator for the semi-discretization step. Global truncation error and stability are discussed through numerical examples.
