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

$\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)

$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