 # Extending discrete exterior calculus to a fractional derivative

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.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 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 Hirani2003 , DEC has been used in computer graphics applications, such as texture mapping Dominitz2010 , direction field design Vaxman2016 , fluid simulation Elcott2007 , and meshing mullen2011hot ; geometry processing applications, such as symmetry detection ben2010discrete , spectral mesh processing zhang2010spectral , mesh parameterization hormann2007mesh , and matching zaharescu2009surface ; and computational simulation, such as geometric mechanics leok2004foundations , Lie advection mullen2011discrete , stress fields in continuum mechanics kanso2007geometric , Hamiltonian PDEs bridges2006numerical , and incompressible fluids Elcott2007 .

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 that cannot be taken arbitrarily small. Such PDEs have terms involving fractional derivative operators , where is a non-integer number; the definition of these operators will be discussed later. Some recent examples of FDE-based modeling include:

• Roop Roop2006 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;

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

• Farquhar et al. Farquhar2018 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.

• Paquet and Viktor PV2013

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,

-forms are used to generalize the description of scalar fields (0-forms), vector fields (1-forms), and tensor fields (

-forms, ) defined over a smooth manifold (see e.g. AMR2012 for a formal treatment). The continuous exterior derivative operator takes -forms to ()-forms, generalizing the notions of grad, curl, and div from vector calculus. In DEC, discrete -forms are defined over a simplicial complex, with discrete -forms represented by values assigned to each -simplex in the complex; see Figure 1. The discrete exterior derivative operator takes a discrete -form to a discrete ()-form by summing the values on the boundary of a simplex with signs according to an orientation convention.

For example, the operator that could be applied to the 0-form in Figure 0(a) is a (# of edges)(# of vertices)= matrix, where the entries are either or according to edge-vertex incidence:

 D0 =⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣−110010−100−1100−10100−11⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦.

We observe that a discrete fractional derivative operator in this context should still take discrete -forms to discrete ()-forms, but must somehow incorporate more than just vertex incidence information, according to the weight 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).

• 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).

• 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 of dimension by using -dimensional simplices , where . These simplices are defined by , where the order that we give specifies the orientation of the -simplex.

On each -simplex, we can attribute some value, and this represents a discrete -form, also called a cochain. To move from -forms to -forms, we use the discrete exterior derivative , which is a rectangular matrix that has a number of columns equal to the number of -simplices in the complex and a number of rows equal to the number of -simplices in the complex. Specifically, 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 only has nonzero entries of 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 where are part of an edge (see Hirani2003 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 Figure 2: Fractional powers of order s=0, 0.2, 0.4, 0.6, 0.8, 1, for the left-sided Caputo derivative of f(x)=−10x3+10x2.

Fractional calculus has many different definitions of a fractional derivative (e.g. Riemann-Liouville, Caputo, Grunwald-Letnikov Herrmann2010 ). Central to all of them is the idea that the “th” derivative of a function for should give back a function in between the 0th (identity) and 1st derivative, varying continuously with ; 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.

A function is said to be -times continuously differentiable if we can take derivatives of and is a continuous function. In this case, we denote it as .

The Gamma function is given by

 Γ(z)=∫∞0e−ttz−1dt (1)

where . For natural numbers, the Gamma function acts like a factorial. That is, whenever . 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.

The left-sided fractional Caputo derivative of order of a function is

 Ds[a,x]f(x)=⎧⎪ ⎪⎨⎪ ⎪⎩1Γ(n−s)∫xaf(n)(τ)dτ(x−τ)s+1−ns∉Nf(s)(x)s∈N (2)

where (see Caputo1966 ). Note that often in the fractional calculus literature, the Caputo derivative is denoted as

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

A two-sided fractional Caputo derivative of order of a function is

 Dsf(x)=1Γ(n−s)(∫xaf(n)(τ)(x−τ)s+1−ndτ+∫bxf(n)(τ)(τ−x)s+1−ndτ). (3)

This can be written in a more compact form as

 Dsf(x)=1Γ(n−s)∫[a,b]∖{x}f(n)(τ)|x−τ|s+1−ndτ. (4)

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 on has the property that

 Ds[a,x]f(x)=0 (5)

for any .

###### Example 2.

Let , and . Then we have that

 Ds[0,x]f(x)=Γ(q+1)xq−sΓ(q+1−s). (6)

This result mimics what we would expect in the case of , i.e. . We see some results for the left-sided fractional Caputo derivative for in Figure 2 for a few different powers of .

By changing the subscript notation to (or more generally ), we can extend the one-variable definition of the Caputo derivative to functions of multiple variables by writing

 ∂s(∂x)s∣∣∣(ax,x)f(x,y)

with the understanding that this is just

 1Γ(1−s)∫xaxfx(τ,y)dτ(x−τ)s (7)

if .

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 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 Roop2006 , 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 mesh, he had to store a full 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:

1. That the operator maps a discrete -form with values at the vertices of a mesh to a discrete -form , taking values at the edges of the mesh. Similarly, should map a discrete -form to a discrete -form taking values at the faces.

2. maps a constant -form to at each edge of the mesh.

3. We want that as the fractional order goes to , we recover .

Note: the Caputo derivative does not yield the identity operator in the limit , 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 -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:

 Dsf(x)=1Γ(1−s)∫[a,b]∖{x}f′(τ)dτ|x−τ|s.

whenever .

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 , 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 -form , yielding a discrete -form . The discrete -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 .

Next, we discretize the integral with a kernel of . This is similar to discretizing a convolution of against the convolution kernel of . This translates to computing a weighted sum

 |E|−1∑i=01||x−ti||sD0α

where runs over each of the edges in the mesh, is the edge of interest, and represents the distance between the barycenters of the two edges.

However, this is undefined when . To fix this, we will denote a weighting vector

 (wx)i=1||x−ti||s (8)

whenever and

 (wx)i=Csmaxj≠k(1||tk−tj||s) (9)

when . Note that the constant 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 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 is a reasonable choice. Finally, we weight by . This yields the value of .

 (Ds0α)(x)=1Γ(1−s)wx⋅D0α.

We can generalize this definition as follows:

[Fractional Discrete Exterior Derivative] Let be a discrete -form on a connected simplicial complex . Then the -fractional discrete derivative of at a simplex is

 (Dspα)(σp+1)=1Γ(1−s)wσp+1⋅Dpα (10)

where is a specific simplex, is the weighting vector from above, replacing edges and with -simplices as necessary, and is the typical order discrete exterior derivative of .

## 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, will represent the number of -simplices in the mesh.

1. This definition does enact the property that correctly maps from -forms to -forms. To see this explicitly, see the equivalent matrix form of the definition that we give below.

2. We know from DEC that the dot product

 wσp+1⋅Dpα

will result in whenever takes a constant value for all -simplices in the mesh.

3. Currently, this definition does not yield a recovery of when we take . 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 be the matrix made out of using the vectors as the columns. This creates a matrix that we can use to say

 Dspα=1Γ(1−s)WDnα. (11)

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 .

The natural way that we can force the property of recovering in the limit of is to piecewise define based off the fractional order. If (or more generally, ), then define as in Equation 8, and in the case that (or ), define

 (wxk)i=δik, (12)

where is the Kronecker delta function, taking a value of if and when . One reason for doing this piecewise definition is that . 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 to recover in the limit , we would need that . In section 5, we give some evidence that as , our approximation gets better. This indicates that the matrix is indeed becoming diagonally dominant so that the fractional contribution gets less as we get close to . Hence we believe that we do achieve property .

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 . However, our definition does not satisfy an extension to . To see this, let be the corresponding weighting matrix to and notice that:

 Dsp+1Dspα (13) =1Γ(1−s)Wp+1Dp+1WpDpα. (14)

In general, and cannot commute, and the distances that are involved in will not force from any specific entries.

Part of the expectation for comes from the smooth setting, where we have . 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 Tarasov .

Finally, we must address the use of a norm in the definition. The DEC framework only assumes the existence of a local metric Hirani2003 , 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 -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 and , we know the distance where 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 , where the subscript is used to differentiate the minimum distance between any pair of vertices on the mesh from the local metric .

We extend distance between -simplices to build a distance between two -simplices. Let , be two simplices, and be the set of vertices attached to and respectively, with . Then we define the distance

 ||σp−ηp||

computing

 mini,jdm(vi,yj)+l(σp)+l(ηp)

for all pairs of with , where and 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 uses a metric that is induced by Euclidean metric of , 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 -form .

1. Create matrix.

2. Create matrix .

3. Apply the product of and to .

4. Multiply by the scaling factor .

In step we will have to do a computation that finds the distance between any two given -simplices before can be created.

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

 Ds[0,x]xq=Γ(q+1)xq−sΓ(q+1−s).

We used the left-sided fractional Caputo derivative to find this result, and assumed that with a domain of interest of . However, we based our fractional discrete exterior derivative on the 2-sided fractional Caputo derivative. In the following examples, we follow the convention that represents a left-sided Caputo derivative, while the notation represents a two-sided Caputo derivative.

###### Example 3.

As a specific example, consider on and . Then we have that

 D.5x3=Γ(4)x2.5Γ(3.5)−3Γ(3.5)(1−x)12(.75+x+2x2) (15)

as the analytic solution. Figure 3: The comparison of D.50x3 with D.5x3, for uniform meshes of various element widths. The element widths tested are given in the legend.

As we see in this example, according to Figure 3 using our definition yields a similar values. The plots of the weighted discrete -forms that result from the fractional discrete exterior derivative applied to 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 , and plot it as a constant value for the rest of the edge.

###### Example 4.

Let , , and consider the function

 f(x)=−10x3+10x2,

with a fractional Caputo derivative of

 D.5f(x)=−10Γ(4)Γ(72)x52+10Γ(3)Γ(52)x32−30Γ(72)(1−x)12(34+x+2x2)+20Γ(52)(1−x)12(x+12). (16)

Since the function we are working with has an easily calculated closed form for its fractional Caputo derivative and a discrete -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.

Edges error
2 1.5619
4 0.9933
8 .6778
16 .4759
32 .3363
64 .2378
128 .1681
256 .1188
512 .0839
1024 .0593

The error results here give the indication that our error is going to 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. Figure 4: The comparison of D.50(−10x3+10x2) with D.5(−10x3+10x2). We use varying step sizes denoted in the legend.

Since this example has a simple analytic solution to the fractional Caputo derivative we can test the effect of the value of .

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 , our error can fluctuate within this range. Part of this error is expected because of the constant based on . Figure 5: We compare the fractional exponent s against the L∞ error for the function f(x)=−10x3+10x2. We use a step size varying step sizes given in the legend.

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

The Mittag-Leffler two parameter function is

 Ea,b(z)=∞∑k=0zkΓ(ak+b) (17)

where and . In general, we can think of this function as being a generalization of the exponential function. In fact, for the values of , 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 . Explicitly, we get

 Ds[0,x]ex=x1−sE1,2−s(x) (18)

whenever Mariya .

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 . In the left-sided fractional Caputo derivative, we have an integration domain of . To account for this, we need to modify the weighting matrix . Specifically, each entry of represents the weight assigned to one simplex in relation to another simplex, say and . Then, in the -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 and . If , we leave as it is and otherwise we set . The results of this comparison are given in Figure 6. Figure 6: Comparison of the Caputo derivative of ex versus D.50ex.

In Figure 6, we see the difference in trajectories for the Caputo derivative of . The code used to generate the Mittag-Leffler function is from Garrappa2015 . 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 norm. We see the results of this analysis in Figure 7. Figure 7: Error analysis for Ds0ex. We use the L∞ norm to find max errors at varying step sizes in comparison to x.5E1,1.5(x), the closed form of the left-sided Caputo derivative of ex.

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

###### Example 5.

Let . Then the 2-sided fractional gradient field on the region of order is

 ⎛⎜ ⎜⎝−2x32Γ(52)−2(1−x)12(x+12)Γ(52),2y32Γ(52)+2(1−y)12(y+12)Γ(52)⎞⎟ ⎟⎠. (19)

Using this, we consider a triangulated plot of 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 -form vector that represents the value of at the vertices of the mesh. After applying the fractional discrete exterior derivative, we get a discrete -form , 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 (Hirani2003, , Definition 3.3.4). The construction is as follows: compute the barycentric coordinates , , locally on each triangle, use these to compute the three Whitney 1-forms () 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 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. Figure 9: The triangles are colored to represent the relative L2 error for f(x,y)=−x2+y2. We get error values ranging between .0627 and 4.4227, with an average of 1.2126.
###### Example 6.

We do one last example. Let

 f(x,y)=(x−.1)2+(y−.1)2 (20)

be defined on . In this case, we know that has a minimum at . The 2-sided fractional Caputo gradient field of is

 D12f=(2x32−(1−x)12(12−2x)Γ(52)−x12−(1−x)125Γ(32),2y32−(1−y)12(12−2y)Γ(52)−y12−(1−y)125Γ(32)) (21)

This fractional gradient field is plotted in Figure 10. We also give a plot for relative error for 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 . 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. Figure 10: The 2-sided fractional Caputo gradient field for the function f(x,y)=(x−.1)2+(y−.1)2. Figure 11: The triangles are colored to represent the relative L2 error for f(x)=(x−.1)2+(y−.1)2. Errors here range from .2369 to 2.7504, with an average error of 1.3561.

## 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 - and -d examples. Our definition requires an interpretation of distance between two arbitrary -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 . 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. Tarasov ; CS2001 . This includes a definition of a fractional exterior derivative, with the form

 ds=n∑i=1∂s∂xis[dxi]s.

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 , as is shown in CS2001 . 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, . The Laplacian operator has a nice definition in DEC where . Hence studying the fractional Laplacian in the DEC setting may reduce to studying , with 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 Lischke2018 .

As discussed in the introduction, there was also a shape analysis application by Paquet and Viktor in PV2013 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

 RLDs[a,x]f(x)=1Γ(n−s)dndxn∫xaf(t)(x−t)n+s−1dt, (22)

where is some domain of interest, and .

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

 1Γ(1−s)DpWα (23)

where is a weighting matrix similar to before, but instead of acting on discrete -forms, it acts on discrete -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.

## References

• (1) A. N. Hirani, Discrete Exterior Calculus, Ph.D. thesis, Caltech (2003).
• (2) A. Dominitz, A. Tannenbaum, Texture mapping via optimal mass transport, IEEE transactions on visualization and computer graphics 16 (3) (2010) 419–433.
• (3) A. Vaxman, M. Campen, O. Diamanti, D. Panozzo, D. Bommes, K. Hildebrandt, M. Ben-Chen, Directional field synthesis, design, and processing, in: Computer Graphics Forum, Vol. 35, Wiley Online Library, 2016, pp. 545–572.
• (4) S. Elcott, Y. Tong, E. Kanso, P. Schröder, M. Desbrun, Stable, circulation-preserving, simplicial fluids, ACM Transactions on Graphics (TOG) 26 (1) (2007) 4.
• (5) P. Mullen, P. Memari, F. de Goes, M. Desbrun, HOT: Hodge-optimized triangulations, in: ACM Transactions on Graphics (TOG), Vol. 30, ACM, 2011, p. 103.
• (6) M. Ben-Chen, A. Butscher, J. Solomon, L. Guibas, On discrete killing vector fields and patterns on surfaces, in: Computer Graphics Forum, Vol. 29, Wiley Online Library, 2010, pp. 1701–1711.
• (7) H. Zhang, O. Van Kaick, R. Dyer, Spectral mesh processing, in: Computer graphics forum, Vol. 29, Wiley Online Library, 2010, pp. 1865–1894.
• (8) K. Hormann, B. Lévy, A. Sheffer, Mesh parameterization: Theory and practice, presentation at SIGGRAPH (2007).
• (9)

A. Zaharescu, E. Boyer, K. Varanasi, R. Horaud, Surface feature detection and description with applications to mesh matching, in: Computer Vision and Pattern Recognition, 2009. CVPR 2009. IEEE Conference on, IEEE, 2009, pp. 373–380.

• (10) M. Leok, Foundations of computational geometric mechanics, Ph.D. thesis, California Institute of Technology (2004).
• (11) P. Mullen, A. McKenzie, D. Pavlov, L. Durant, Y. Tong, E. Kanso, J. E. Marsden, M. Desbrun, Discrete Lie advection of differential forms, Foundations of Computational Mathematics 11 (2) (2011) 131–149.
• (12) E. Kanso, M. Arroyo, Y. Tong, A. Yavari, J. G. Marsden, M. Desbrun, On the geometric character of stress in continuum mechanics, Zeitschrift für angewandte Mathematik und Physik 58 (5) (2007) 843–856.
• (13) T. J. Bridges, S. Reich, Numerical methods for Hamiltonian PDEs, Journal of Physics A: mathematical and general 39 (19) (2006) 5287.
• (14) J. P. Roop, Computational aspects of FEM approximation of fractional advection dispersion equations on bounded domains in R2, Journal of Computational and Applied Mathematics 193 (1) (2006) 243–268.
• (15) O. Ozgen, M. Kallmann, L. E. Ramirez, C. F. Coimbra, Underwater cloth simulation with fractional derivatives, ACM Transactions on Graphics (TOG) 29 (3) (2010) 1–9.
URL http://graphics.ucmerced.edu/papers/10-tog-fracdef.pdf
• (16) M. E. Farquhar, T. J. Moroney, Q. Yang, I. W. Turner, K. Burrage, Computational modelling of cardiac ischaemia using a variable-order fractional laplacian, arXiv:1809.07936 (2018).
• (17) E. Paquet, H. L. Viktor, Isometrically invariant description of deformable objects based on the fractional heat equation, in: R. Wilson, E. Hancock, A. Bors, W. Smith (Eds.), Computer Analysis of Images and Patterns, Springer Berlin Heidelberg, Berlin, Heidelberg, 2013, pp. 135–143.
• (18) R. Abraham, J. E. Marsden, T. Ratiu, Manifolds, tensor analysis, and applications, Vol. 75, Springer Science & Business Media, 2012.
• (19) R. Herrmann, Fractional Calculus An Introduction for Physicists, 2010.
• (20) M. Caputo, Linear models of dissipation whose Q is almost frequency independent-II, Geophysical Journal International 13 (5) (1967) 529–539.
• (21) V. E. Tarasov, Fractional dynamics: applications of fractional calculus to dynamics of particles, fields and media, Springer Science & Business Media, 2011.
• (22) M. K. Ishteva, Properties and Applications of the Caputo Fractional Operator (2005).
• (23) R. Garrappa, Numerical evaluation of two and three parameter Mittag-Leffler functionsarXiv:1503.06569, doi:10.1137/140971191.
• (24) K. Cottrill-Shepherd, M. Naber, Fractional differential forms, Journal of Mathematical Physics 42 (5) (2001) 2203–2212.
• (25) A. Lischke, G. Pang, M. Gulian, F. Song, C. Glusa, X. Zheng, Z. Mao, W. Cai, M. M. Meerschaert, M. Ainsworth, G. E. Karniadakis, What Is the Fractional Laplacian? (2018) 1–80arXiv:1801.09767.
URL http://arxiv.org/abs/1801.09767