Extending discrete exterior calculus to a fractional derivative

TL;DR

This paper extends discrete exterior calculus to fractional derivatives, enabling more accurate numerical solutions for fractional PDEs involving non-local interactions, by defining a fractional discrete derivative based on DEC principles.

Contribution

It introduces a Caputo-like fractional discrete derivative within DEC, bridging the gap between fractional calculus and geometric discretization methods.

Findings

The fractional discrete derivative retains key properties of the continuous derivative.

Numerical experiments demonstrate the effectiveness of the proposed discretization.

The method provides a new computational tool for fractional PDEs.

Abstract

Fractional partial differential equations (FDEs) are used to describe phenomena that involve a "non-local" or "long-range" interaction of some kind. Accurate and practical numerical approximation of their solutions is challenging due to the dense matrices arising from standard discretization procedures. In this paper, we begin to extend the well-established computational toolkit of Discrete Exterior Calculus (DEC) to the fractional setting, focusing on proper discretization of the fractional derivative. We define a Caputo-like fractional discrete derivative, in terms of the standard discrete exterior derivative operator from DEC, weighted by a measure of distance between $p$-simplices in a simplicial complex. We discuss key theoretical properties of the fractional discrete derivative and compare it to the continuous fractional derivative via a series of numerical experiments.

Equations

$${\mathbb{D}}_0$$

$$\Gamma (z)={\int }_0^{\infty }{e}^{-t}{t}^{z-1}dt$$

$$D_{[a,x]}^{s}f(x)=\{\begin{array}{ll}\frac{1}{\Gamma (n-s)}{\int }_a^x\frac{f^{(n)}(\tau )d\tau }{(x-\tau {)}^{s+1-n}}& s\notin \mathbb{N}\\ f^{(s)}(x)& s\in \mathbb{N}\end{array}$$

$${}_a^C{D}_x^{s}f(x).$$

$$D^{s}f(x)=\frac{1}{\Gamma (n-s)}\left({\int }_a^x\frac{f^{(n)}(\tau )}{(x-\tau {)}^{s+1-n}}d\tau +{\int }_x^b\frac{f^{(n)}(\tau )}{(\tau -x{)}^{s+1-n}}d\tau \right).$$

$$D^{s}f(x)=\frac{1}{\Gamma (n-s)}{\int }_{[a,b]\backslash \{x\}}\frac{f^{(n)}(\tau )}{|x-\tau {|}^{s+1-n}}d\tau .$$

$$D_{[a,x]}^{s}f(x)=0$$

$$D_{[0,x]}^{s}f(x)=\frac{\Gamma (q+1)x^{q-s}}{\Gamma (q+1-s)}.$$

$${\frac{\partial^{s}}{(\partial x)^{s}}}\bigg{\rvert}_{(a_{x},x)}f(x,y)$$

$$\frac{1}{\Gamma (1-s)}{\int }_{a_x}^x\frac{f_x(\tau ,y)d\tau }{(x-\tau {)}^s}$$

$$D^{s}f(x)=\frac{1}{\Gamma (1-s)}{\int }_{[a,b]\backslash \{x\}}\frac{f^{\prime }(\tau )d\tau }{|x-\tau {|}^s}.$$

$$\sum _{i=0}^{|E|-1}\frac{1}{||x-t_i|{|}^s}{\mathbb{D}}_0\alpha$$

$$(w_x{)}_i=\frac{1}{||x-t_i|{|}^s}$$

$$(w_x{)}_i=C_s{\text{max}}_{j\ne k}\left(\frac{1}{||t_k-t_j|{|}^s}\right)$$

$$\left({\mathbb{D}}_0^s\alpha \right)(x)=\frac{1}{\Gamma (1-s)}{w}_x\cdot {\mathbb{D}}_0\alpha .$$

$$\left({\mathbb{D}}_p^s\alpha \right)({\sigma }^{p+1})=\frac{1}{\Gamma (1-s)}{w}_{{\sigma }^{p+1}}\cdot {\mathbb{D}}_p\alpha$$

$$w_{{\sigma }^{p+1}}\cdot {\mathbb{D}}_p\alpha$$

$${\mathbb{D}}_p^s\alpha =\frac{1}{\Gamma (1-s)}W{\mathbb{D}}_n\alpha .$$

$$(w_{x_k}{)}_i={\delta }_{ik},$$

$${\mathbb{D}}_{p+1}^s{\mathbb{D}}_p^s\alpha$$

$$=\frac{1}{\Gamma (1-s)}{W}_{p+1}{\mathbb{D}}_{p+1}{W}_p{\mathbb{D}}_p\alpha .$$

$$||{\sigma }^p-{\eta }^p||$$

$$\underset{i,j}{\min }{d}_m(v_i,y_j)+l({\sigma }^p)+l({\eta }^p)$$

$$D_{[0,x]}^s{x}^q=\frac{\Gamma (q+1)x^{q-s}}{\Gamma (q+1-s)}.$$

$$D^{.5}{x}^3=\frac{\Gamma (4)x^{2.5}}{\Gamma (3.5)}-\frac{3}{\Gamma (3.5)}(1-x{)}^{\frac{1}{2}}(.75+x+2x^2)$$

$$f(x)=-10x^3+10x^2,$$

$$\begin{array}{rl}{D}^{.5}f(x)=& -\frac{10\Gamma (4)}{\Gamma \left(\frac{7}{2}\right)}{x}^{\frac{5}{2}}+\frac{10\Gamma (3)}{\Gamma \left(\frac{5}{2}\right)}{x}^{\frac{3}{2}}-\\ & \frac{30}{\Gamma \left(\frac{7}{2}\right)}(1-x{)}^{\frac{1}{2}}\left(\frac{3}{4}+x+2x^2\right)+\frac{20}{\Gamma \left(\frac{5}{2}\right)}(1-x{)}^{\frac{1}{2}}(x+\frac{1}{2}).\end{array}$$

$$E_{a,b}(z)=\sum _{k=0}^{\infty }\frac{z^k}{\Gamma (ak+b)}$$

$$D_{[0,x]}^s{e}^x=x^{1-s}{E}_{1,2-s}(x)$$

$$\left(\frac{-2x^{\frac{3}{2}}}{\Gamma \left(\frac{5}{2}\right)}-\frac{2(1-x{)}^{\frac{1}{2}}(x+\frac{1}{2})}{\Gamma \left(\frac{5}{2}\right)},\frac{2y^{\frac{3}{2}}}{\Gamma \left(\frac{5}{2}\right)}+\frac{2(1-y{)}^{\frac{1}{2}}(y+\frac{1}{2})}{\Gamma \left(\frac{5}{2}\right)}\right).$$

$$f(x,y)=(x-.1{)}^2+(y-.1{)}^2$$

Full text

1 Introduction

Discrete exterior calculus (DEC) is a computational toolkit based on the principle of creating discrete analogues of objects and operators from smooth differential geometry. Since the original thesis work of Hirani [1], DEC has been used in computer graphics applications, such as texture mapping [2], direction field design [3], fluid simulation [4], and meshing [5]; geometry processing applications, such as symmetry detection [6], spectral mesh processing [7], mesh parameterization [8], and matching [9]; and computational simulation, such as geometric mechanics [10], Lie advection [11], stress fields in continuum mechanics [12], Hamiltonian PDEs [13], and incompressible fluids [4].

Absent from these applications is a consideration of how to apply the theory in the context of fractional partial differential equations (FDEs). A PDE is called “fractional” if it involves a non-local interaction of some kind, such as a time-dependent problem with “memory” or a spatial variable whose value at a point depends on values within some radius $\epsilon$ that cannot be taken arbitrarily small. Such PDEs have terms involving fractional derivative operators $({\partial}/{\partial x})^{s}$, where $s$ is a non-integer number; the definition of these operators will be discussed later. Some recent examples of FDE-based modeling include:

Roop [14] studied models of “anomalous diffusion,” in which a fractional advection-dispersion equation gave a more accurate description of fluid flow through porous media than a standard diffusion equation;

• 2.

Ozgen et al. [15] simulated cloth movement underwater using an FDE with a drag term, which gave visibly better results than standard techniques;

• 3.

Farquhar et al. [16] modeled electric currents through a heart with tissue affected by ischaemia using a mix of standard and fractional versions of the monodomain equation; they vary the fractional order across the computational domain according to the local level of blood supply.

• 4.

Paquet and Viktor [17] studied shape analysis by using a fractional deRham operator. They used an eigenvalue decomposition to compute the fractional power of their operator, and found that the nonlocal nature of the fractional operator better allowed them to capture the shape of an object.

While many additional applications of FDEs can be found in the literature, an accurate, generalized, and efficient approach to their discretization remains elusive. Our long term goal is to provide such an approach via the tools and framework provided by discrete exterior calculus. As a starting point, in this paper, we introduce certain fundamental operators that would come up in a discrete approximation of a solution to a FDE. Specifically, we focus on techniques for discretizing fractional derivative operators in accordance with DEC principles and conventions.

A brief review of key concepts from smooth and discrete exterior calculus helps motivate our approach. In continuous exterior calculus, $p$-forms are used to generalize the description of scalar fields (0-forms), vector fields (1-forms), and tensor fields ($p$-forms, $p\geq 2$) defined over a smooth manifold (see e.g. [18] for a formal treatment). The continuous exterior derivative operator takes $p$-forms to ($p+1$)-forms, generalizing the notions of grad, curl, and div from vector calculus. In DEC, discrete $p$-forms are defined over a simplicial complex, with discrete $p$-forms represented by values assigned to each $p$-simplex in the complex; see Figure 1. The discrete exterior derivative operator $\mathbb{D}_{p}$ takes a discrete $p$-form to a discrete ($p+1$)-form by summing the values on the boundary of a $p+1$ simplex with signs $\pm 1$ according to an orientation convention.

For example, the $\mathbb{D}_{0}$ operator that could be applied to the 0-form in Figure 1(a) is a (# of edges)$\times$(# of vertices)=$5\times 4$ matrix, where the entries are either [math] or $\pm 1$ according to edge-vertex incidence:

[TABLE]

We observe that a discrete fractional derivative operator in this context should still take discrete $p$-forms to discrete ($p+1$)-forms, but must somehow incorporate more than just vertex incidence information, according to the weight $s$ of the fractional derivative that is involved. Toward that goal, we present the following results:

Summarize relevant background from DEC and fractional calculus (Section 2).

• 2.

Give a computable definition for the fractional derivative that adheres to the structure laid out in DEC (Section 3) and investigate its properties (Section 4).

• 3.

Test our definition on 1D and 2D examples (Section 5) and discuss its effectiveness and possible future applications (Section 6).

2 Background and Notation

In this section, we will give some background information to both DEC and fractional calculus.

2.1 Discrete Exterior Calculus

To use DEC, we build up a simplicial complex $K$ of dimension $N$ by using $p$-dimensional simplices $\sigma^{p}$, where $p\in\{0,1,\ldots,N\}$. These simplices are defined by $\sigma^{p}=[v_{0},v_{1}\ldots,v_{p}]$, where the order that we give $v_{0},v_{1},\ldots v_{p}$ specifies the orientation of the $p$-simplex.

On each $p$-simplex, we can attribute some value, and this represents a discrete $p$-form, also called a cochain. To move from $p$-forms to $(p+1)$-forms, we use the discrete exterior derivative $\mathbb{D}_{p}$, which is a rectangular matrix that has a number of columns equal to the number of $p$-simplices in the complex and a number of rows equal to the number of $(p+1)$-simplices in the complex. Specifically, $\mathbb{D}_{p}$ is the transpose of the boundary operator matrix.

The discrete exterior derivative is a local operator that does not make use of any metric information. The matrix $\mathbb{D}_{p}$ only has nonzero entries of $\pm 1$ for simplices that are adjacent. However, other DEC operators (such as the Hodge star) make use of metric information that is given by a local metric $d(v_{0},v_{1})$ where $v_{0},v_{1}$ are part of an edge $[v_{0},v_{1}]\in K$ (see [1] a discussion on this). These values give information about the distances between vertices on the primal mesh. As we will see later, this will become important when extending the discrete exterior derivative to a fractional setting.

2.2 Fractional Calculus

Fractional calculus has many different definitions of a fractional derivative (e.g. Riemann-Liouville, Caputo, Grunwald-Letnikov [19]). Central to all of them is the idea that the “$s$th” derivative of a function for $s\in[0,1]$ should give back a function in between the 0th (identity) and 1st derivative, varying continuously with $s$; see Figure 2. In our work, we will focus on the Caputo derivative. Before we give the definition of the Caputo derivative itself, first we give some background information.

Definition**.**

A function $f$ is said to be $n$-times continuously differentiable if we can take $n$ derivatives of $f$ and $f^{(n)}$ is a continuous function. In this case, we denote it as $f\in C^{n}[a,b]$.

Definition**.**

The Gamma function is given by

[TABLE]

where $Re(z)>0$.

For natural numbers, the Gamma function acts like a factorial. That is, $\Gamma(n+1)=n!$ whenever $n\in\mathbb{N}$. For our purposes, the Gamma function will typically just be a constant that we can compute when needed. Now we have the fractional Caputo derivative.

Definition**.**

The left-sided fractional Caputo derivative of order $s\in\mathbb{R}^{+}$ of a function $f\in C^{n}[a,b]$ is

[TABLE]

where $n=\lceil s\rceil$ (see [20]). Note that often in the fractional calculus literature, the Caputo derivative is denoted as

[TABLE]

We choose to not use this notation because the majority of this work will focus only on the Caputo derivative. The only time we deviate from this is in Section 6, where we will briefly discuss the Riemann-Liouville derivative. In that case, we will make sure it is clear that we are discussing a fractional derivative operator that is not the Caputo derivative.

The subscript $[a,x]$ serves two purposes. The first is that it defines the domain of integration for the fractional derivative, and the second is that it is used to denote which variable this is with respect to. We will also make use of a 2-sided definition of a fractional Caputo derivative.

Definition**.**

A two-sided fractional Caputo derivative of order $s\in\mathbb{R}^{+}\backslash\mathbb{Z}$ of a function $f\in C^{n}[a,b]$ is

[TABLE]

This can be written in a more compact form as

[TABLE]

To understand the behavior of these operators, we have some examples in the continuous setting of the left-sided fractional Caputo derivative.

Example 1**.**

Due to taking a derivative before integrating, any constant function $f(x)=c$ on $[0,1]$ has the property that

[TABLE]

for any $s\in\mathbb{R}$.

Example 2**.**

Let $f(x)=x^{q}$, $q>0$ and $s\in(0,1)$. Then we have that

[TABLE]

This result mimics what we would expect in the case of $s=1$, i.e. $D^{1}_{[0,x]}x^{q}=qx^{q-1}$. We see some results for the left-sided fractional Caputo derivative for $f(x)=-10x^{3}+10x^{2}$ in Figure 2 for a few different powers of $s$.

By changing the subscript notation to $[a_{x},x]$ (or more generally $[a_{i},x_{i}]$), we can extend the one-variable definition of the Caputo derivative to functions of multiple variables by writing

[TABLE]

with the understanding that this is just

[TABLE]

if $s\in(0,1)$.

Regardless of how many variables we choose to work in, part of the difficulty with FDEs becomes obvious from the definitions: the operators that we are using are non-local. Since these derivatives require the integration over some domain $[a,x]$, a small neighborhood of information at a point is not enough to give the value of the fractional derivative at that point. This leads to some numerical difficulties.

Primarily, a major issue is that numeric computations are more challenging. For spatial FDEs, examining a large neighborhood manifests as more computational work to gather and process the surrounding information. In general, for every fractional exponent, we may need to examine the values defined on the entire domain to accurately compute the derivative. For fractional derivatives that involve time, the computation necessitates storing historical information for all previous values. This represents a process with memory, but there is a tradeoff in that storing all history in practice would be computationally infeasible.

In our paper we will be focusing on spatial FDEs, but this manifests as its own set of computational problems. Consider using a matrix to encode the derivative operator. Due to the non-local nature of a fractional derivative, such matrices lose a lot of the sparsity and structure that they would normally have. For example, in [14], Roop investigated the use of a finite element scheme to solve the fractional advection-dispersion equation. In order to use quadratic basis elements on a $64\times 64$ mesh, he had to store a full $16641\times 16641$ stiffness matrix with no special structure.

Spatial FDEs are a primary factor in choosing to use a 2-sided fractional Caputo operator. Using a left-sided operator in a temporal setting can be interpreted as meaning the process has a memory effect. A left-sided operator in a spatial setting would ignore non-local effects that occur on the right side with no geometrical justification. Thus the 2-sided definition we have given will be used in our discretization for defining a fractional discrete exterior derivative.

We will see a similar issue of having to store large, dense matrices later in this work. Our main goal at this stage is the definition of a fractional discrete exterior derivative. Addressing how to effectively and accurately do this computations without losing the information the operator is trying to provide will be considered in future work.

3 Defining a fractional discrete exterior derivative

Combining fractional calculus and DEC comes with its own set of challenges. Part of the design of DEC is that it is a set of tools that work in a local setting, with no need to understand global distances.

Our goal is to create a definition that imitates the fractional Caputo derivative. While doing this, our definition for fractional discrete exterior derivative should also act in a similar fashion to the normal discrete exterior derivative map. Thus we have some properties that we seek to preserve:

That the operator $\mathbb{D}_{0}^{s}$ maps a discrete [math]-form $f$ with values at the vertices of a mesh to a discrete $1$-form $\mathbb{D}_{0}^{s}f$, taking values at the edges of the mesh. Similarly, $\mathbb{D}_{1}^{s}$ should map a discrete $1$-form to a discrete $2$-form taking values at the faces. 2. 2.

$\mathbb{D}_{0}^{s}$ maps a constant [math]-form to [math] at each edge of the mesh. 3. 3.

We want that as the fractional order goes to $1$, we recover $\mathbb{D}_{k}$.

Note: the Caputo derivative does not yield the identity operator in the limit $s\rightarrow 0$, so it is not a property that we try to emulate.

We mention this here before illustrating the design of our definition. In choosing to use DEC to tackle fractional differential equations, we must make certain design choices that are going to lead to complicating portions of DEC that are normally straightforward. We believe that the extra complications are worth it to make use of the framework that DEC gives us in solving differential equations. In particular, the ability to computing distances between $p$-simplices is a non-trivial requirement for the computation. We will discuss this more in section 4.

3.1 The Fractional Discrete Exterior Derivative

As we build our definition, we will start with the Caputo derivative definition and replace smooth operators with DEC operators wherever possible. For convenience, we restrict ourselves to the setting where the two-sided fractional Caputo derivative is:

[TABLE]

whenever $s\in(0,1)$.

Before we begin discretizing, notice that the Caputo derivative operator can be decomposed into a sequence of standard operations working inside out. First we differentiate $f$, then we integrate a weighted version of the derivative. Thus we will end up discretizing this in the same sequence of steps.

To discretize, we first apply a discrete exterior derivative to a [math]-form $\alpha$, yielding a discrete $1$-form $\mathbb{D}_{0}\alpha$. The discrete $1$-form should take values at each edge, so for the remaining steps we will focus on what the fractional operator computes at a specific edge $x$.

Next, we discretize the integral with a kernel of $\dfrac{1}{|x-t|^{s}}$. This is similar to discretizing a convolution of $\mathbb{D}_{0}\alpha$ against the convolution kernel of $\dfrac{1}{|x-t|^{s}}$. This translates to computing a weighted sum

[TABLE]

where $t_{i}$ runs over each of the edges in the mesh, $x$ is the edge of interest, and $||x-t_{i}||$ represents the distance between the barycenters of the two edges.

However, this is undefined when $t_{i}=x$. To fix this, we will denote a weighting vector

[TABLE]

whenever $x\neq t_{i}$ and

[TABLE]

when $x=t_{i}$. Note that the constant $C_{s}$ is used to give more weight to the current edge than the other edges in the mesh. Choosing how to properly define this entry in the vector $w_{x}$ is an interesting problem by itself, and our definition is sensitive to the choice of this value. We found in our numerical experiments that that $C_{s}=\frac{2s}{1-s}$ is a reasonable choice. Finally, we weight by $\frac{1}{\Gamma(1-s)}$. This yields the value of $(\mathbb{D}_{0}^{s}\alpha)(x)$.

[TABLE]

We can generalize this definition as follows:

Definition** (Fractional Discrete Exterior Derivative).**

Let $\alpha$ be a discrete $p$-form on a connected simplicial complex $M$. Then the $(p,s)$-fractional discrete derivative of $\alpha$ at a $(p+1)$ simplex $\sigma^{p+1}$ is

[TABLE]

where $\sigma^{p+1}$ is a specific $p+1$ simplex, $w_{\sigma^{p+1}}$ is the weighting vector from above, replacing edges $x_{k}$ and $t_{i}$ with $(p+1)$-simplices as necessary, and $\mathbb{D}_{p}\alpha$ is the typical $p^{th}$ order discrete exterior derivative of $\alpha$.

4 Properties of the fractional discrete exterior derivative

Before we get to testing the definition on some examples, we proceed with a brief discussion on the properties of the above subsection 3.1. In this discussion, $|P|$ will represent the number of $p$-simplices in the mesh.

This definition does enact the property that $\mathbb{D}_{p}^{s}$ correctly maps from $p$-forms to $(p+1)$-forms. To see this explicitly, see the equivalent matrix form of the definition that we give below. 2. 2.

We know from DEC that the dot product

[TABLE]

will result in $\mathbb{D}_{p}^{s}\alpha=0$ whenever $\alpha$ takes a constant value for all $p$-simplices in the mesh. 3. 3.

Currently, this definition does not yield a recovery of $\mathbb{D}_{p}$ when we take $s\rightarrow 1$. This can be fixed, and we will show what this may look like below.

Note that we can rewrite the definition of the fractional discrete exterior derivative operator using matrices. To do this, let $W$ be the matrix made out of using the vectors $w_{x_{k}}$ as the columns. This creates a $|P|\times|P|$ matrix that we can use to say

[TABLE]

The choice to start this way was made to make it easier to see the action at a specific simplex. However, this form makes it clear that $\mathbb{D}_{p}^{s}:\Omega_{d}^{p}\rightarrow\Omega_{d}^{p+1}$.

The natural way that we can force the property of recovering $\mathbb{D}_{p}$ in the limit of $s\rightarrow 1$ is to piecewise define $w_{x_{k}}$ based off the fractional order. If $s\in(0,1)$ (or more generally, $s\in\mathbb{R}^{+}\backslash\mathbb{N}$), then define $w_{x_{k}}$ as in Equation 8, and in the case that $s=1$ (or $s\in\mathbb{N}$), define

[TABLE]

where $\delta_{ik}$ is the Kronecker delta function, taking a value of [math] if $i\neq k$ and $1$ when $i=k$. One reason for doing this piecewise definition is that $\Gamma(0)=0$. Therefore both the fractional Caputo derivative and the fractional discrete exterior derivative cannot be defined here by the usual definition or the new definition, respectively. Furthermore, in order for $\mathbb{D}_{p}^{s}$ to recover $\mathbb{D}_{p}$ in the limit $s\rightarrow 1$, we would need that $W\rightarrow\Gamma\left(1-s\right)I$. In section 5, we give some evidence that as $s\rightarrow 1$, our approximation gets better. This indicates that the $W$ matrix is indeed becoming diagonally dominant so that the fractional contribution gets less as we get close to $s=1$. Hence we believe that we do achieve property $3$.

There are other properties that we chose not to enforce that may have been natural given the setting of DEC. For example, the composition property $\mathbb{D}_{p+1}\circ\mathbb{D}_{p}=0$. However, our definition does not satisfy an extension to $\mathbb{D}_{p+1}^{s}\circ\mathbb{D}_{p}^{s}=0$. To see this, let $W_{k}$ be the corresponding weighting matrix to $\mathbb{D}_{k}^{s}$ and notice that:

[TABLE]

In general, $W_{p}$ and $\mathbb{D}_{p}$ cannot commute, and the distances that are involved in $W_{k}$ will not force [math] from any specific entries.

Part of the expectation for $\mathbb{D}_{p+1}\circ\mathbb{D}_{p}=0$ comes from the smooth setting, where we have $d^{p+1}\circ d^{p}=0$. In other words, all exact forms in the smooth setting are closed. While a similar property would be good to have, we do not have a clear reason for expecting it to be true after the weighting process occurs. For additional formal discussion on the smooth theory of fractional exterior calculus, see chapter 12 of [21].

Finally, we must address the use of a norm $||\cdot||$ in the definition. The DEC framework only assumes the existence of a local metric [1], which can be used to measure distance only between adjacent simplices. We next discuss how to extend the local metric to be able to compute the distance between two arbitrary $p$-simplices.

To do this, we first need a well-defined way to find the distance between two arbitrary vertices on our mesh. We make use of the fact that between any adjacent vertices $u$ and $v$, we know the distance $d(u,v)$ where $d$ is the local metric on the mesh. Then we can apply an all pairs-shortest path algorithm to find the minimum distance between each pair of vertices on the path. This gives a distance between any pair of vertices, not just adjacent ones, which we denote $d_{m}(u,v)$, where the subscript $m$ is used to differentiate the minimum distance between any pair of vertices on the mesh from the local metric $d$.

We extend distance between [math]-simplices to build a distance between two $p$-simplices. Let $\sigma^{p}$, $\eta^{p}$ be two simplices, $\{u_{k}\}$ and $\{v_{k}\}$ be the set of vertices attached to $\sigma^{p}$ and $\eta^{p}$ respectively, with $k=0,1,\ldots,p$. Then we define the distance

[TABLE]

computing

[TABLE]

for all pairs of $i,j$ with $i,j=0,1,\ldots,p$, where $l(\sigma^{p})$ and $l(\eta^{p})$ is the length between the barycenter and the boundary of each of the respective simplices.

As an aside, if we know that the embedding of the simplicial complex $K$ uses a metric that is induced by Euclidean metric of $\mathbb{R}^{N}$, then it is easier to implement the weighting matrix by using the Euclidean distance between circumcenters or barycenters of the two simplices. That is, if we know that there exists a global metric that correctly embeds the complex, we can use the global metric.

5 Numerical experiments

With our definition in hand, we test it in some controlled scenarios. First though, we give a brief outline for how one would go about implementing the fractional discrete exterior derivative as defined above, assuming we already have a discrete $p$-form $\alpha$.

Create $\mathbb{D}_{p}$ matrix. 2. 2.

Create matrix $W$. 3. 3.

Apply the product of $W$ and $\mathbb{D}_{p}$ to $\alpha$. 4. 4.

Multiply by the scaling factor $\displaystyle\frac{1}{\Gamma(1-s)}$.

In step $2$ we will have to do a computation that finds the distance between any two given $(p+1)$-simplices before $W$ can be created.

Now we do an example. Recall from Equation 6 the form of the Caputo derivative acting on power functions

[TABLE]

We used the left-sided fractional Caputo derivative to find this result, and assumed that $s\in(0,1)$ with a domain of interest of $[0,1]$. However, we based our fractional discrete exterior derivative on the 2-sided fractional Caputo derivative. In the following examples, we follow the convention that $D^{s}_{[0,x]}$ represents a left-sided Caputo derivative, while the notation $D^{s}$ represents a two-sided Caputo derivative.

Example 3**.**

As a specific example, consider $x^{3}$ on $[0,1]$ and $s=.5$. Then we have that

[TABLE]

as the analytic solution.

As we see in this example, according to Figure 3 using our definition yields a similar values. The plots of the weighted discrete $1$-forms that result from the fractional discrete exterior derivative applied to $x^{3}$ on the mesh are obtained using the stairs plot in Matlab. This plotting method results in a staircase function where a given step starts at a barycenter of an edge and lasts until the next step. For each of these steps, we use the computed value in the vector $\mathbb{D}^{s}_{0}x^{3}$, and plot it as a constant value for the rest of the edge.

Example 4**.**

Let $s=0.5$, $[a,b]=[0,1]$, and consider the function

[TABLE]

with a fractional Caputo derivative of

[TABLE]

Since the function we are working with has an easily calculated closed form for its fractional Caputo derivative and a discrete $1$-form can be thought of as piecewise constant on edges, we are able to symbolically integrate to get our error values. The results of the error analysis are given in the table below.

[TABLE]

The error results here give the indication that our error is going to [math] as the step size decreases. In Figure 4, we compare the closed-form version against another piecewise linear estimate that we get as the result of applying the fractional discrete exterior derivative. In the case of different step sizes, we see that we do get closer to the true trajectory of the fractional derivative.

Since this example has a simple analytic solution to the fractional Caputo derivative we can test the effect of the value of $s\in(0,1)$.

In FDEs, it is commonplace to test the error against the fractional powers. Some applications have a clear reasoning, whether experimental or physical, for desiring a fractional power to be some specific value. Other applications treat the fractional power as free parameter. Our results in Figure 5 show that while we can change the fractional power as we desire in $(0,1)$, our error can fluctuate within this range. Part of this error is expected because of the $\frac{1}{\Gamma(1-s)}$ constant based on $s$.

Finally, we want to check how our definition compares to the left-sided fractional Caputo derivative of $e^{x}$. To that end, we give the following definition and then a brief discussion of the modification of $W$ to account for a left-sided fractional derivative instead of a 2-sided fractional derivative.

Definition**.**

The Mittag-Leffler two parameter function is

[TABLE]

where $a,b\in\mathbb{R}^{+}$ and $z\in\mathbb{C}$.

In general, we can think of this function as being a generalization of the exponential function. In fact, for the values of $a=1,b=1$, we recover the exponential function exactly. The Mittag-Leffler function can be used to give a nice looking form for the result of the fractional Caputo derivative of $e^{x}$. Explicitly, we get

[TABLE]

whenever $s\in(0,1)$ [22].

Our discretization of the fractional discrete exterior derivative is based off the two-sided fractional Caputo derivative. However, we want to be able to compare against $D_{[0,x]}^{s}e^{x}$. In the left-sided fractional Caputo derivative, we have an integration domain of $[0,x]$. To account for this, we need to modify the weighting matrix $W$. Specifically, each entry of $W$ represents the weight assigned to one simplex in relation to another simplex, say $\sigma_{i}$ and $\sigma_{j}$. Then, in the $1$-d case, to determine if the entry should be zero to account for the left-sided version of the derivative we compare the x-coordinates of the barycenters of $\sigma_{i}$ and $\sigma_{j}$. If $x_{i}<x_{j}$, we leave $W_{i,j}$ as it is and otherwise we set $W_{i,j}=0$. The results of this comparison are given in Figure 6.

In Figure 6, we see the difference in trajectories for the Caputo derivative of $e^{x}$. The code used to generate the Mittag-Leffler function is from [23]. Our definition undershoots the value of the Caputo derivative, but does seem to exhibit the correct behavior. To do an accurate error analysis in this case, we choose the $L_{\infty}$ norm. We see the results of this analysis in Figure 7.

The following examples are efforts in testing the fractional discrete exterior derivative in $2$-d. We start with a basic function

Example 5**.**

Let $f(x,y)=-x^{2}+y^{2}$. Then the 2-sided fractional gradient field on the region $M=[0,1]\times[0,1]$ of order $s\in(0,1)$ is

[TABLE]

Using this, we consider a triangulated plot of $M$ with the gradient vector field drawn at the barycenters of the triangles. This field is plotted in Figure 8.

To test the fractional discrete exterior derivative on this example and visualize a vector field requires some extra work. We first create the [math]-form vector $\alpha$ that represents the value of $f$ at the vertices of the mesh. After applying the fractional discrete exterior derivative, we get a discrete $1$-form $\mathbb{D}^{.5}\alpha$, i.e. a cochain storing a value for each edge of the mesh with edge orientation implied by the global vertex ordering. Visualizing and measuring the error of this data type requires some care.

We first construct the Whitney map for 1-cochains, based on [1, Definition 3.3.4]. The construction is as follows: compute the barycentric coordinates $\lambda_{1}$, $\lambda_{2}$, $\lambda_{3}$ locally on each triangle, use these to compute the three Whitney 1-forms ($\lambda_{i}\nabla\lambda_{j}-\lambda_{j}\nabla\lambda_{i}$) associated to each edge, then weight each Whitney 1-form by the cochain value on the corresponding edge. The result is a vector field, defined piecewise over the entire mesh, that recovers the cochain values on each edge (when projected onto an edge). We evaluate this global vector field at the barycenter of each triangle for both qualitative and quantitative error analysis.

To give a sense of what the fractional Caputo gradient field on $M$ looks like, see Figure 8. In Figure 9, we look at relative error. Notice that originally, the function has a saddle type critical point (i.e. the gradient of the function is zero) at the origin. This critical point, along with some boundary effects from our definition, lead to a high relative error near the origin.

Example 6**.**

We do one last example. Let

[TABLE]

be defined on $M=[0,1]\times[0,1]$. In this case, we know that $f(x,y)$ has a minimum at $(.1,.1)$. The 2-sided fractional Caputo gradient field of $f$ is

[TABLE]

This fractional gradient field is plotted in Figure 10. We also give a plot for relative error for $f$ in Figure 11. In this case, we again have a critical point that affects our results. For this function, the critical point exists in the domain $M$. One point of interest is that the fractional Caputo derivative does not preserve the location of critical points, which is what causes a large amount of the error in Figure 11. In either setting, understanding the behavior of the fractional discrete exterior derivative on functions of multiple variables with critical points is a topic that we intend to investigate further in the future.

6 Discussion and conclusions

In this paper we gave a new definition of the fractional discrete exterior derivative that satisfied three key properties. Our definition was based on the fractional Caputo derivative, and we tested on some $1$- and $2$-d examples. Our definition requires an interpretation of distance between two arbitrary $p$-simplices on the mesh.

There are some limitations in the fractional discrete exterior derivative definition. First, it complicates the process of taking a discrete exterior derivative, which is normally represented by a single matrix with a nice structure. While this does mean that the process of taking a fractional discrete exterior derivative is more involved, this is to be expected in comparison to a normal discrete exterior derivative. Specifically, the non-local nature of the fractional derivative means that we have to dedicate time in the computation to determining distances between simplices that are far apart. The other limitation is that our definition does not act on critical points the same way that the fractional Caputo derivative does. This means that near critical points, we end up with a much higher error term than we would like.

For future work, we have many avenues to investigate. First, we would like to study the way that the fractional discrete exterior derivative responds to tweaking the diagonal of the matrix $W$. Furthermore, we would like to be able to understand the behavior better at the critical points of a function.

There is also a connection that we have not yet discussed. In the past couple of decades, there has been work on defining an extension of exterior calculus to a fractional exterior calculus setting, e.g. [21, 24]. This includes a definition of a fractional exterior derivative, with the form

[TABLE]

We intentionally chose not to make use of definition this when defining our operator. Directly discretizing this fractional exterior derivative would not fit in the “discrete first” approach that DEC takes, and discretizing it in a way analogous to how the discrete exterior derivative is discretized would require a notion of a fractional boundary. One possible strength of discretizing with this formula is the possibility of ensuring $d^{s}\circ d^{s}=0$, as is shown in [24]. Still, fractional exterior calculus as a whole is a relatively new field, and it is not clear which properties are essential to preserve when considering their discretization.

Beyond just investigating more of the properties of the definition, we want to extend into a more DEC related setting. The results on an interval and a rectangular region subdivided into triangles are promising, however one strength of DEC is its ability to work on manifold-like simplicial complexes. With a definition of the fractional discrete exterior derivative, we should be able to extend current work using DEC to FDEs on simplicial complexes that represent more intricate geometric objects.

In the vein of solving FDEs, one of the prominent operators in FDEs is the fractional Laplacian, $(-\Delta)^{s}$. The Laplacian operator has a nice definition in DEC where $\Delta=*d*d$. Hence studying the fractional Laplacian in the DEC setting may reduce to studying $\star\mathbb{D}_{p-1}^{s}\star\mathbb{D}_{0}^{s}$, with $\star$ being some sort of fractional discrete Hodge star operator. A good overview of the work that has been done on the fractional Laplacian was written by Lischke et. al in [25].

As discussed in the introduction, there was also a shape analysis application by Paquet and Viktor in [17] using the fractional de Rham operator. In their approach, they take a fractional power by using an eigenvalue decomposition, and raising the eigenvalues to the fractional power. In the fractional Laplacian literature, this is referred to as the spectral method for computing the fractional power of the operator. This would be interesting to test against when we are ready to use our method to try discretizing the fractional Laplacian. In general, we cannot use the spectral method for the fractional discrete exterior derivative because it is a non-square matrix.

We would also be interested in exploring what the other definitions of fractional derivatives yield when discretized through DEC. For example, the Riemann-Liouville fractional derivative is given by

[TABLE]

where $[a,b]$ is some domain of interest, $s\in\mathbb{R}^{+}$ and $n=\lceil s\rceil$.

The main difference between this definition and the Caputo definition is the order of the differentiation and integration. Because the integration occurs before taking any derivatives, the class of functions we act on is slightly less restrictive. We suspect that the Riemann-Liouville definition leads to a fractional discrete exterior derivative operator that has the form

[TABLE]

where $W$ is a weighting matrix similar to before, but instead of acting on discrete $(p+1)$-forms, it acts on discrete $p$-forms.

Having a working version of the fractional discrete exterior derivative derived from the Riemann-Liouville derivative will be useful. Many applications of FDEs require either the Caputo derivative or the Riemann-Liouville derivative, and they in general do not yield the same results. Enabling the computation of both forms will allow for more tests that can be used to compare and evaluate fractional derivatives.

We are also interested in using these operators to study fractional equivalents to problems that have been previously investigated with DEC. For example, recent work has been done on spatial Darcy flow problems to numerically compare results against experimental data.

Acknowledgements

This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, under Award Number(s) DE-SC-0019039.