## 4 Method for Assembling Numerical Schemes
### PDE Classification  
To find a numerical scheme that produces stable results, knowing the qualitative behavior of the system is oftentimes a necessity. In particular, this means that the specific capabilities that a chosen numerical scheme possesses need to reflect the properties that the system presents. We illustrate this by example. We once again investigate the simple advection equation (Equation (1)), which is known to be first-order hyperbolic. Such systems are prone to either preserving or even amplifying discontinuities given in the initial condition; thus, the capability of accurately representing these should be incorporated into the choice of numerical schemes. Suitable candidates would then be a finite volume or discontinuous Galerkin method. However, the finite difference method using a centered stencil or the continuous Galerkin method would give suboptimal results. The strong imposition of continuity in the domain would then yield spurious oscillations that affect stability. We hence require the system of PDEs to be classified firsthand. We follow the popular taxonomy of second-order PDEs that can, for example, be found in the book by [PERSON] [47], but follow a more general method for determining the appropriate class [48]. That is, we define a singular governing equation in the form of a PDE to be either elliptic, parabolic, or hyperbolic, depending on the shape that its characteristic quadric takes in space. More formally, consider a differential operator \(L\) of the form  
\[L[u]=\sum_{i,j=1}^{n}a_{ij}\frac{\partial^{2}u}{\partial x_{i}\partial x_{j}}, \tag{17}\]  
where \(x_{i}\) are the dependent variables and \(a_{ij}\) is the matrix forming the coefficients of the highest spatial derivatives. Considering the eigenvalues \(\lambda_{i}\) of \(a_{ij}\), \(L\) is called  
* Elliptic, if all \(\lambda_{i}\) are either positive or negative;
* Parabolic, if at least one eigenvalue is zero and all others are either positive or negative;
* Hyperbolic, if at least one eigenvalue is positive and at least one is negative. The characterization of first-order differential operators is more straightforward, however. It can be shown that first-order PDEs with constant, real coefficients are always hyperbolic. This condition is met for most cases relevant to engineering or physical applications. More precisely, a first-order PDE is hyperbolic if the resulting Cauchy problem is uniquely solvable. In the case of real, constant coefficients, the polynomial equation for each variable has to admit \(n\) solutions for an equation of order \(n\) while keeping all other variables constant. In the present case, this is trivially true. We apply this classification for each governing equation of the independent variables for a given multiphysics problem. In practice, one may oftentimes identify the class by the differential operators that frequently appear in a given PDE. For example, a PDE that only has a Laplacian as a spatial differential operator--such as the Laplace equation \(\Delta u=0\) or the heat equation \(\partial_{t}u-\Delta u=0\)--exhibits dissipative behavior and is prototypical for elliptic and parabolic PDEs. Oftentimes, one can easily identify a differential operator as parabolic if it has an elliptic operator in its spatial derivatives and an additional temporal derivative, as is exactly the case for the heat equation. Both the abovementioned classes of PDE are dissipative, with the reason being that PDEs of second order can only have discontinuous derivatives along their characteristics. Since elliptic differential operators lack any characteristics, they strictly admit smooth solutions in that sense [48]. Thus, we associate this qualitatively dissipative behavior with elliptic and parabolic PDEs, as defined above. However, the advection equation (Equation (1)) only has the gradient as a spatial differential operator, representing purely convective behavior.