Many problems in physics are described by differential equations which in general can only be solved numerically. As a result one obtains an approximative solution whose error can be reduced to the cost of a higher numerical effort. For special classes of differential equations in space the finite element method is one of the most powerful tools regarding its numerical properties as well as the fact that it can be applied to arbitrary-shaped three-dimensional domains. However, following its original derivation the method can only be applied to even-ordered differential equations in space. The scope of this contribution is to broaden the finite element method towards the solution of first-order differential equations which may result directly from the underlying problem or as the state-space representation of a higher-order differential equation. One may think of the epidemic models exemplary towards the actual spreading of the new corona virus.
Due to the meaningfulness of the finite element method many commercial and non-commercial programs exist. In order to obtain a finite-element formulation of a problem which is given in terms of a differential equation, one can apply the method of weighted residuals. Thereby, the so-called weak form is deduced from a partial integration of the weighted residuum. As a consequence, one order of derivative is shifted from the field variable to the weight function. As soon as the orders of derivative of both variables coincide, a finite element ansatz in conjunction with Galerkin’s method results in symmetric element matrices. Thus, the assembled system matrices as well become symmetric which is a key property of the finite element method. Therefore, it is restricted to problems which are governed by even-ordered differential equations.
In order to overcome this restriction, one of the authors took a first-order differential equation and applied fractional partial integration of order 1/2 to the weighted residuum. As a consequence, the field variable and the weight function were operated by derivative of order 1/2. However, due to the occurrence of left and right fractional derivatives the resulting system matrix was non-symmetric and thus he failed to succeed . For this reason, in the following a different approach is applied which makes use of fractional powers of operators   . In particular, the procedure leads to a positive self-adjoint operator and hence to a symmetric system matrix. The overall goal is to establish a method that can be applied to any spatial first-order differential equation which results in conjunction with a finite element ansatz in symmetric system matrices. Since the main properties of the classical finite element method still hold with this approach the infrastructure of existing codes can be used for its implementation. A link between fractional powers of differential operators and fractional derivatives is given in .
In Section 2, a linear operator equation is derived from a general linear first-order partial differential equation and in Section 3, the polar decomposition of the resulting differential operator (applying a fractional power of order 1/2) is used to deduce a related finite element scheme with a symmetric (but dense) system matrix. The resulting method is applied to the barometric equation in Section 4. In particular, we derive the related operator formulation, determine the linear algebraic equation by an eigenvalue analysis of an associated Sturm-Liouville operator and introduce a numerical scheme to approximate the occurring integrals and solve the algebraic equation. Finally, in Section 5 we draw conclusions and give a perspective on future work.
2. Transformation of a Linear First-Order Differential Equation into an Operator Equation
In the following, we sketch the transformation of a linear first-order differential equation into an operator equation. In the case of the barometric equation, we give full detail later. Let , , , , , , . We consider a linear hyperbolic first-order partial differential equation
for , , , and initial data
which are given on the hyperplane , that is assumed to be non-characteristic. Further, let be such that
By introducing the new “unknown” function as
we obtain homogeneous initial data
and a transformed differential equation
Hence, we arrive at an operator equation
is a linear operator in and
The operator Equation (4) has a unique solution
which may be approximated as described in the next section.
3. Transformation of the Operator Equation into a Suitable form for the Application of Finite Elements
In the following, we introduce the finite element method for an approximation of the solution (5) of (4). Therefore, we use that the polar decomposition of A , which is uniquely associated with A, i.e.
where is a partial isometry with initial space and end space . Herein, denotes the range of an operator. Since A is bijective, we have
and V is a unitary transformation, i.e.
We note that (4), since,
is equivalent to the equation
The operator in X, is densely-defined, linear and positive self-adjoint and bijective. Therefore (6) has the unique solution
Now (6) can be solved by the usual finite element methods. For this purpose, let , be linearly independent elements of . Further, we denote by the orthogonal projection onto . Then
for every , where denotes the scalar product in X. In the following, we are going to solve the corresponding “approximate” system
where . The vector is an element of . Hence, there are uniquely determined such that
As a consequence,
for every . Hence, we arrive at the system of linear equations
, where the real matrix
is symmetric and positive definite. As a consequence,
for every .
4. Application to the Barometric Equation
4.1. Operator Equation
Let , , , be such that
for every and
Further, let be such that
We define a new “unknown” function by
for every . Hence, we arrive at the transformation of the system of Equations, (8) and (9), into an operator equation of the form (4), where
and denotes the derivative operator with domain
in , where (4) has the unique solution (5). In particular, if ,
for every , then
for every , implying that
4.2. Derivation of the Finite Element Formulation
where , is bijective, we need to calculate the corresponding operators
First, we note that
Hence if , where
Hence, is a densely-defined, linear and self-adjoint extension of the densely-defined, linear, symmetric and essentially self-adjoint operator , defined by
for every . Since is, in particular, closed, it follows that
and hence that
In the following, we are going to calculate and , using an approach via the theory of Sturm-Liouville operators    . For this purpose, we need the eigenvalues and eigenfunctions of . If , the solutions of
for every are given by
for every , where . The boundary condition
for every , where , and the boundary condition
gives a non-trivial u iff
for some . For such , it follows that
and hence that, if is defined by
for every , where , then
is a Hilbert basis of X. For
it follows that
and hence that
for every . Hence, it follows that
and for almost all ,
where we used that
is convergent in X, as a consequence of the estimations
, and as consequence of the existence of
We compute for satisfying
where we used the product-to-sum identity
and the series representation
found in ( , p. 1073, Formula (19)). As a consequence,
for every , where
for almost all . The function is depicted in Figure 1.
Further, from the latter, we conclude that
Figure 1. Graph of for .
for every and .
Now, let , . We choose piecewise linear, continuous shape functions of the form
and zero otherwise, for every ,
and zero otherwise. These “hat functions” fulfill
In the following, we want to compute , . Therefore, we prove that , such that we can represent in the Hilbert basis . We compute for ,
and, analogously as
As a consequence,
is summable, and therefore
is summable. The latter implies that
for every and . From (11), (12), we know that for satisfying
for almost all , . Analogously, we obtain
for almost all . A graphical representation of the shape functions and the expressions for is given in Figures 2-6.
We note that for every , the corresponding vanishes outside the interval
Also are linearly independent, since if are such that
for every . Hence, following the approach in Section 0, we arrive at the system of linear equations
Figure 2. Graphs of and for , , .
Figure 3. Graphs of and for , , .
Figure 4. Graphs of and for , , .
Figure 5. Graphs of and for , , .
Figure 6. Graphs of and for , , .
for almost all ,
for every and
is a symmetric and positive definite -matrix. As a consequence,
for every .
4.3. Numerical Implementation
To solve the system of linear Equations (14), we have to approximate the integrals (15) and (16). Thereby, the weak singularities of the functions and have to be taken into account. Accordingly, we split the integrals at the critical points and introduce a Gauss-Jacobi quadrature. In general, a Gauss-Jacobi quadrature is an approximation of an integral over the interval of a continuous function F weighted by an algebraic function with (possibly) weak singularities at the boundaries of the integration interval. It has the form
with , nodes and weights such that polynomials of degree are integrated exactly. The details on determining the nodes and weights may be found e.g. in (  Ch. 2.7). In the following, we derive the quadrature formulas for (15) and (16). As has a weak singularity at , we filter the integrand in (15) by a term , where and obtain the representation
In Figure 7, the smoothing effect of in is shown for two examples.
Hence, a linear transformation of the integrals leads to
Figure 7. Graphs of and for , , , , (left) and , (right).
Finally, we obtain the quadrature formula
In a similar fashion, the integral in (16) can be approximated, where the integrand has a singularity at . As vanishes outside the interval , we obtain
for , where
In Figure 8, the effect of in is shown for two examples.
Again, using a linear transformation results in the quadrature formula
for . Analogously, we obtain
Using the above quadrature formulas, the solution of (14) is approximated for parameters
In Figure 9 and Figure 10, the results in comparison to those from classical integration methods to solve (8) are shown, in particular, a forward difference, backward difference and a classical finite element method as derived in .
Furthermore, the convergence behavior of the schemes mentioned above is studied in Figure 11. Thereby, for the method introduced in this article, a steep decrease of the mean relative error similar as for the classical finite element method may be observed. However, the mean relative error appears to be related to the number N of nodes in the Gaussian quadrature. Hence, in convergence, Gauss points are chosen for the fractional finite element method.
The present paper investigates a finite element method to solve first-order differential equations. Thereby, a fractional power of a differential operator is used to obtain a symmetric system matrix in order to solve the problem with common finite element software. The method is applied to a simple first-order ordinary
Figure 8. Graphs of and for , , , (left) and , (right).
Figure 9. Numerical solutions of (8) using several integration schemes.
Figure 10. Relative error of numerical solutions of (8) using several integration schemes.
Figure 11. Mean relative error of numerical solutions of (8) using several integration schemes depending on the step size/number of elements.
differential equation and the related numerical scheme shows good results compared to classical integration schemes. A generalization to more complex problems is described in Section 2. However, not for all densely-defined, linear and closed operators A, it is possible to explicitly calculate the polar decomposition. For this reason, in future, it will be tried to derive an abstract decomposition for such operators, of the form ; in this connection we note that , is densely-defined, linear and self-adjoint, which avoids the introduction of the intermediate (higher order) operator . This should lead to a simplification of the process in future and pave the way to apply the method to more general, higher dimensional problems.
 Pazy, A. (1983) Semigroups of Linear Operators and Applications to Partial Differential Equations. Applied Mathematical Sciences, Vol. 44, Springer, Berlin.
 Ashyralyev, A. (2009) A Note on Fractional Derivatives and Fractional Powers of Operators. Journal of Mathematical Analysis and Applications, 357, 232-236.