Fractional Operators

The Dunford-Taylor representation of regularly accretive operators (here Laplacian) $$ (-\Delta)^{-s} = \frac{1}{2\pi i}\int_{\mathcal C} z^{-s} (z+\Delta)^{-1}dz, \qquad 0 < s < 1, $$ leads to numerical method for the approximation of the spectral and integral fractional laplacian.

Spectral Fractional Laplacian

Let \(\Omega\) be a bounded Lipschitz domain. The spectral fractional Laplacian \((-\Delta)^s\) with \(s\in (0,1)\) is defined for \(v\in \{w\in L^2(\Omega) : (-\Delta)^s w\in L^2(\Omega) \}\) by $$(-\Delta)^s v = \sum_{j=0}^\infty \lambda_j^s (v,\psi_j) \psi_j,$$ where \((\cdot,\cdot)\) denotes the \(L^2\) inner product and \(\{\psi_j\}\) is an orthonormal basis of eigenfunctions of \(-\Delta\) correspondin g to eigenvalues \(\{\lambda_j\}\).

Direct Approximation of the Inverse

We denote by \(\mathcal H^s(\Omega) := (L^2(\Omega), H^1_0(\Omega))_{s,2}\) the interpolation space (real method). Given \( f \in L_2(\Omega)\), we seek to approximate \(u \in \mathcal H^s(\Omega)\) satisfying $$(-\Delta)^{s} u = f \quad \textrm{in } \Omega.$$ The proposed approximation method relies on the so-called Balakrishnan formula $$ u = (-\Delta)^{-s} f = \frac{c_s}{2} \int_{0}^{\infty} \mu^{-s} (\mu I-\Delta)^{-1}f \,d \mu = c_s \int_{-\infty}^{\infty} e^{2st} (I-e^{2t}\Delta)^{-1}f \,d t,$$ where \(c_s:={2\sin(\pi s)}/{\pi}\), and consists in two steps:

  • Apply a sinc quadrature scheme $$ u \approx u^{k} := c_s k\sum_{j = -N_{-}}^{N_{+}} e^{2st_j} (I-e^{2t_j}\Delta)^{-1} f, $$ where \(k\) is the quadrature step size, \(t_j = jk\) are the quadrature points, \(N_{+} := \left\lceil \frac{2\pi^2}{sk^2} \right\rceil, \;\text{and}\; N_{-} := \left\lceil \frac{4\pi^2}{(1-s) k^2} \right\rceil\), with \(\lceil\cdot\rceil\) denoting the round up to the nearest integer;
  • Use standard finite element method to approximate \( w_j:=(I-e^{2t_j}\Delta)^{-1} f\) with \(w_{h,j}\in V_h\) for each \(j\) on the same mesh by solving $$ (w_{h,j}, v_h)_\Omega + e^{2t_j} (\nabla w_{h,j}, \nabla v_h)_\Omega = (f,v_h)_\Omega, \qquad \forall v \in V_h, $$ where \((\cdot, \cdot)_\Omega\) denotes the \(L^2\) inner product on \(\Omega\), and \(V_h \subset H^1_0(\Omega)\) is the finite element space consisting in continuous piecewise linear functions.

The resulting \(N_- + N_+ + 1\) finite element problems are mutually independent, making the parallel implementation straightforward. In addition, for each subproblem it only requires the implementation of a classical finite element solver for the diffusion-reaction problem.

The discrepancy between the exact solution and its approximation consists of the sinc quadrature error, which is exponentially convergent, and the finite element approximation error, which is optimal up to a logarithmic factor depending on the smoothness of the data \(f\). It is worth mentioning that the error analysis does not require the domain to be convex.

Approximation of the Scalar Product $a_s(u,\varphi)=\int_\Omega (-\Delta)^{s/2}u (-\Delta)^{s/2}\varphi$

For \(u,\varphi\in \mathcal H^s(\Omega)\), the inner product also has an integral representation $$a_s(u,\varphi):=\int_\Omega (-\Delta)^{s/2} u (-\Delta)^{s/2}\varphi = \frac{c_s}{2} \int_{-\infty}^{\infty} e^{st} \int_\Omega \Big((-\Delta)(e^tI - {\Delta})^{-1} u\Big)\varphi dt.$$ The numerical approximation to the inner product follows essentially the same two steps as in the spectral fractional operator approximation.

The error estimation is achieved by analyzing the consistency error between \(a_s(u,\varphi)\) and its two-step numerical approximation. The consistency error is the sum of the exponentially convergent sinc quadrature error and the finite element error. A Strang's lemma type argument is used to derive error estimates.

Integral Fractional Laplacian

The integral fractional Laplacian is defined using the Fourier transform. For \(v\) in the Schwartz space, the integral fractional Laplacian is a pseudo-differential operator with symbol \(|\zeta|^{2s}\), $$\mathscr{F} ((-\Delta)^s v)(\zeta) = |\zeta|^{2s}\mathscr{F}(v)(\zeta),$$ with \(\mathscr{F}\) denoting the Fourier transform. The definition can be extended to the Sobolev space \(\mathcal H^s(\mathbb R^d)\) by a density argument.

Approximation of integral fractional Laplacian

Consider the following model problem involving the integral fractional Laplacian: given \(f\in L^2(\Omega)\), find \(u\in \mathcal H^s(\Omega)\) satisfying $$ (-\Delta)^s \tilde{u} = f \qquad\;\text{in}\; \Omega, $$ where \(\tilde u\) denotes the extension by zero outside the domain \(\Omega\). Unlike the spectral fractional Laplacian, we could not find a useful integral representation of \(\tilde u\). Instead, the approximation relies on a formulation of the scalar product similar to the spectral case $$ \int_{\mathbb R^d} (-\Delta)^{s/2} \tilde u (-\Delta)^{s/2}\tilde\varphi = \frac{c_s}{2} \int_{-\infty}^{\infty} e^{st} \int_\Omega \Big( (-\Delta)(e^tI-\Delta)^{-1}\tilde u \Big) \varphi~ dt, $$ with \(u,\varphi\) in \(\mathcal H^s(\Omega)\). This integral representation suggests the following approximation strategy similar to the one for the spectral Laplacian above:

  • Apply a sinc quadrature scheme with stepping \(k\) resulting in $$ \tilde a_s(u,\varphi) \approx \frac{c_s}{2} k \sum_{j = -N_-}^{N_+} \int_\Omega \Big(e^{st_j}(-\Delta)(e^{t_j}I-\Delta)^{-1}\tilde u \Big) \varphi, $$ where \(\{t_j\}=\{jk\}\), \(j = -N_-, \cdots, N_+\), for \(N_-,N_+ \propto 1/k^2\) but depending on the regularity of the solution \(u\). Notice that at each quadrature point, the solution to an elliptic problem in \(\mathbb R^d\) is required;
  • For a truncation parameter \(M\), truncate the problems in \(\mathbb R^d\) to a ball of radius depending on the quadrature point \(j\);
  • Approximate using a finite element method the resulting elliptic PDEs supplemented with vanishing Dirichlet conditions.

The discrepancy between the exact solution and its fully discrete approximation consists of three parts: an exponentially convergent sinc quadrature error, an exponentially convergent domain truncation error, and a finite element approximation error. The algorithm does not require \(\Omega\) to be convex and works on three dimensional domain as illustrated in the figure below.

Numerical approximation of the solution to the fractional integral Laplacian for \(s = 0.3\) and \(f = 1\) in the unit ball of \(\mathbb R^3\) (darker = 0.0, whiter = 0.7). The lines represent the isosurfaces \(\{u(x) = k/10\}\) for \(k = 0,\dots,7\) and reflect the singular behavior of the solution towards the domain boundary.

Collaborators and Relevant Literature

  • BONITO, A. AND LEI, W.; Approximation of the spectral fractional powers of the Laplace-Beltrami Operator. Numerical Mathematics: Theory, Methods and Applications, 15(4), 1193--1218, 2022.
  • BONITO, A. AND GUIGNARD, D. AND ZHANG, A.; Reduced basis approximations of the solutions to fractional diffusion problems. J. Numer. Math., 28(3), 147--160, 2020.
  • BONITO, A. AND LEI, W. AND PASCIAK, J.E; Numerical Approximation of the Integral Fractional Laplacian. Numer. Math, 142(2), 235--278, 2019
  • BONITO, A. AND LEI, W. AND PASCIAK, J.E; On Sinc Quadrature Approximations of Fractional Powers of Regularly Accretive Operators. J. Numer. Math, accepted.
  • BONITO, A. AND BORTHAGARAY, J.P. AND NOCHETTO, R.H AND OTÁROLA, E. AND SALGADO , A.J; Numerical Methods for Fractional Diffusion. Comput. Vis. Sci., 19(5), 19--46, 2018.
  • BONITO, A. AND LEI, W. AND PASCIAK, J.E; Numerical Approximation of space-time fractional parabolic equations. Comput. Methods Appl. Math., 17(4), 679--705, 2017.
  • BONITO, A. AND LEI, W. AND PASCIAK, J.E; The Approximation of Parabolic Equations Involving Fractional Powers of Elliptic Operators, J. Comput. Appl. Math., 315, 32--48, 2017.
  • BONITO, A. AND PASCIAK, J.E; Numerical Approximation of Fractional Powers of Regularly Accretive Operators, IMA J. Numer. Anal., 37(3), 1245--1273, 2017.
  • BONITO, A. AND PASCIAK, J.E; Numerical Approximation of Fractional Powers of Elliptic Operators, Math. Comp., 84 (295), 2137--2162, 2015.

The electrical convection, or electroconvection in short, is the phenomenon caused by the application of a sufficiently large voltage difference across a thin-layer fluid such that the unbalanced charges are pushed by the applied electric field to form pairs of vortices of flows.

In the physical experiments by Tsai et al. (2007), at the room temperature the liquid crystal is at smectic-A stage so that the cigar-shaped molecules are arranged into layers. Within each layer, the molecules are aligned perpendicular to the plane of the film. Under such arrangement, the molecules are free to move along the film while movements normal to the layers are prevented. This leads to a two dimensional motion.

This thin-layered film is placed in between two concentric electrodes with the boundaries touching the electrodes. The inner electrode is charged with positive voltage \(V\) and the outer one is grounded. Inside the film, the positive charges are gathering near the inner boundary, while negative charges gather around the outer boundary. Such unstable equilibrium can be easily broke by a small perturbation on the charge distribution. As a consequence, counter-rotating vortex flows emerge when \(V\) exceeds an threshold value \(V_c\).

As in Tsai et al. (2007), we consider two settings: the electrodes extend to infinity in the normal direction or the electrodes have negligible thickness, see Figure 1. The former corresponds to the simulation reported in (Tsai 2007, Tsai et al. 2007) while the later is related to the analysis by Constantin et. al. (2016). A model reduction to the fluid domain introduces nonlocal representation for the electric potential \(\Psi:\Omega \to \mathbb R\) either involving the spectral fractional Laplacian (infinite electrodes): $$ (-\Delta)^{1/2} \Psi = q \quad \textrm{in } \Omega, $$ or the integral fractional Laplacian (slim electrodes): $$ (-\Delta)^{1/2} \tilde\Psi = q + g \quad \textrm{in } \Omega. $$ Here \(q:\Omega \to \mathbb R\) is the liquid surface charge density, and \(g:\Omega \to \mathbb R\) is an axis-symmetric function accounting for the voltage difference.

Experimental settings. The liquid film is located in between two concentric electrodes. Left: the two electrodes extend to infinity (infinite case) and Right: the two electrodes have negligible thickness (slim case). Notice that in both cases, the outer electrode is extended to infinity in the plane (not pictured).

Numerical Simulations

The electric potential \(V\) and the aspect ratio between the inner radius and outer radius of the film are critical parameters for electroconvection to develop.

Comparison between the infinite and slim electrodes models

The difference between the two models is strikingly significant. To illustrate this claim, we consider an aspect ratio of \(0.33\), a Rayleigh number of \(\mathcal R=800\), a Prandtl number of \(\mathcal P = 10\) and provide in Figure 2 a comparison of the evolution of kinetic and circulation energies $$ E_k := \frac{1}{2}\int_{\Omega} |\mathbf u|^2 \;d\mathbf x, \qquad E_{curl} := \int_{\Omega} |\nabla\times \mathbf u|^2 \;d\mathbf x, $$ where \(\mathbf u\) denotes the velocity field of the fluid.

Comparison of the evolution the energies versus time for the slim and infinite electrodes configurations with same parameter configuration. The slim configuration energies are significantly larger and remain large during the entire evolution indicating a sustained electroconvection phenomena.
Slim model - Effect of the aspect ratio

We investigate the effect of the geometry.

The velocity field of the fluid (left) and the surface charge density distribution (right) at developped convection.
Aspect ratio = 0.33, Rayleigh = 100, Prandtl = 10
Aspect ratio = 0.56, Rayleigh = 100, Prandtl = 10

Collaborators and Relevant Literature

  • BONITO, A. AND WEI, P.; Electroconvection of Thin Liquid Crystals: Model Reduction and Numerical Simulations. J. Comput. Phys, 405, 2020.

The general quasi-geostrophic model, first introduced by J. G. Charney in the 1940s, has been very successful for the study of oceanic and atmospheric dynamics in the mid-to-high latitude region of the earth where the Coriolis effect is significant. The characteristic property of the quasi-geostrophic system is the conservation of potential vorticity along the geostrophic flow. If in addition to the potential vorticity conservation, the assumption on the potential temperature conservation reduces the model to the surface quasi-geostrophic equation (Constantin et al. 1994,Held et al. 1995)

The governing equation is a transport equation with a fractional dissipation term: $${\partial_t \theta} + \mathbf{u}\cdot \nabla \theta + \kappa(-\Delta)^{s}\theta = 0, \quad \textbf{u} = \nabla^{\bot}\psi, \quad (-\Delta)^{1/2} \psi= -\theta \quad\text{ all in }\; \Omega.$$ Here \(\Omega\) is a square domain sitting on the \(xy\) plane, the potential temperature \(\theta = \theta(\mathbf{x},t)\) is a scalar function defined on \(\Omega\times[0,\texttt T]\) with \(\texttt T>0\), \(\nabla^{\bot}\psi := (-\partial_y \psi, \partial_x \psi)^T\), and \(\psi:\Omega\times[0,\texttt T]\to \mathbb R\) is the geopotential. The fractional power \(s\) belongs to \([1/2,1)\), and the coefficient \(\kappa \ge 0\) is a constant. Both \(\psi\) and \(\theta\) are supplemented with periodic boundary condition, and further fixed by an averaging constraint. Both fractional operators are spectral fractional Laplacians. The non-local property of the fractional term is the main difficulty in both analysis and numerical approximation.

The expression of the velocity \(\textbf{u}\) corresponds to the motion along the iso-bars of \(\psi\). This is due to the fact that the gradient pressure is balanced by the Coriolis force in the gradient direction. For instance, for a high-pressure area in the southern hemisphere, the pressure gradient is pointing outward, while the Coriolis force is perpendicular to the direction of the movement to the left side, leading to a counter-clockwise rotation to balance the two forces. The term \(\kappa(-\Delta)^{s}\theta\) is the so-called Ekman pumping, representing a balance between the vertical component of the Coriolis force and vertical frictional force. Typically, \(\kappa > 0\) and \(s = 1/2\).

Initial buoyancy \(\theta(x,y) = \sin(x)\sin(y) + \cos(y)\) containing a saddle structure (nessecary for singularity to develop). Viscid-limit case \(\kappa=0.001\) and \(s = 0.5\).
Image Image Image Image Image
Vortex simulation. Initially \(\theta(x,y) = \exp\left(-(x-\pi)^2-16(y-\pi)^2\right)\) the buoyancy turns into a spinning vortex with thin filaments. Small vertices appear. It resembles cyclonic circulations within the atmosphere.
Image Image Image Image Image Image
Double vortex simulation. A sharp layer develop between the two vortices around $t=8$. Although the intensity of the layer separating the two vortices is reducing over time, the vortices do not merge.
Image Image Image Image Image Image
Freely decaying turbulance. White noise initial condition organizes into larger and larger vortices.

Collaborators and Relevant Literature

  • BONITO, A. AND NAZAROV, M.; Numerical Simulations of Surface-Quasi Geostrophic Flows on Periodic Domains. SIAM Journal of Scientific Computing, 43(2), B405-B430, 2021.

This material is based upon work partially supported by the National Science Foundation. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.