Finite Element Approach for the Solution of First-Order Differential Equations

The finite element method has established itself as an efficient numerical procedure for the solution of arbitrary-shaped field problems in space. Basi-cally, the finite element method transforms the underlying differential equation into a system of algebraic equations by application of the method of weighted residuals in conjunction with a finite element ansatz. However, this procedure is restricted to even-ordered differential equations and leads to symmetric system matrices as a key property of the finite element method. This paper aims in a generalization of the finite element method towards the solution of first-order differential equations. This is achieved by an approach which replaces the first-order derivative by fractional powers of operators making use of the square root of a Sturm-Liouville operator. The resulting procedure incorporates a finite element formulation and leads to a symmetric but dense system matrix. Finally, the scheme is applied to the barometric equation where the results are compared with the analytical solution and other numerical approaches. It turns out that the resulting numerical scheme shows excellent convergence properties.


Introduction
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 How to cite this paper: Schmidt, A., Beyer, H.R., Hinze, M. and Vandoros, E.N. (2020) Finite Element Approach for the Solution of First-Order Differential Equations. Journal of Applied Mathematics and Physics, 8, 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 [1]. For this reason, in the following a different approach is applied which makes use of fractional powers of operators [2] [3] [4]. 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 [5].
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.

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 which are given on the hyperplane { } 0 t × Ω , that is assumed to be non-characteristic. Further, let ( ) ∈ Ω By introducing the new "unknown" function ( ) ∈ Ω and a transformed differential equation Hence, we arrive at an operator equation .
The operator Equation (4) has a unique solution which may be approximated as described in the next section.

Transformation of the Operator Equation into a Suitable form for the Application of Finite Elements
( ) is equivalent to the equation The operator in X, ( )

* 2
A A is densely-defined, linear and positive self-adjoint and bijective. Therefore (6) has the unique solution ( ) (6) can be solved by the usual finite element methods. For this purpose, , where . . denotes the scalar product in X. In the following, we are going to solve the corresponding "approximate" system ( ) . Hence, we arrive at the system of linear equations is symmetric and positive definite. As a consequence,

Operator Equation
for every x I ∈ and for every 0 x > . Hence, we arrive at the transformation of the system of Equations, (8) and (9), into an operator equation of the form (4), where (4) has the unique solution (5). In particular, if 0 q = ,

Derivation of the Finite Element Formulation
Since , is bijective, we need to calculate the corresponding operators First, we note that : , for every x I ∈ , where k ∈  , then where we used the product-to-sum identity (11) and the series representation found in ( [10], p. 1073, Formula (19)). As a consequence, x y x y a a G x y x y x y a a for almost all ( ) 2 , x y I ∈ . The function ( ) 1 2 , G x y is depicted in Figure 1.
Further, from the latter, we conclude that   In the following, we want to compute ( ) , such that we can represent ( )           1, , x y a G x y x y a

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 ( )   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 [1].
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, 5 N n = Gauss points are chosen for the fractional finite element method.

Conclusion
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   . 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 ( )