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

André Schmidt^1*, Horst R. Beyer², Matthias Hinze¹, Evangelos N. Vandoros¹

¹Institute for Nonlinear Mechanics, Stuttgart, Germany.
²Instituto Tecnológico Superior de Uruapan, Uruapan, México.

DOI: 10.4236/jamp.2020.810155                

Abstract

The finite element method has established itself as an efficient numerical procedure for the solution of arbitrary-shaped field problems in space. Basically, 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.

Keywords

Finite Element Method, First-Order Differential Equations, Fractional Powers of Operators

1. 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 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.

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 $a>0$, $I:=\left(0,a\right)$, $n\in \mathbb{N}$, $\Omega \subset {ℝ}^n$, $a_k\in C\left(\Omega ,ℂ\right)$, $k=0,\cdots ,n$, $u\in C^1\left(I\times \Omega \mathrm{,}ℂ\right)$. We consider a linear hyperbolic first-order partial differential equation

${\partial }_{t}u\left(t,x\right)+{\displaystyle \underset{k=1}{\overset{n}{\sum }}}{a}_k\left(x\right){\partial }_{k}u\left(t,x\right)+a_0\left(x\right)u\left(t,x\right)=0$ (1)

for $t\in I$, $x\in \Omega$, ${\partial }_t=\frac{\partial }{\partial t}$, ${\partial }_k=\frac{\partial }{\partial x_k}$ and initial data

$u\left(t_0\mathrm{,}x\right)=u_0\left(x\right)\mathrm{,}x\in \Omega \mathrm{,}{u}_0\mathrm{<mo>\; :\; </mo>}\Omega \to ℂ\mathrm{,}$ (2)

which are given on the hyperplane $\left\{t_0\right\}\times \Omega$, that is assumed to be non-characteristic. Further, let $U\in C^1\left(I\times \Omega \mathrm{,}ℂ\right)$ be such that

$U\left(t_0\mathrm{,}x\right)=u_0\left(x\right)\mathrm{,}x\in \Omega \mathrm{.}$

By introducing the new “unknown” function $f\in C^1\left(I\times \Omega \mathrm{,}ℂ\right)$ as

$f=u-U,$

we obtain homogeneous initial data

$f\left(t_0\mathrm{,}x\right)=\mathrm{0,}x\in \Omega$

and a transformed differential equation

${\partial }_{t}f+{\displaystyle \underset{k=1}{\overset{n}{\sum }}}{a}_k{\partial }_{k}f+a_{0}f=-{\partial }_{t}U-{\displaystyle \underset{k=1}{\overset{n}{\sum }}}{a}_k{\partial }_{k}U-a_{0}U.$ (3)

Hence, we arrive at an operator equation

$Af=g,$ (4)

where

$A:={\partial }_t+{\displaystyle \underset{k=1}{\overset{n}{\sum }}}{a}_k{\partial }_k+a_0$

is a linear operator in $X=L_{ℂ}^2\left(I\times \Omega \right)$ and

$g:=-\left({\partial }_t+{\displaystyle \underset{k=1}{\overset{n}{\sum }}}{a}_k{\partial }_k+a_0\right)U.$

The operator Equation (4) has a unique solution

$f=A^{-1}g,$ (5)

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 [6], which is uniquely associated with A, i.e.

$A=V{\left(A^{\mathrm{<mo>\; *\; </mo>}}A\right)}^{1/2}\mathrm{,}$

where $V\in L\left(X\mathrm{,}X\right)$ is a partial isometry with initial space $\overline{\text{Ran}\left(A^{\mathrm{<mo>\; *\; </mo>}}\right)}$ and end space $\overline{\text{Ran}\left(A\right)}$. Herein, $\text{Ran}(\cdot )$ denotes the range of an operator. Since A is bijective, we have

$\overline{\text{Ran}\left(A^{\mathrm{<mo>\; *\; </mo>}}\right)}=\overline{\text{Ran}\left(A\right)}$

and V is a unitary transformation, i.e.

$V^{\mathrm{<mo>\; *\; </mo>}}V={\text{id}}_X\mathrm{.}$

We note that (4), since,

${\left(A^{*}A\right)}^{1/2}f=V^{*}V{\left(A^{*}A\right)}^{1/2}f=V^{*}Af=V^{*}g,$

is equivalent to the equation

${\left(A^{\mathrm{<mo>\; *\; </mo>}}A\right)}^{1/2}f=V^{\mathrm{<mo>\; *\; </mo>}}g\mathrm{.}$ (6)

The operator in X, ${\left(A^{\mathrm{<mo>\; *\; </mo>}}A\right)}^{1/2}$ is densely-defined, linear and positive self-adjoint and bijective. Therefore (6) has the unique solution

$f={\left[{\left(A^{\mathrm{<mo>\; *\; </mo>}}A\right)}^{1/2}\right]}^{-1}{V}^{\mathrm{<mo>\; *\; </mo>}}g\mathrm{.}$

Now (6) can be solved by the usual finite element methods. For this purpose, let $n\in {\mathbb{N}}^{\mathrm{<mo>\; *\; </mo>}}$, ${\phi }_1\mathrm{,}\cdots \mathrm{,}{\phi }_n$ be linearly independent elements of $D\left({\left(A^{\mathrm{<mo>\; *\; </mo>}}A\right)}^{1/2}\right)$. Further, we denote by $P_n\in L\left(X\mathrm{,}X\right)$ the orthogonal projection onto $\mathcal{L}\left(\left\{{\phi }_1\mathrm{,}\cdots \mathrm{,}{\phi }_n\right\}\right)$. Then

$\begin{array}{c}\langle {\phi }_k|V^{\mathrm{<mo>\; *\; </mo>}}g\rangle =\langle {\phi }_k|{\left(A^{\mathrm{<mo>\; *\; </mo>}}A\right)}^{1/2}f\rangle =\langle {\left(A^{\mathrm{<mo>\; *\; </mo>}}A\right)}^{1/2}{\phi }_k|f\rangle \\ =\langle {\left(A^{\mathrm{<mo>\; *\; </mo>}}A\right)}^{1/2}{\phi }_k|P_{n}f\rangle +\langle {\left(A^{\mathrm{<mo>\; *\; </mo>}}A\right)}^{1/2}{\phi }_k|\left(1-P_n\right)f\rangle \mathrm{,}\end{array}$

for every $k\in \left\{\mathrm{1,}\cdots \mathrm{,}n\right\}$, where $\langle .|.\rangle$ denotes the scalar product in X. In the following, we are going to solve the corresponding “approximate” system

$\langle {\phi }_k|V^{\mathrm{<mo>\; *\; </mo>}}g\rangle =\langle {\left(A^{\mathrm{<mo>\; *\; </mo>}}A\right)}^{1/2}{\phi }_k|P_{n}f\rangle \mathrm{,}$ (7)

where $k\in \left\{\mathrm{1,}\cdots \mathrm{,}n\right\}$. The vector $P_{n}f$ is an element of $\mathcal{L}\left(\left\{{\phi }_1\mathrm{,}\cdots \mathrm{,}{\phi }_n\right\}\right)$. Hence, there are uniquely determined ${\alpha }_1\mathrm{,}\cdots \mathrm{,}{\alpha }_n$ such that

$P_{n}f={\displaystyle \underset{l=1}{\overset{n}{\sum }}}{\alpha }_l{\phi }_l.$

As a consequence,

$\begin{array}{c}\langle {\phi }_k|V^{*}g\rangle =\langle {\left(A^{*}A\right)}^{1/2}{\phi }_k|P_{n}f\rangle ={\displaystyle \underset{l=1}{\overset{n}{\sum }}}{\alpha }_l\langle {\left(A^{*}A\right)}^{1/2}{\phi }_k|{\phi }_l\rangle \\ ={\displaystyle \underset{l=1}{\overset{n}{\sum }}}\langle {\phi }_k|{\left(A^{*}A\right)}^{1/2}{\phi }_l\rangle {\alpha }_l,\end{array}$

for every $k\in \left\{\mathrm{1,}\cdots \mathrm{,}n\right\}$. Hence, we arrive at the system of linear equations

${\displaystyle \underset{l=1}{\overset{n}{\sum }}}{M}_{kl}{\alpha }_l=\langle {\phi }_k|V^{*}g\rangle ,$

$k=1,\cdots ,n$, where the real $n\times n$ matrix

$M:={\left(\langle {\phi }_k|{\left(A^{*}A\right)}^{1/2}{\phi }_l\rangle \right)}_{k,l=1,\cdots ,n}$

is symmetric and positive definite. As a consequence,

${\alpha }_l={\displaystyle \underset{k=1}{\overset{n}{\sum }}}{\left(M^{-1}\right)}_{lk}\langle {\phi }_k|V^{*}g\rangle ,$

for every $l\in \left\{\mathrm{1,}\cdots \mathrm{,}n\right\}$.

4. Application to the Barometric Equation

4.1. Operator Equation

Let $a>0$, $I:=\left(0,a\right)$, $q\in C\left(I\mathrm{,}ℂ\right)$, $p\in C^1\left(I\mathrm{,}ℂ\right)$ be such that

$p^{\prime }\left(x\right)+cp\left(x\right)=q\left(x\right),$ (8)

for every $x\in I$ and

$\underset{x\to 0}{\lim }p\left(x\right)=p_0.$ (9)

Further, let $P\in C^1\left(I\mathrm{,}ℂ\right)$ be such that

$\underset{x\to 0}{\lim }P\left(x\right)=p_0\mathrm{.}$

We define a new “unknown” function $f\in C^1\left(I\mathrm{,}ℂ\right)$ by

$f\left(x\right)={\text{e}}^{cx}\left(p\left(x\right)-P\left(x\right)\right)\mathrm{.}$

Then

$\underset{x\to 0}{\lim }f\left(x\right)=0$

and

$\begin{array}{c}{f}^{\prime }\left(x\right)={\text{e}}^{cx}\left[p^{\prime }\left(x\right)-P^{\prime }\left(x\right)\right]+cf\left(x\right)\\ ={\text{e}}^{cx}\left[p^{\prime }\left(x\right)-P^{\prime }\left(x\right)\right]+c{\text{e}}^{cx}\left[p\left(x\right)-P\left(x\right)\right]\\ ={\text{e}}^{cx}\left[p^{\prime }\left(x\right)+cp\left(x\right)\right]-{\text{e}}^{cx}\left[P^{\prime }\left(x\right)+cP\left(x\right)\right]\\ ={\text{e}}^{cx}q\left(x\right)-{\text{e}}^{cx}\left[P^{\prime }\left(x\right)+cP\left(x\right)\right]\\ ={\text{e}}^{cx}\left[q\left(x\right)-P^{\prime }\left(x\right)-cP\left(x\right)\right],\end{array}$

for every $x>0$. Hence, we arrive at the transformation of the system of Equations, (8) and (9), into an operator equation of the form (4), where

$A:=\overline{D_{\text{I}}},g\left(x\right):={\text{e}}^{cx}\left[q\left(x\right)-P^{\prime }\left(x\right)-cP\left(x\right)\right],$

and $D_{\text{I}}$ denotes the derivative operator with domain

$D\left(D_{\text{I}}\right)=\left\{h\in C^1\left(\bar{I}\mathrm{,}ℂ\right)\mathrm{<mo>\; |\; </mo>}\underset{x\to 0}{\lim }h\left(x\right)=0\right\}$ (10)

in $X=L_{ℂ}^2\left(I\right)$, where (4) has the unique solution (5). In particular, if $q=0$,

$P\left(x\right):=p_0\left(1-\frac{cx}{1+ca}\right),$

for every $x\in I$, then

$\begin{array}{c}\left(P^{\prime }+cP\right)\left(x\right)=-\frac{p_{0}c}{1+ca}+p_{0}c\left(1-\frac{cx}{1+ca}\right)\\ =-\frac{p_{0}c}{1+ca}+p_{0}c\frac{1+ca-cx}{1+ca}\\ =\frac{p_0{c}^2}{1+ca}\left(a-x\right)\end{array}$

and hence

$g\left(x\right)=-\frac{p_0{c}^2}{1+ca}\left(a-x\right){\text{e}}^{cx}\mathrm{,}$

for every $x\in I$, implying that

$g\in D\left(A^{\mathrm{<mo>\; *\; </mo>}}\right)\supset \left\{h\in C^1\left(\bar{I}\mathrm{,}ℂ\right)\mathrm{<mo>\; |\; </mo>}\underset{x\to a}{\lim }h\left(x\right)=0\right\}\mathrm{.}$

4.2. Derivation of the Finite Element Formulation

Since

$A=\overline{D_{\text{I}}}\mathrm{,}$

where $I:=\left(0,a\right)$, is bijective, we need to calculate the corresponding operators

${\left(A^{\mathrm{<mo>\; *\; </mo>}}A\right)}^{1/2}$ and $V^{\mathrm{<mo>\; *\; </mo>}}={\left(A^{\mathrm{<mo>\; *\; </mo>}}A\right)}^{1/2}{A}^{-1}$.

First, we note that

$A^{\mathrm{<mo>\; *\; </mo>}}A={\overline{D_{\text{I}}}}^{\mathrm{<mo>\; *\; </mo>}}\overline{D_{\text{I}}}=\overline{D_{\text{I}\mathrm{,}\text{ad}}}\overline{D_{\text{I}}}\mathrm{.}$

Hence if $f\in D_0$, where

$D_0\mathrm{<mo>\; :\; </mo>}=\left\{f\in C^2\left(\bar{I}\mathrm{,}ℂ\right)\mathrm{<mo>\; :\; </mo>}\underset{x\to 0}{\lim }f\left(x\right)=\mathrm{0,}\underset{x\to a}{\lim }{f}^{\prime }\left(x\right)=0\right\}\subset D\left(A^{\mathrm{<mo>\; *\; </mo>}}A\right)\mathrm{,}$

then

$A^{\mathrm{<mo>\; *\; </mo>}}Af=\overline{D_{\text{I}\mathrm{,}\text{ad}}}{f}^{\prime }=-f^{″}\mathrm{.}$

Hence, $A^{\mathrm{<mo>\; *\; </mo>}}A$ is a densely-defined, linear and self-adjoint extension of the densely-defined, linear, symmetric and essentially self-adjoint operator $A_0$, defined by

$A_0\mathrm{<mo>\; :\; </mo>}{D}_0\to X\mathrm{,}$

and

$A_{0}f\mathrm{<mo>\; :\; </mo>}=-f^{″}\mathrm{,}$

for every $f\in D_0$. Since $A^{\mathrm{<mo>\; *\; </mo>}}A$ is, in particular, closed, it follows that

$A^{\mathrm{<mo>\; *\; </mo>}}A\supset {\bar{A}}_0$

and hence that

$A^{\mathrm{<mo>\; *\; </mo>}}A={\bar{A}}_0\mathrm{.}$

In the following, we are going to calculate ${\left(A^{\mathrm{<mo>\; *\; </mo>}}A\right)}^{1/2}$ and $V^{\mathrm{<mo>\; *\; </mo>}}$, using an approach via the theory of Sturm-Liouville operators [6] [7] [8] [9]. For this purpose, we need the eigenvalues and eigenfunctions of $A^{\mathrm{<mo>\; *\; </mo>}}A$. If $\lambda >0$, the solutions $u\in C^2\left(\bar{I}\mathrm{,}ℂ\right)$ of

$-u^{″}\left(x\right)-\lambda u\left(x\right)=\mathrm{0,}$

for every $x\in I$ are given by

$u\left(x\right)=\alpha sin\left(\sqrt{\lambda }x\right)+\beta cos\left(\sqrt{\lambda }x\right)\mathrm{,}$

for every $x\in I$, where $\alpha \mathrm{,}\beta \in ℂ$. The boundary condition

$\underset{x\to 0}{\lim }u\left(x\right)=0\mathrm{,}$

gives

$u\left(x\right)=\alpha sin\left(\sqrt{\lambda }x\right)\mathrm{,}$

for every $x\in I$, where $\alpha \in ℂ$, and the boundary condition

$\underset{x\to a}{\lim }{u}^{\prime }\left(x\right)=0$

gives a non-trivial u iff

$cos\left(\sqrt{\lambda }a\right)=0\mathrm{,}$

i.e. iff

$\lambda ={\lambda }_k:=\frac{{\pi }^2}{a^2}{\left(k+\frac{1}{2}\right)}^2,$

for some $k\in \mathbb{N}$. For such $k\in \mathbb{N}$, it follows that

$\begin{array}{c}{\displaystyle {\int }_0^a}{\sin }^2\left(\sqrt{\lambda }x\right)\text{d}x={\displaystyle {\int }_0^a}\frac{1}{2}\left[1-\cos \left(2\sqrt{\lambda }x\right)\right]\text{d}x\\ =\frac{a}{2}-\frac{1}{2}{\displaystyle {\int }_0^a}\cos \left(2\sqrt{\lambda }x\right)\text{d}x\\ =\frac{a}{2}-\frac{1}{2}{\left[\frac{\sin \left(2\sqrt{\lambda }x\right)}{2\sqrt{\lambda }}\right]}_0^a=\frac{a}{2}\end{array}$

and hence that, if $e_k\in C^{\infty }\left(I\mathrm{,}ℝ\right)$ is defined by

$e_k\left(x\right)\mathrm{<mo>\; :\; </mo>}=\sqrt{\frac{2}{a}}sin\left(\left(k+\frac{1}{2}\right)\frac{\pi x}{a}\right)\mathrm{,}$

for every $x\in I$, where $k\in \mathbb{N}$, then

$e_0\mathrm{,}{e}_1\mathrm{,}{e}_2\mathrm{,}\cdots$

is a Hilbert basis of X. For

$f\in \left\{g\in C^2\left(\bar{I}\mathrm{,}ℂ\right)\mathrm{<mo>\; :\; </mo>}\underset{x\to 0}{\lim }g\left(x\right)=0\wedge \underset{x\to a}{\lim }{g}^{\prime }\left(a\right)=0\right\}\mathrm{,}$

it follows that

$A^{*}Af=A^{*}A{\displaystyle \underset{k=0}{\overset{\infty }{\sum }}}\langle e_k|f\rangle e_k={\displaystyle \underset{k=0}{\overset{\infty }{\sum }}}\langle e_k|f\rangle A^{*}A{e}_k={\displaystyle \underset{k=0}{\overset{\infty }{\sum }}}{\lambda }_k\langle e_k|f\rangle e_k$

and hence that

$\langle e_k|-f^{″}\rangle =\langle e_k|A^{*}Af\rangle ={\lambda }_k\langle e_k|f\rangle ,\langle e_k|f\rangle =\frac{1}{{\lambda }_k}\langle e_k|-f^{″}\rangle ,$

for every $k\in \mathbb{N}$. Hence, it follows that

$\begin{array}{c}\left(A^{*}{A}^{1/2}\right)f={\left(A^{\mathrm{<mo>\; *\; </mo>}}A\right)}^{1/2}{\displaystyle \underset{k=0}{\overset{\infty }{\sum }}}\langle e_k|f\rangle e_k={\displaystyle \underset{k=0}{\overset{\infty }{\sum }}}\langle e_k|f\rangle {\left(A^{\mathrm{<mo>\; *\; </mo>}}A\right)}^{1/2}{e}_k\\ ={\displaystyle \underset{k=0}{\overset{\infty }{\sum }}}\langle e_k|f\rangle \sqrt{{\lambda }_k}{e}_k={\displaystyle \underset{k=0}{\overset{\infty }{\sum }}}\frac{1}{{\lambda }_k}\langle e_k|-f^{″}\rangle \sqrt{{\lambda }_k}{e}_k\\ ={\displaystyle \underset{k=0}{\overset{\infty }{\sum }}}\frac{1}{\sqrt{{\lambda }_k}}\langle e_k|-f^{″}\rangle e_k\end{array}$

and for almost all $x\in I$,

$\begin{array}{c}\left[\left(A^{*}{A}^{1/2}\right)f\right]\left(x\right)={\displaystyle \underset{k=0}{\overset{\infty }{\sum }}}\frac{1}{\sqrt{{\lambda }_k}}\langle e_k|-f^{″}\rangle e_k\left(x\right)\\ =\langle {\displaystyle \underset{k=0}{\overset{\infty }{\sum }}}\frac{1}{\sqrt{{\lambda }_k}}{e}_k\left(x\right)e_k|-f^{″}\rangle \\ ={\displaystyle {\int }_0^a}\left[{\displaystyle \underset{k=0}{\overset{\infty }{\sum }}}\frac{1}{\sqrt{{\lambda }_k}}{e}_k\left(x\right)e_k\left(y\right)\right]\left(-f^{″}\left(y\right)\right)\text{d}y,\end{array}$

where we used that

${\left({\displaystyle \underset{k=0}{\overset{n}{\sum }}}\frac{1}{\sqrt{{\lambda }_k}}{e}_k\left(x\right)e_k\right)}_{n\in \mathbb{N}}$

is convergent in X, as a consequence of the estimations

${\left|\frac{1}{\sqrt{{\lambda }_k}}{e}_k\left(x\right)\right|}^2=\frac{1}{{\lambda }_k}{\left|e_k\left(x\right)\right|}^2\le \frac{2}{a}\cdot \frac{a^2}{{\pi }^2}\frac{1}{{\left(k+\frac{1}{2}\right)}^2}=\frac{8a}{{\pi }^2}\frac{1}{{\left(2k+1\right)}^2}\mathrm{,}$

$k\in {\mathbb{N}}^{\mathrm{<mo>\; *\; </mo>}}$, and as consequence of the existence of

${\displaystyle \underset{k=0}{\overset{\infty }{\sum }}}\frac{1}{{\left(2k+1\right)}^2}.$

We compute for $x\mathrm{,}y\in I$ satisfying $x\ne y$

$\begin{array}{l}{\displaystyle \underset{k=0}{\overset{\infty }{\sum }}}\frac{1}{\sqrt{{\lambda }_k}}{e}_k\left(x\right)e_k\left(y\right)\\ =\frac{2}{\pi }{\displaystyle \underset{k=0}{\overset{\infty }{\sum }}}\frac{1}{k+\frac{1}{2}}\sin \left(\left(k+\frac{1}{2}\right)\frac{\pi x}{a}\right)\sin \left(\left(k+\frac{1}{2}\right)\frac{\pi y}{a}\right)\\ =\frac{1}{\pi }{\displaystyle \underset{k=0}{\overset{\infty }{\sum }}}\frac{1}{k+\frac{1}{2}}\left[\cos \left(\left(k+\frac{1}{2}\right)\frac{\pi \left(x-y\right)}{a}\right)-\cos \left(\left(k+\frac{1}{2}\right)\frac{\pi \left(x+y\right)}{a}\right)\right]\\ =\frac{2}{\pi }{\displaystyle \underset{k=0}{\overset{\infty }{\sum }}}\frac{1}{2k+1}\left[\cos \left(\left(2k+1\right)\frac{\pi \left(x-y\right)}{2a}\right)-\cos \left(\left(2k+1\right)\frac{\pi \left(x+y\right)}{2a}\right)\right]\end{array}$

$\begin{array}{l}=\frac{2}{\pi }{\displaystyle \underset{k=0}{\overset{\infty }{\sum }}}\frac{1}{2k+1}\left[\cos \left(\left(2k+1\right)\frac{\pi \left|x-y\right|}{2a}\right)-\cos \left(\left(2k+1\right)\frac{\pi \left(x+y\right)}{2a}\right)\right]\\ =\frac{1}{\pi }\left\{\ln \left[\cot \left(\frac{\pi \left|x-y\right|}{4a}\right)\right]-\ln \left[\cot \left(\frac{\pi \left(x+y\right)}{4a}\right)\right]\right\}\\ =\frac{1}{\pi }\ln \left[\frac{\cot \left(\frac{\pi \left|x-y\right|}{4a}\right)}{\cot \left(\frac{\pi \left(x+y\right)}{4a}\right)}\right],\end{array}$

where we used the product-to-sum identity

$2sin\left(\alpha \right)sin\left(\beta \right)=cos\left(\alpha -\beta \right)-cos\left(\alpha +\beta \right)\mathrm{,}\alpha \mathrm{,}\beta \in ℝ$ (11)

and the series representation

$\frac{1}{2}\ln \left(\cot \left(\frac{x}{2}\right)\right)={\displaystyle \underset{k=0}{\overset{\infty }{\sum }}}\frac{1}{2k+1}\cos \left(\left(2k+1\right)x\right)$ (12)

found in ( [10], p. 1073, Formula (19)). As a consequence,

$\left[{\left(A^{\mathrm{<mo>\; *\; </mo>}}A\right)}^{1/2}f\right]\left(x\right)=-{\displaystyle {\int }_0^a}{G}_{\frac{1}{2}}\left(x\mathrm{,}y\right)f^{″}\left(y\right)\text{d}y\mathrm{,}$

for every $x\in I$, where

$G_{\frac{1}{2}}\left(x\mathrm{,}y\right)\mathrm{<mo>\; :\; </mo>}=\frac{1}{\pi }ln\left[\frac{cot\left(\frac{\pi \left|x-y\right|}{4a}\right)}{cot\left(\frac{\pi \left(x+y\right)}{4a}\right)}\right]=\frac{1}{\pi }ln\left[\frac{tan\left(\frac{\pi \left(x+y\right)}{4a}\right)}{tan\left(\frac{\pi \left|x-y\right|}{4a}\right)}\right]\mathrm{,}$

for almost all $\left(x\mathrm{,}y\right)\in I^2$. The function $G_{\frac{1}{2}}\left(x\mathrm{,}y\right)$ is depicted in Figure 1.

Further, from the latter, we conclude that

Figure 1. Graph of $G_{\frac{1}{2}}$ for $a=1$.

$\left(V^{\mathrm{<mo>\; *\; </mo>}}\text{}g\right)\left(x\right)=-{\displaystyle {\int }_0^a}{G}_{\frac{1}{2}}\left(x\mathrm{,}y\right)g^{\prime }\left(y\right)\text{d}y\mathrm{,}$

for every $g\in D\left(A^{\mathrm{<mo>\; *\; </mo>}}\right)$ and $x\in I$.

Now, let $n\in \mathbb{N}\backslash \left\{\mathrm{0,1}\right\}$, $l:=\frac{a}{n}$. We choose piecewise linear, continuous shape functions ${\phi }_1\mathrm{,}\cdots \mathrm{,}{\phi }_n$ of the form

${\phi }_j\left(x\right)=(\begin{array}{ll}\frac{x-\left(j-1\right)l}{l}\hfill & \text{for}x\in \left[\left(j-1\right)l,jl\right]\hfill \\ \frac{\left(j+1\right)l-x}{l}\hfill & \text{for}x\in \left[jl,\left(j+1\right)l\right]\hfill \end{array}$

and zero otherwise, for every $j\in \left\{\mathrm{1,}\cdots \mathrm{,}n-1\right\}$,

${\phi }_n\left(x\right)=\frac{x-\left(n-1\right)l}{l}\text{for}x\in \left[\left(n-1\right)l,nl\right],$

and zero otherwise. These “hat functions” ${\phi }_1\mathrm{,}\cdots \mathrm{,}{\phi }_n$ fulfill

${\phi }_j\left(kl\right)={\delta }_{jk}\mathrm{.}$ (13)

In the following, we want to compute ${\left(A^{\mathrm{<mo>\; *\; </mo>}}A\right)}^{1/2}{\phi }_j$, $j=\mathrm{1,}\cdots \mathrm{,}n$. Therefore, we prove that ${\phi }_j\in D\left({\left(A^{\mathrm{<mo>\; *\; </mo>}}A\right)}^{1/2}\right)$, such that we can represent ${\left(A^{\mathrm{<mo>\; *\; </mo>}}A\right)}^{1/2}{\phi }_j$ in the Hilbert basis ${\left\{e_k\right\}}_{k=0,1,\cdots }$. We compute for $k\in \mathbb{N}$, $j\in \left\{\mathrm{1,}\cdots \mathrm{,}n-1\right\}$

$\begin{array}{c}\langle e_k|{\phi }_j\rangle =\sqrt{\frac{2}{a}}{\displaystyle {\int }_0^a}\sin \left(\sqrt{{\lambda }_k}x\right){\phi }_j\left(x\right)\text{d}x\\ =\sqrt{\frac{2}{a}}{\displaystyle {\int }_{\left(j-1\right)l}^{jl}}\sin \left(\sqrt{{\lambda }_k}x\right)\frac{x-\left(j-1\right)l}{l}\text{d}x\\ +\sqrt{\frac{2}{a}}{\displaystyle {\int }_{jl}^{\left(j+1\right)l}}\sin \left(\sqrt{{\lambda }_k}x\right)\frac{\left(j+1\right)l-x}{l}\text{d}x\\ =\frac{1}{l}\sqrt{\frac{2}{a}}{\left[\frac{\sin \left(\sqrt{{\lambda }_k}x\right)}{{\lambda }_k}-\left(x-\left(j-1\right)l\right)\frac{\cos \left(\sqrt{{\lambda }_k}x\right)}{\sqrt{{\lambda }_k}}\right]}_{\left(j-1\right)l}^{jl}\\ +\frac{1}{l}\sqrt{\frac{2}{a}}{\left[-\frac{\sin \left(\sqrt{{\lambda }_k}x\right)}{{\lambda }_k}-\left(\left(j+1\right)l-x\right)\frac{\cos \left(\sqrt{{\lambda }_k}x\right)}{\sqrt{{\lambda }_k}}\right]}_{jl}^{\left(j+1\right)l}\end{array}$

$\begin{array}{l}=\frac{1}{l}\sqrt{\frac{2}{a}}\left(\frac{\sin \left(\sqrt{{\lambda }_k}jl\right)}{{\lambda }_k}-\frac{\sin \left(\sqrt{{\lambda }_k}\left(j-1\right)l\right)}{{\lambda }_k}-l\frac{\cos \left(\sqrt{{\lambda }_k}jl\right)}{\sqrt{{\lambda }_k}}\right)\\ -\frac{1}{l}\sqrt{\frac{2}{a}}\left(\frac{\sin \left(\sqrt{{\lambda }_k}\left(j+1\right)l\right)}{{\lambda }_k}-\frac{\sin \left(\sqrt{{\lambda }_k}jl\right)}{{\lambda }_k}-l\frac{\cos \left(\sqrt{{\lambda }_k}jl\right)}{\sqrt{{\lambda }_k}}\right)\\ =\frac{1}{l}\sqrt{\frac{2}{a}}\frac{1}{{\lambda }_k}\left(\sin \left(\sqrt{{\lambda }_k}jl\right)-\sin \left(\sqrt{{\lambda }_k}\left(j-1\right)l\right)\right)\\ -\frac{1}{l}\sqrt{\frac{2}{a}}\frac{1}{{\lambda }_k}\left(\sin \left(\sqrt{{\lambda }_k}\left(j+1\right)l\right)-\sin \left(\sqrt{{\lambda }_k}jl\right)\right)\end{array}$

and, analogously as $\sqrt{{\lambda }_k}nl=\frac{\pi }{a}\left(k+\frac{1}{2}\right)a=\pi \left(k+\frac{1}{2}\right)$

$\begin{array}{c}\langle e_k|{\phi }_n\rangle =\frac{1}{l}\sqrt{\frac{2}{a}}\left(\frac{sin\left(\sqrt{{\lambda }_k}nl\right)}{{\lambda }_k}-\frac{sin\left(\sqrt{{\lambda }_k}\left(n-1\right)l\right)}{{\lambda }_k}-l\frac{cos\left(\sqrt{{\lambda }_k}nl\right)}{\sqrt{{\lambda }_k}}\right)\\ =\frac{1}{l}\sqrt{\frac{2}{a}}\frac{1}{{\lambda }_k}\left(sin\left(\sqrt{{\lambda }_k}nl\right)-sin\left(\sqrt{{\lambda }_k}\left(n-1\right)l\right)\right)\mathrm{.}\end{array}$

As a consequence,

${\left({\left|\sqrt{{\lambda }_k}\langle e_k|{\phi }_j\rangle \right|}^2\right)}_{k\in \mathbb{N}}\mathrm{,}j=\mathrm{1,}\cdots \mathrm{,}n$

is summable, and therefore

${\left(\sqrt{{\lambda }_k}\langle e_k|{\phi }_j\rangle e_k\right)}_{k\in \mathbb{N}},j=1,\cdots ,n$

is summable. The latter implies that

${\phi }_j\in D\left({\left(A^{*}A\right)}^{1/2}\right),j=1,\cdots ,n$

and that

${\left(A^{*}A\right)}^{1/2}{\phi }_j={\displaystyle \underset{k=0}{\overset{\infty }{\sum }}}\sqrt{{\lambda }_k}\langle e_k|{\phi }_j\rangle e_k,j=1,\cdots ,n.$

Further, for $y\in I$

$\begin{array}{l}\sqrt{{\lambda }_k}\langle e_k|{\phi }_j\rangle e_k\left(y\right)\\ =\frac{1}{l}\frac{2}{a}\frac{1}{{\lambda }_k}\left(\sin \left(\sqrt{{\lambda }_k}jl\right)-\sin \left(\sqrt{{\lambda }_k}\left(j-1\right)l\right)\right)\sin \left(\sqrt{{\lambda }_k}y\right)\\ -\frac{1}{l}\frac{2}{a}\frac{1}{{\lambda }_k}\left(\sin \left(\sqrt{{\lambda }_k}\left(j+1\right)l\right)-\sin \left(\sqrt{{\lambda }_k}jl\right)\right)\sin \left(\sqrt{{\lambda }_k}y\right)\\ =\frac{4}{\pi l}\frac{1}{2k+1}\left(\sin \left(\left(2k+1\right)\frac{jl\pi }{2a}\right)-\sin \left(\left(2k+1\right)\frac{\left(j-1\right)l\pi }{2a}\right)\right)\\ \times \sin \left(\left(2k+1\right)\frac{\pi y}{2a}\right)\\ -\frac{4}{\pi l}\frac{1}{2k+1}\left(\sin \left(\left(2k+1\right)\frac{\left(j+1\right)l\pi }{2a}\right)-\sin \left(\left(2k+1\right)\frac{jl\pi }{2a}\right)\right)\\ \times \sin \left(\left(2k+1\right)\frac{\pi y}{2a}\right),\end{array}$

for every $k\in \mathbb{N}$ and $j\in \left\{\mathrm{1,}\cdots \mathrm{,}n\right\}$. From (11), (12), we know that for $x\mathrm{,}y\in I$ satisfying $x\ne y$

${\displaystyle \underset{k=0}{\overset{\infty }{\sum }}}\frac{1}{2k+1}\sin \left(\left(2k+1\right)\frac{\pi x}{2a}\right)\sin \left(\left(2k+1\right)\frac{\pi y}{2a}\right)=\frac{1}{4}\ln \left[\frac{\cot \left(\frac{\pi \left|x-y\right|}{4a}\right)}{\cot \left(\frac{\pi \left(x+y\right)}{4a}\right)}\right].$

Hence,

$\begin{array}{l}\left[{\left(A^{\mathrm{<mo>\; *\; </mo>}}A\right)}^{1/2}{\phi }_j\right]\left(y\right)={\displaystyle \underset{k=0}{\overset{\infty }{\sum }}}\sqrt{{\lambda }_k}\langle e_k|{\phi }_j\rangle e_k\left(y\right)\\ =\frac{1}{\pi l}\left\{\ln \left[\frac{\cot \left(\frac{\pi \left|jl-y\right|}{4a}\right)}{\cot \left(\frac{\pi \left[jl+y\right]}{4a}\right)}\right]-\ln \left[\frac{\cot \left(\frac{\pi \left|\left(j-1\right)l-y\right|}{4a}\right)}{\cot \left(\frac{\pi \left[\left(j-1\right)l+y\right]}{4a}\right)}\right]\right\}\\ -\frac{1}{\pi l}\left\{\ln \left[\frac{\cot \left(\frac{\pi \left|\left(j+1\right)l-y\right|}{4a}\right)}{\cot \left(\frac{\pi \left[\left(j+1\right)l+y\right]}{4a}\right)}\right]-\ln \left[\frac{\cot \left(\frac{\pi \left|jl-y\right|}{4a}\right)}{\cot \left(\frac{\pi \left[jl+y\right]}{4a}\right)}\right]\right\}\\ =\frac{1}{\pi l}\{\ln \left[\frac{\tan \left(\frac{\pi \left[jl+y\right]}{4a}\right)\tan \left(\frac{\pi \left|\left(j-1\right)l-y\right|}{4a}\right)}{\tan \left(\frac{\pi \left|jl-y\right|}{4a}\right)\tan \left(\frac{\pi \left[\left(j-1\right)l+y\right]}{4a}\right)}\right]\\ -\ln \left[\frac{\tan \left(\frac{\pi \left[\left(j+1\right)l+y\right]}{4a}\right)\tan \left(\frac{\pi \left|jl-y\right|}{4a}\right)}{\tan \left(\frac{\pi \left|\left(j+1\right)l-y\right|}{4a}\right)\tan \left(\frac{\pi \left[jl+y\right]}{4a}\right)}\right]\},\end{array}$

for almost all $y\in I$, $j\in \left\{\mathrm{1,}\cdots \mathrm{,}n\right\}$. Analogously, we obtain

$\left[{\left(A^{*}A\right)}^{1/2}{\phi }_n\right]\left(y\right)=\frac{1}{\pi l}\ln \left[\frac{\tan \left(\frac{\pi \left[nl+y\right]}{4a}\right)\tan \left(\frac{\pi \left|\left(n-1\right)l-y\right|}{4a}\right)}{\tan \left(\frac{\pi \left|nl-y\right|}{4a}\right)\tan \left(\frac{\pi \left[\left(n-1\right)l+y\right]}{4a}\right)}\right],$

for almost all $y\in I$. A graphical representation of the shape functions ${\phi }_j$ and the expressions ${\left(A^{*}A\right)}^{1/2}{\phi }_j$ for $j=1,\cdots ,5$ is given in Figures 2-6.

We note that for every $j\in \left\{\mathrm{1,}\cdots \mathrm{,}n\right\}$, the corresponding ${\phi }_j$ vanishes outside the interval

$\left(\left(j-1\right)l\mathrm{,}\left(j+1\right)l\right)\mathrm{.}$

Also ${\phi }_1\mathrm{,}\cdots \mathrm{,}{\phi }_n$ are linearly independent, since if ${\alpha }_1\mathrm{,}\cdots \mathrm{,}{\alpha }_n\in ℝ$ are such that

${\displaystyle \underset{k=1}{\overset{n}{\sum }}}{\alpha }_k{\phi }_k=0,$

then

$0={\displaystyle \underset{k=1}{\overset{n}{\sum }}}{\alpha }_k{\phi }_k\left(jl\right)={\alpha }_j,$

for every $j\in \left\{\mathrm{1,}\cdots \mathrm{,}n\right\}$. Hence, following the approach in Section 0, we arrive at the system of linear equations

${\displaystyle \underset{m=1}{\overset{n}{\sum }}}{M}_{km}{\alpha }_m={\displaystyle {\int }_0^a}{\phi }_k\left(x\right)V^{*}g\left(x\right)\text{d}x,$ (14)

Figure 2. Graphs of ${\phi }_1$ and ${\left(A^{*}A\right)}^{1/2}{\phi }_1$ for $a=1$, $n=5$, $l=1/5$.

Figure 3. Graphs of ${\phi }_2$ and ${\left(A^{*}A\right)}^{1/2}{\phi }_2$ for $a=1$, $n=5$, $l=1/5$.

Figure 4. Graphs of ${\phi }_3$ and ${\left(A^{*}A\right)}^{1/2}{\phi }_3$ for $a=1$, $n=5$, $l=1/5$.

Figure 5. Graphs of ${\phi }_4$ and ${\left(A^{*}A\right)}^{1/2}{\phi }_4$ for $a=1$, $n=5$, $l=1/5$.

Figure 6. Graphs of ${\phi }_5$ and ${\left(A^{*}A\right)}^{1/2}{\phi }_5$ for $a=1$, $n=5$, $l=1/5$.

$k=\mathrm{1,}\cdots \mathrm{,}n$, where

${\displaystyle {\int }_0^a}{\phi }_k\left(x\right)V^{\mathrm{<mo>\; *\; </mo>}}g\left(x\right)\text{d}x=-{\displaystyle {\int }_0^a}{\phi }_k\left(x\right){\displaystyle {\int }_0^a}{G}_{\frac{1}{2}}\left(x\mathrm{,}y\right)g^{\prime }\left(y\right)\text{d}y\text{d}x\mathrm{,}$ (15)

with

$G_{\frac{1}{2}}\left(x,y\right)=\frac{1}{\pi }\ln \left[\frac{\tan \left(\frac{\pi \left(x+y\right)}{4a}\right)}{\tan \left(\frac{\pi \left|x-y\right|}{4a}\right)}\right],$

for almost all $\left(x\mathrm{,}y\right)\in I^2$,

$g\left(x\right)=-\frac{p_0{c}^2}{1+ca}\left(a-x\right){\text{e}}^{cx},$

for every $x\in I$ and

$M\mathrm{<mo>\; :\; </mo>}={\left({\displaystyle {\int }_0^a}{\phi }_k\left(x\right){\left(A^{*}A\right)}^{1/2}{\phi }_m\left(x\right)\text{d}x\right)}_{k\mathrm{,}m=\mathrm{1,}\cdots \mathrm{,}n}$ (16)

is a symmetric and positive definite $n\times n$ -matrix. As a consequence,

${\alpha }_j={\displaystyle \underset{k=1}{\overset{n}{\sum }}}{\left(M^{-1}\right)}_{jk}{\displaystyle {\int }_0^a}{\phi }_k\left(x\right)V^{\mathrm{<mo>\; *\; </mo>}}g\left(x\right)\text{d}x\mathrm{,}$ (17)

for every $j\in \left\{\mathrm{1,}\cdots \mathrm{,}n\right\}$.

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 ${\left(A^{*}A\right)}^{1/2}{\phi }_j$ and $G_{\frac{1}{2}}$ 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 $\left[-\mathrm{1,1}\right]$ 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

${\displaystyle {\int }_{-1}^1}F\left(s\right){\left(1-s\right)}^{\beta }{\left(1+s\right)}^{\gamma }\text{d}s\approx {\displaystyle \underset{n=0}{\overset{N}{\sum }}}F\left(s_n^{\left(\beta ,\gamma \right)}\right)w_n^{\left(\beta ,\gamma \right)}$ (18)

with $\beta ,\gamma >-1$, nodes $s_n^{\left(\beta ,\gamma \right)}$ and weights $w_n^{\left(\beta ,\gamma \right)}$ such that polynomials of degree $2N-1$ are integrated exactly. The details on determining the nodes and weights may be found e.g. in ( [11] Ch. 2.7). In the following, we derive the quadrature formulas for (15) and (16). As $G_{\frac{1}{2}}\left(x\mathrm{,}y\right)$ has a weak singularity at $x=y$, we filter the integrand in (15) by a term ${\left|x-y\right|}^{-\alpha }$, where $\alpha \in \left(-\mathrm{1,0}\right)$ and obtain the representation

$-{\displaystyle {\int }_0^a}{\phi }_k\left(x\right){\displaystyle {\int }_0^a}{G}_{\frac{1}{2}}\left(x\mathrm{,}y\right)g^{\prime }\left(y\right)\text{d}y\text{d}x={\displaystyle {\int }_0^a}{\displaystyle {\int }_0^a}{\tilde{G}}_k\left(x\mathrm{,}y\right){\left|x-y\right|}^{\alpha }\text{d}y\text{d}x\mathrm{,}$

where

${\tilde{G}}_k\left(x,y\right):=-{\phi }_k\left(x\right)G_{\frac{1}{2}}\left(x,y\right)g^{\prime }\left(y\right){\left|x-y\right|}^{-\alpha }.$ (19)

In Figure 7, the smoothing effect of ${\left|x-y\right|}^{-\alpha }$ in ${\tilde{G}}_k$ is shown for two examples.

Hence, a linear transformation of the integrals leads to

$\begin{array}{l}{\displaystyle {\int }_0^a}{\displaystyle {\int }_0^a}{\tilde{G}}_k\left(x\mathrm{,}y\right){\left|x-y\right|}^{\alpha }\text{d}y\text{d}x\\ ={\displaystyle {\int }_0^a}\left({\displaystyle {\int }_0^x}{\tilde{G}}_k\left(x,y\right){\left(x-y\right)}^{\alpha }\text{d}y+{\displaystyle {\int }_x^a}{\tilde{G}}_k\left(x,y\right){\left(y-x\right)}^{\alpha }\text{d}y\right)\text{d}x\\ ={\displaystyle {\int }_0^a}\frac{x}{2}{\displaystyle {\int }_{-1}^1}{\tilde{G}}_k\left(x,\frac{x}{2}\left(1+s\right)\right){\left(\frac{x}{2}\left(1-s\right)\right)}^{\alpha }\text{d}s\text{d}x\\ +{\displaystyle {\int }_0^a}\frac{a-x}{2}{\displaystyle {\int }_{-1}^1}{\tilde{G}}_k\left(x,\frac{a-x}{2}s+\frac{a+x}{2}\right){\left(\frac{a-x}{2}\left(1+s\right)\right)}^{\alpha }\text{d}s\text{d}x\end{array}$

$\begin{array}{l}={\displaystyle {\int }_{-1}^1}{\displaystyle {\int }_0^a}{\left(\frac{x}{2}\right)}^{1+\alpha }{\tilde{G}}_k\left(x,\frac{x}{2}\left(1-s\right)\right)\text{d}x{\left(1-s\right)}^{\alpha }\text{d}s\\ +{\displaystyle {\int }_{-1}^1}{\displaystyle {\int }_0^a}{\left(\frac{a-x}{2}\right)}^{1+\alpha }{\tilde{G}}_k\left(x,\frac{a-x}{2}s+\frac{a+x}{2}\right)\text{d}x{\left(1+s\right)}^{\alpha }\text{d}s\\ ={\displaystyle {\int }_{-1}^1}\frac{a}{2}{\displaystyle {\int }_{-1}^1}{\left(\frac{a}{4}\left(1+t\right)\right)}^{1+\alpha }{\tilde{G}}_k\left(\frac{a}{2}\left(1+t\right),\frac{a}{4}\left(1+t\right)\left(1+s\right)\right)\text{d}t{\left(1-s\right)}^{\alpha }\text{d}s\\ +{\displaystyle {\int }_{-1}^1}\frac{a}{2}{\displaystyle {\int }_{-1}^1}{\left(\frac{a}{4}\left(1-t\right)\right)}^{1+\alpha }{\tilde{G}}_k\left(\frac{a}{2}\left(1+t\right),\frac{a}{4}\left(1-t\right)s+\frac{a}{4}\left(3+t\right)\right)\text{d}t{\left(1+s\right)}^{\alpha }\text{d}s\\ ={\displaystyle {\int }_{-1}^1}{\displaystyle {\int }_{-1}^1}{\tilde{G}}_{k,1}\left(s,t\right){\left(1+t\right)}^{1+\alpha }\text{d}t{\left(1-s\right)}^{\alpha }\text{d}s+{\displaystyle {\int }_{-1}^1}{\displaystyle {\int }_{-1}^1}{\tilde{G}}_{k,2}\left(s,t\right){\left(1-t\right)}^{1+\alpha }\text{d}t{\left(1+s\right)}^{\alpha }\text{d}s,\end{array}$

Figure 7. Graphs of ${\tilde{G}}_k\left(x\mathrm{,}y\right){\left|x-y\right|}^{\alpha }$ and ${\tilde{G}}_k\left(x\mathrm{,}y\right)$ for $a=1$, $n=5$, $l=1/5$, $k=1$, $y=l$ (left) and $k=3$, $y=a$ (right).

where

${\tilde{G}}_{k,1}\left(s,t\right)=2{\left(\frac{a}{4}\right)}^{2+\alpha }{\tilde{G}}_k\left(\frac{a}{2}\left(1+t\right),\frac{a}{4}\left(1+t\right)\left(1+s\right)\right),$

${\tilde{G}}_{k,2}\left(s,t\right)=2{\left(\frac{a}{4}\right)}^{2+\alpha }{\tilde{G}}_k\left(\frac{a}{2}\left(1+t\right),\frac{a}{4}\left(1-t\right)s+\frac{a}{4}\left(3+t\right)\right).$

Finally, we obtain the quadrature formula

$\begin{array}{l}{\displaystyle {\int }_0^a}{\displaystyle {\int }_0^a}{\tilde{G}}_k\left(x\mathrm{,}y\right){\left|x-y\right|}^{\alpha }\text{d}y\text{d}x\\ \approx {\displaystyle \underset{i=0}{\overset{N}{\sum }}}{\displaystyle \underset{j=0}{\overset{N}{\sum }}}{\tilde{G}}_{k,1}\left(s_i^{\left(\alpha ,0\right)},t_j^{\left(0,1+\alpha \right)}\right)w_i^{\left(\alpha ,0\right)}{w}_j^{\left(0,1+\alpha \right)}\\ +{\displaystyle \underset{i=0}{\overset{N}{\sum }}}{\displaystyle \underset{j=0}{\overset{N}{\sum }}}{\tilde{G}}_{k,2}\left(s_i^{\left(0,\alpha \right)},t_j^{\left(1+\alpha ,0\right)}\right)w_i^{\left(0,\alpha \right)}{w}_j^{\left(1+\alpha ,0\right)}.\end{array}$ (20)

In a similar fashion, the integral in (16) can be approximated, where the integrand has a singularity at $x=kl$. As ${\phi }_k$ vanishes outside the interval $\left\{\left(k-1\right)l\mathrm{,}\left(k+1\right)l\right\}$, we obtain

$\begin{array}{c}{M}_{km}={\displaystyle {\int }_0^a}{\phi }_k\left(x\right){\left(A^{*}A\right)}^{1/2}{\phi }_m\left(x\right)\text{d}x\\ ={\displaystyle {\int }_{\left(k-1\right)l}^{kl}}{\phi }_k\left(x\right){\left(A^{*}A\right)}^{1/2}{\phi }_m\left(x\right)\text{d}x+{\displaystyle {\int }_{kl}^{\left(k+1\right)l}}{\phi }_k\left(x\right){\left(A^{*}A\right)}^{1/2}{\phi }_m\left(x\right)\text{d}x\\ ={\displaystyle {\int }_{\left(k-1\right)l}^{kl}}{\Phi }_{km}\left(x\right){\left(kl-x\right)}^{\alpha }\text{d}x+{\displaystyle {\int }_{kl}^{\left(k+1\right)l}}{\Phi }_{km}\left(x\right){\left(x-kl\right)}^{\alpha }\text{d}x,\end{array}$

for $k\in \left\{\mathrm{1,}\cdots ,n-1\right\}$, where

${\Phi }_{km}\left(x\right)={\phi }_k\left(x\right){\left(A^{*}A\right)}^{1/2}{\phi }_m\left(x\right){\left|x-kl\right|}^{-\alpha }\mathrm{.}$ (21)

In Figure 8, the effect of ${\left|x-kl\right|}^{\alpha }$ in ${\Phi }_{km}$ is shown for two examples.

Again, using a linear transformation results in the quadrature formula

$\begin{array}{c}{M}_{km}={\left(\frac{l}{2}\right)}^{1+\alpha }{\displaystyle {\int }_{-1}^1}{\Phi }_{km}\left(\frac{l}{2}s+\left(k-\frac{1}{2}\right)l\right){\left(1-s\right)}^{\alpha }\text{d}s\\ +{\left(\frac{l}{2}\right)}^{1+\alpha }{\displaystyle {\int }_{-1}^1}{\Phi }_{km}\left(\frac{l}{2}s+\left(k+\frac{1}{2}\right)l\right){\left(1+s\right)}^{\alpha }\text{d}s\\ \approx {\left(\frac{l}{2}\right)}^{1+\alpha }{\displaystyle \underset{i=0}{\overset{N}{\sum }}}{\Phi }_{km}\left(\frac{l}{2}{s}_i^{\left(\alpha ,0\right)}+\left(k-\frac{1}{2}\right)l\right)w_i^{\left(\alpha ,0\right)}\\ +{\left(\frac{l}{2}\right)}^{1+\alpha }{\displaystyle \underset{i=0}{\overset{N}{\sum }}}{\Phi }_{km}\left(\frac{l}{2}{s}_i^{\left(0,\alpha \right)}+\left(k+\frac{1}{2}\right)l\right)w_i^{\left(0,\alpha \right)}\end{array}$ (22)

for $k\in \left\{\mathrm{1,}\cdots ,n-1\right\}$. Analogously, we obtain

$M_{nm}\approx {\left(\frac{l}{2}\right)}^{1+\alpha }{\displaystyle \underset{i=0}{\overset{N}{\sum }}}{\Phi }_{nm}\left(\frac{l}{2}{s}_i^{\left(\alpha ,0\right)}+\left(n-\frac{1}{2}\right)l\right)w_i^{\left(\alpha ,0\right)}\mathrm{.}$ (23)

Using the above quadrature formulas, the solution of (14) is approximated for parameters

$a=10000\text{m}\mathrm{,}{p}_0=1.013\text{bar}\mathrm{,}c=1.865\times {10}^4{\text{m}}^{-1}$

$n=\mathrm{20,}N=\mathrm{50,}\alpha =-0.5\mathrm{.}$

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

5. 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

Figure 8. Graphs of ${\Phi }_{km}\left(x\right){\left|x-kl\right|}^{\alpha }$ and ${\Phi }_{km}\left(x\right)$ for $a=1$, $n=5$, $l=1/5$, $k=m=1$ (left) and $k=3$, $m=5$ (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 $A=V{\left(A^{1/2}\right)}^{\mathrm{<mo>\; *\; </mo>}}{A}^{1/2}$ ; in this connection we note that ${\left(A^{1/2}\right)}^{\mathrm{<mo>\; *\; </mo>}}{A}^{1/2}$, is densely-defined, linear and self-adjoint, which avoids the introduction of the intermediate (higher order) operator $A^{\mathrm{<mo>\; *\; </mo>}}A$. 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.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Schmidt, A. (2017) Derivation of the Weak Form of a First-Order Differential Equation by Application of Fractional Calculus. Proceedings of ICEDyn 2017, Ericeira, Portugal, 3-5 July 2107.
[2]	Balakrishnan, A.V. (1960) Fractional Powers of Closed Operators and the Semigroups Generated by Them. Pacific Journal of Mathematics, 10, 419-437. https://doi.org/10.2140/pjm.1960.10.419
[3]	Pazy, A. (1983) Semigroups of Linear Operators and Applications to Partial Differential Equations. Applied Mathematical Sciences, Vol. 44, Springer, Berlin. https://doi.org/10.1007/978-1-4612-5561-1
[4]	Engel, K.-J. and Nagel, R. (1999) One-Parameter Semigroups for Linear Evolution Equations. Graduate Texts in Mathematics, Vol. 194, Springer, Berlin.
[5]	Ashyralyev, A. (2009) A Note on Fractional Derivatives and Fractional Powers of Operators. Journal of Mathematical Analysis and Applications, 357, 232-236. https://doi.org/10.1007/978-1-4612-5561-1
[6]	Weidmann, J. (1980) Linear Operators in Hilbert Spaces. Springer, New York. https://doi.org/10.1007/978-1-4612-6027-1
[7]	Jörgens, K. and Rellich, F. (1976) Eigenwerttheorie Gewöhnlicher Differentialgleichungen. Springer, Berlin. https://doi.org/10.1007/978-3-642-66132-7
[8]	Reed, M. and Simon, B. (1975) Methods of Modern Mathematical Physics: Fourier Analysis, Self-Adjointness. Vol. 2, Academic Press, New York.
[9]	Reed, M. and Simon, B. (1978) Methods of Modern Mathematical Physics: Analysis of Operators. Vol. 4, Academic Press, New York.
[10]	Bronstein, I.N., Semendjajew, K.A., Musiol, G. and Mühlig, H. (2008) Taschenbuch der Mathematik. 7th Edition, Verlag Harri Deutsch, Frankfurt.
[11]	Davis, P.J. and Rabinowitz, P. (1984) Methods of Numerical Integration. Computer Science and Applied Mathematics, 2nd Edition, Academic Press, New York.