# Cantilevered, Rectangular Plate Dynamics
by Finite Difference Methods

Benjamin Brown
University of Maryland, Baltimore County
Baltimore, MD
<EMAIL_ADDRESS>
###### Abstract

In this technical note, we consider a dynamic linear, cantilevered rectangular
plate. The evolutionary PDE model is given by the fourth order plate dynamics
(via the spatial biharmonic operator) with clamped-free-free-free boundary
conditions. We additionally consider damping/dissipation terms, as well as
non-conservative lower order terms arising in various applications. Dynamical
numerical simulations are achieved by way of a finite difference spatial
approximation with a MATLAB time integrator. The rectangular geometry allows
the use of standard 2D spatial finite differences, while the high spatial
order of the problem and mixed clamped-free type boundary conditions present
challenges. Dynamic energies are also computed. The relevant code is
presented, with discussion of the model and context.

##  1 Introduction

The purpose of this note is to describe a dynamic simulation method for a
linear 2D cantilevered plate via finite differences. We specifically determine
the ghost points required for the enforcement of the cantilever boundary
conditions (mixed clamped-free with corners).

###  1.1 Motivation

The primary motivation for studying the particular PDE model below derives
from the application of piezoelectric energy harvesting [14, 13]. In this
application, a cantilevered beam or rectangular plate is placed in a
surrounding flow of air. In the axial configuration (whereby the flow runs
from clamped end to free end of the structure), a cantilever is particularly
prone to an aeroelastic instability termed flutter, even for low flow speeds.
After the onset of instability, the beam settles into a limit cycle
oscillation. If a piezo-electric device (for instance a thin layer, or patch)
is affixed to the deflecting structure, power is generated and can be
harvested. The energy is produced by the system’s sustained large, flapping
motion which causes a current to be produced via the piezoelectric effect.
Such energy-harvesting configurations have recently been shown to be viable
[14]. As an alternative energy source, scaled harvester systems could be
effective in providing power to remote locations, or small surges of power to
supplement a traditional grid.

To effectively and efficiently harvest energy from the flow-induced
oscillations of an elastic cantilever, one must have a viable distributed
parameter system describing the large deflections of the 1D or 2D cantilevered
structure. In particular, to produce a limit cycle oscillation in the post-
flutter regime (i.e., after the onset of the flow-induced bifurcation), the
structural model must permit large deflections and incorporate a nonlinear
restoring force. Recent work has 1D inextensible cantilevers, including
modeling [8, 7], and more recently, a well-posedness analysis of solutions
[5]. More recently, the modeling work was extended to (several) 2D
inextensible cantilevered plates [6]. These plate models are very involved,
including nonlinear inertial and stiffness terms, as well as Lagrange
multipliers for the enforcement of relevant constraints. In this note, we
focus on the linearization of the model in [6], and provide a fast numerical
solver for a plate with given forcing (through the boundary and the right hand
side). We accommodate damping effects, and crude approximations of the flow
via a linear piston-theoretic approximation of potential flow (see, e.g.,
[10]). This model is can be used viably in the pre-onset regime, as well as to
compute a flow-induced bifurcation point via the non-conservative flow
parameters. Additionally, the code is general enough to compute dynamic
responses to arbitrary edge loading, as well as dynamic and spatially
distributed right hand sides.

Much numerical work has been done on linear cantilevered models in 1D [11, 10,
7], however much less mathematical literature exists for 2D cantilevered plate
models, likely owing to the challenges associated to the high order of the
problem and requisitely mixed (with corners) boundary conditions111The mixed
nature of the boundary conditions, as well as the corner angles, are known to
affect the expected solution regularity—and also numerical convergence—near
corners [3].. We call attention to one older technical note that discusses
finite difference stencils for rectangular plates with different boundary
conditions [2]—the analysis there is entirely static. It is clear, in general,
that 2D models (and associated simulations), are necessary for the engineering
applications mentioned above. We resolve the PDE dynamics below using finite
difference methods with careful calculations of ghost nodes that allow
algebraic corner conditions to be properly solved.

###  1.2 Setup and Model

An open rectangular domain, $\Omega$ with boundary
$\Gamma=\overline{\Gamma_{E}}\cup\overline{\Gamma_{N}}\cup\overline{\Gamma_{W}}\cup\overline{\Gamma_{S}}$,
will be used for this model, and will use cardinal coordinates for naming, as
can be seen in the following diagram:

The linear cantilevered plate model, with the rectangular domain $\Omega$, is
described by the following equations:

$\begin{cases}w_{tt}+D\Delta^{2}w+k_{0}w_{t}-k_{1}\Delta
w_{t}+a_{1}w_{x}+a_{2}w_{y}=f(x,t)~{}~{}~{}&\text{in }\Omega\\\\[5.0pt]
w=\frac{\partial w}{\partial n}=0~{}~{}~{}&\text{on }\Gamma_{W}\\\\[5.0pt] \nu
w_{xx}+w_{yy}=g_{1}(x,y),~{}~{}w_{yyy}+(2-\nu)w_{xxy}=h_{1}(x,y)~{}~{}~{}&\text{on
}\Gamma_{N}\\\\[5.0pt] \nu
w_{xx}+w_{yy}=g_{2}(x,y),~{}~{}w_{yyy}+(2-\nu)w_{xxy}=h_{2}(x,y)~{}~{}~{}&\text{on
}\Gamma_{S}\\\\[5.0pt] w_{xx}+\nu
w_{yy}=g_{3}(x,y),~{}~{}w_{xxx}+(2-\nu)w_{yyx}=h_{3}(x,y)~{}~{}~{}&\text{on
}\Gamma_{E}\\\\[5.0pt] w(t=0)=w_{0}(x),~{}w_{t}(t=0)=w_{1}(x)\end{cases}$ (1)

The quantity $\nu\in(0,1/2)$ represents the Poisson ratio [6]. Above, the
function $w:\overline{\Omega}\times[0,T)\to\mathbf{R}$ represents the
transverse displacement of an thin, isotropic, homogeneous plate. We do not
explicitly incorporate piezoelectric effects here, and, assuming a
piezoelectric material is attached to the surface of the plate, we incorporate
the combined the combined stiffness effect through the parameter $D>0$. The
damping coefficients are $k_{0},k_{1}\geq 0$, where the former measures weak
or frictional damping, and the latter is so-called square-root type damping
and is of a stronger nature. (See the discussion in [11] for further
discussion about the interpretation and physical meaning of this type of
damping.) In reality, the term $-k_{1}\Delta w_{t}$ term may have questionable
meaning for the cantilever, but is left in the implementation as it is
unproblematic mathematically; we set $k_{1}=0$ for all examples shown
throughout this note. The flow-velocity parameters are $a_{1},a_{2}$
corresponding the respective $x$ and $y$ components of the surrounding
airflow. A linear combination of these, for instance, represents an off-axis
flow; pure axial flow would take $a_{1}>0$ and $a_{2}=0$. The functions
$g_{i}(x,y)$ correspond to edge-loading by way of boundary moments, while the
boundary shears are given by $h_{i}(x,y)$ on each of the non-clamped
boundaries in the model. The function $f(x,t)$ provides a distributed loading
across the surface of $\Omega$. In practice, for regular solutions, these
boundary and interior data must satisfy natural compatibility conditions, in
particular at/near corner points.

Due to the fourth order bi-harmonic operator, each point in the mesh created
on the plate must be capable of having a 13-point stencil applied, including
the edges and corners, such as in Figure 1. Owing to this, ghost points must
be added around the boundary of the plate. To wit, a ghost point or ghost node
is an artificial value for a point that created outside of $\overline{\Omega}$
to enforce the boundary conditions, permitting a finite difference stencil to
be applied when the stencil overlaps boundary nodes. The particular challenges
in the resolution and regularity of a given solution at corner points
manifests themselves in algebraic issues for finite difference methods. Ghost
points are calculated through using the boundary conditions, as well as known
values from the interior, both of which are obtained through the data in the
problem (e.g., $f$, $g_{i}$, and/or $h_{i}$). Indeed, the ghost node
calculations are a central contribution of this note.

$O$$N$$NN$$S$$SS$$WW$$W$$E$$EE$$NW$$NE$$SW$$SE$ Figure 1: The 13-point stencil
used in the finite difference method to approximate a fourth derivative

##  2 Ghost Points

Figure 2: Diagram showing the order of ghost points being calculated around a
free-free corner. Begin by calculating the first row, which are the green
nodes. Then the second row will be calculated, represented by yellow nodes.
Finally the system of seven equations can be solved for the corner nodes,
shown in red.

As alluded to above, the inherent mixed boundary conditions in the cantilever
problem introduce complications into the calculation of the ghost points. To
determine ghost node values, we will begin with the clamped edge, because it
is simple and only has one unknown. This will provide information to supply to
more stencils in latter calculations of nodal values.

The nomenclature used in this section makes use of cardinal values that
correspond to the current displacement of the point in the stated direction,
with $O$ representing the ‘origin’. For example, $E$ refers to the mesh point
to the right of $O$.

###  2.1 Clamped Edge

On this edge, the conditions require that $w=0$ and $w_{x}=0$. Using the
finite difference approximation on this point:

$E=W,$

for all $O$ along the clamped boundary.

###  2.2 Free Edge

All of the following calculations will be made for the $\Gamma_{N}$ boundary,
but they apply to the other edges with a simple rotation of coordinates as
needed.

The conditions on the three free edges require specified moments and shears:

$\begin{cases}\nu
w_{xx}+w_{yy}=g_{1}(x,y),~{}~{}w_{yyy}+(2-\nu)w_{xxy}=h_{1}(x,y)~{}~{}~{}&\text{on
}\Gamma_{N}\\\\[5.0pt] \nu
w_{xx}+w_{yy}=g_{2}(x,y),~{}~{}w_{yyy}+(2-\nu)w_{xxy}=h_{2}(x,y)~{}~{}~{}&\text{on
}\Gamma_{S}\\\\[5.0pt] w_{xx}+\nu
w_{yy}=g_{3}(x,y),~{}~{}w_{xxx}+(2-\nu)w_{yyx}=h_{3}(x,y)~{}~{}~{}&\text{on
}\Gamma_{E}\\\\[5.0pt] \end{cases}$

Because a 13-point stencil will be needed for the final calculation, two rows
of ghost points must be determined.

####  2.2.1 First Row on Edge: Avoiding Corners

The second order boundary condition will be used for the first row. The free-
free corner will be avoided for this calculation, as it has has special
calculations, shown later. Using the finite difference method to approximate
this condition yields the following: expression:

${\frac{\nu\,\left(W-2\,O+E\right)}{{{dx}}^{2}}}+{\frac{N-2\,O+S}{{{dy}}^{2}}}=g_{1}(x,y).$

Notice that along all the free edges, except when adjacent to a free-free
corner, there is only one unknown value, the ghost value. This can be
explicitly solved for, for example on the southern edge, to give the requisite
value of the ghost point:

$S=-{\frac{\nu\,\left(W-2\,O+E\right){{dy}}^{2}}{{{dx}}^{2}}}-N+2\,O+g_{1}(x,y).$

####  2.2.2 Second Row on Edge: Avoiding Corners

The second boundary condition will be used for the second row of ghost points.
When this stencil is utilized, and the first edge has already been calculated,
the second row ghost values can then be solved directly. This calculations
uses ghost nodes found in the previous step in the first row. Because of this,
there must again be caution about the corner. The nodes to be calculated are
shown in Figure 2 as the yellow nodes. Using the finite difference
approximation on the second boundary condition yields the following
expression:

${\frac{{SS}-2\,S+2\,N-{NN}}{2\,{{dy}}^{3}}}+{\frac{\left(2-\nu\right)\left({SE}-{NE}-2\,S+2\,N+{SW}-{NW}\right)}{2\,{{dx}}^{2}{dy}}}=h_{1}(x,y).$

The ghost point $SS$ can be determined from the above, as it is the only
unknown when applied to the southern free edge. The expression for this value
is:

${SS}=-{\frac{\left(2-\nu\right){{dy}}^{2}\left({SE}-{NE}-2\,S+2\,N+{SW}-{NW}\right)}{{{dx}}^{2}}}+2\,S-2\,N+{NN}+h_{1}(x,y).$

####  2.2.3 Corners

For each of the two free-free corners, there are seven unknown ghost points.
We can make a linear system of the seven boundary conditions that can be
applied to the surrounding points in order to explicitly find each of them. In
this section, for simplicity, the $g_{i}(x,y)$ and $h_{i}(x,y)$ terms will be
set to zero. In a latter section we will elaborate upon the edge loading
conditions.

The relevant finite difference boundary conditions are:

$\begin{cases}{\frac{\nu\,\left(W-2\,O+E\right)}{{{dx}}^{2}}}+{\frac{N-2\,O+S}{{{dy}}^{2}}}=0\\\\[5.0pt]
{\frac{W-2\,O+E}{{{dx}}^{2}}}+{\frac{\nu\,\left(N-2\,O+S\right)}{{{dy}}^{2}}}=0\\\\[5.0pt]
{\frac{{SS}-2\,S+2\,N-{NN}}{2\,{{dy}}^{3}}}+{\frac{\left(2-\nu\right)\left({SE}-{NE}-2\,S+2\,N+{SW}-{NW}\right)}{2\,{{dx}}^{2}{dy}}}=0\\\\[5.0pt]
{\frac{{SSW}-2\,{SW}+2\,{NW}-{NNW}}{2\,{{dy}}^{3}}}+{\frac{\left(2-\nu\right)\left(S-N-2\,{SW}+2\,{NW}+{SWW}-{NWW}\right)}{2\,{{dx}}^{2}{dy}}}=0\\\\[5.0pt]
{\frac{{EE}-2\,E+2\,W-{WW}}{2\,{{dx}}^{3}}}+{\frac{\left(2-\nu\right)\left({SE}-{SW}-2\,E+2\,W+{NE}-{NW}\right)}{2\,{{dy}}^{2}{dx}}}=0\\\\[5.0pt]
{\frac{{SEE}-2\,{SE}+2\,{SW}-{SWW}}{2\,{{dx}}^{3}}}+{\frac{\left(2-\nu\right)\left({SSE}-{SSW}-2\,{SE}+2\,{SW}+E-W\right)}{2\,{{dy}}^{2}{dx}}}=0\\\\[5.0pt]
{SE}-{NE}-{SW}+{NW}=0\\\\[5.0pt] \end{cases}$

Solving this system of boundary conditions yields the expressions for each
ghost point about corner $O$, when $O$ is the corner between $\Gamma_{N}$ and
$\Gamma_{E}$:

$\begin{cases}E=2\,O-W\\\\[5.0pt] N=2\,O-S\\\\[5.0pt]
{NE}={SE}-{SW}+{NW}\\\\[5.0pt]
EE=\frac{(-4O+4W+2SE-2SW)(\nu-2){dx}^{2}}{{dy}^{2}}+(4O-4W+WW)\\\\[5.0pt]
NN=\frac{(-4O+4S+2NW-2SW)(\nu-2){dy}^{2}}{{dx}^{2}}+(4O-4S+SS)\\\\[5.0pt]
\it{NNW}=\frac{(\nu-2)(2O-SWW-2S-2NW+2SW+NWW){dy}^{2}}{{dx}^{2}}+(2NW-2SW+SSW)\\\\[5.0pt]
\it{SEE}=\frac{(\nu-2)(2O-W-2SE+2SW+2SSE-
SSW){dy}^{2}}{{dx}^{2}}+(2SWW+2SE-2SW)\\\\[5.0pt] \end{cases}$

The same approach can be applied to the $\Gamma_{S}$ and $\Gamma_{E}$ corner,
and the results are simply a reflection of these expressions.

##  3 Non-zero Edge Conditions

In the previous section the $g_{i}(x,y)$ and $h_{i}(x,y)$ values were set to
be zero. This was done for simplicity, and clarity of exposition; in this
section we provide examples of ghost node calculations with constant edge
loading. As in the theory of boundary lifts, we can compute the solution from
a given edge loading, while setting all other boundary loads to zero. The
response to several edge loads can thus be reconstructed for the linear
problem through the principle of superposition.

In this section, a derivation of nonzero, constant values will be shown, and
the result of these loads will be determined. Again we shall only solve for
the $\Gamma_{N}$ boundary, but other boundary values will be determined via
rotation.

We first handle the case where $g_{i}(x,y)=G$ and $h_{i}(x,y)=0$. The
conditions for this case become:

$\begin{cases}{\frac{\nu\,\left(W-2\,O+E\right)}{{{dx}}^{2}}}+{\frac{N-2\,O+S}{{{dy}}^{2}}}=G\\\\[5.0pt]
{\frac{W-2\,O+E}{{{dx}}^{2}}}+{\frac{\nu\,\left(N-2\,O+S\right)}{{{dy}}^{2}}}=G\\\\[5.0pt]
{\frac{{SS}-2\,S+2\,N-{NN}}{2\,{{dy}}^{3}}}+{\frac{\left(2-\nu\right)\left({SE}-{NE}-2\,S+2\,N+{SW}-{NW}\right)}{2\,{{dx}}^{2}{dy}}}=0\\\\[5.0pt]
{\frac{{SSW}-2\,{SW}+2\,{NW}-{NNW}}{2\,{{dy}}^{3}}}+{\frac{\left(2-\nu\right)\left(S-N-2\,{SW}+2\,{NW}+{SWW}-{NWW}\right)}{2\,{{dx}}^{2}{dy}}}=0\\\\[5.0pt]
{\frac{{EE}-2\,E+2\,W-{WW}}{2\,{{dx}}^{3}}}+{\frac{\left(2-\nu\right)\left({SE}-{SW}-2\,E+2\,W+{NE}-{NW}\right)}{2\,{{dy}}^{2}{dx}}}=0\\\\[5.0pt]
{\frac{{SEE}-2\,{SE}+2\,{SW}-{SWW}}{2\,{{dx}}^{3}}}+{\frac{\left(2-\nu\right)\left({SSE}-{SSW}-2\,{SE}+2\,{SW}+E-W\right)}{2\,{{dy}}^{2}{dx}}}=0\\\\[5.0pt]
{SE}-{NE}-{SW}+{NW}=0\\\\[5.0pt] \end{cases}$

Solving this system for the conditions, as done previously, provides the
solution:

$\begin{cases}E=2\,O-W+\frac{Gdx^{2}}{\nu+1}\\\\[5.0pt]
N=2\,O-S+\frac{Gdy^{2}}{\nu+1}\\\\[5.0pt] NE=SE-SW+NW\\\\[5.0pt]
EE=\frac{(-4O+4W+2SE-2SW)(\nu-2){dx}^{2}}{{dy}^{2}}+(4O-4W+WW)-{\frac{2G{{dx}}^{4}\nu-4G{{dx}}^{4}}{{{dy}}^{2}\left(\nu+1\right)}}+{\frac{2G{{dx}}^{2}}{\nu+1}}\\\\[5.0pt]
NN=\frac{(-4O+4S+2NW-2SW)(\nu-2){dy}^{2}}{{dx}^{2}}+(4O-4S+SS)-{\frac{2G{{dy}}^{4}\nu-4G{{dy}}^{4}}{{{dx}}^{2}\left(\nu+1\right)}}+{\frac{2G{{dy}}^{2}}{\nu+1}}\\\\[5.0pt]
NNW=\frac{(\nu-2)(2O-SWW-2S-2NW+2SW+NWW){dy}^{2}}{{dx}^{2}}+(2NW-2SW+SSW)+{\frac{G\left(-2+\nu\right){{dy}}^{4}}{{{dx}}^{2}\left(\nu+1\right)}}\\\\[5.0pt]
SEE=\frac{(\nu-2)(2O-W-2SE+2SW+2SSE-
SSW){dy}^{2}}{{dx}^{2}}+(2SWW+2SE-2SW)+{\frac{G\left(-2+\nu\right){{dx}}^{4}}{{{dy}}^{2}\left(\nu+1\right)}}\\\\[5.0pt]
\end{cases}$

Similarly we can consider the case where $g_{i}(x,y)=0$ and $h_{i}(x,y)=H$.
The conditions for this system become:

$\begin{cases}{\frac{\nu\,\left(W-2\,O+E\right)}{{{dx}}^{2}}}+{\frac{N-2\,O+S}{{{dy}}^{2}}}=0\\\\[5.0pt]
{\frac{W-2\,O+E}{{{dx}}^{2}}}+{\frac{\nu\,\left(N-2\,O+S\right)}{{{dy}}^{2}}}=0\\\\[5.0pt]
{\frac{{SS}-2\,S+2\,N-{NN}}{2\,{{dy}}^{3}}}+{\frac{\left(2-\nu\right)\left({SE}-{NE}-2\,S+2\,N+{SW}-{NW}\right)}{2\,{{dx}}^{2}{dy}}}=H\\\\[5.0pt]
{\frac{{SSW}-2\,{SW}+2\,{NW}-{NNW}}{2\,{{dy}}^{3}}}+{\frac{\left(2-\nu\right)\left(S-N-2\,{SW}+2\,{NW}+{SWW}-{NWW}\right)}{2\,{{dx}}^{2}{dy}}}=H\\\\[5.0pt]
{\frac{{EE}-2\,E+2\,W-{WW}}{2\,{{dx}}^{3}}}+{\frac{\left(2-\nu\right)\left({SE}-{SW}-2\,E+2\,W+{NE}-{NW}\right)}{2\,{{dy}}^{2}{dx}}}=H\\\\[5.0pt]
{\frac{{SEE}-2\,{SE}+2\,{SW}-{SWW}}{2\,{{dx}}^{3}}}+{\frac{\left(2-\nu\right)\left({SSE}-{SSW}-2\,{SE}+2\,{SW}+E-W\right)}{2\,{{dy}}^{2}{dx}}}=H\\\\[5.0pt]
{SE}-{NE}-{SW}+{NW}=0\\\\[5.0pt] \end{cases}$

The solutions for this system are:

$\begin{cases}E=2\,O-W\\\\[5.0pt] N=2\,O-S\\\\[5.0pt]
NE={SE}-{SW}+{NW}\\\\[5.0pt]
EE=\frac{(-4O+4W+2SE-2SW)(\nu-2){dx}^{2}}{{dy}^{2}}+(4O-4W+WW)+2Hdx^{3}\\\\[5.0pt]
NN=\frac{(-4O+4S+2NW-2SW)(\nu-2){dy}^{2}}{{dx}^{2}}+(4O-4S+SS)-2Hdx^{3}\\\\[5.0pt]
NNW=\frac{(\nu-2)(2O-SWW-2S-2NW+2SW+NWW){dy}^{2}}{{dx}^{2}}+(2NW-2SW+SSW)-2Hdx^{3}\\\\[5.0pt]
SEE=\frac{(\nu-2)(2O-W-2SE+2SW+2SSE-
SSW){dy}^{2}}{{dx}^{2}}+(2SWW+2SE-2SW)+2Hdx^{3}\\\\[5.0pt] \end{cases}$

##  4 Reduction of Order and ODE Solver

To integrate in time, and produce distributed values at each time step, we
invoke MATLAB’s ode15s solver. This solver requires an abstract differential
equation of the form:

$\dot{y}=F(t,y).$ (2)

To reduce our evolutionary PDE, we introduce states to reduce the order and
create a $2\times 2$ system in the state variables to obtain a first order
formulation. The problem also has many spatial nodes, which each must be
resolved. A matrix expression will be used to obtain one principal variable
$y$ to feed to the solver. We begin with the interior PDE:

$w_{tt}+D\Delta^{2}w+a_{1}w_{x}+a_{2}w_{y}+k_{0}w_{t}-k_{1}(w_{txx}+w_{tyy})=f(x,t),$
(3)

and introduce the state variable $v$, capturing $w_{t}$, or the velocity. We
can also take the derivative of $v$ and using the PDE expression for $w_{tt}$,
we obtain the value of $v_{t}$:

$\begin{cases}v=&w_{t}\\\
v_{t}=&f(x,t)-D\Delta^{2}w-a_{1}w_{x}-a_{2}w_{y}-k_{0}w_{t}+k_{1}(w_{txx}+w_{tyy}).\end{cases}$
(4)

We can write each of these as a matrix:

$y=\begin{bmatrix}w\\\ v\end{bmatrix},y_{t}=\begin{bmatrix}w_{t}\\\
v_{t}\end{bmatrix}$

For each node in the model we must solve for $w$ and $v$, and thus we will
extend each of these in the array with indices corresponding to each point in
the mesh. There are a total of $N$ nodes, determined from user set parameters
as laid out in Section 6.

$y=\begin{bmatrix}w_{1},w_{2},\ldots,w_{N},v_{1},v_{2},\ldots,v_{N}\end{bmatrix}^{T}$

Using these vectors, this problem is now written as a first order ODE, and can
be solved directly with the ode15s solver, where each iteration will be stored
for future animation, post-processing calculations, and other data analyses.
Appendix A contains MATLAB code used to simulate the plate dynamics.

##  5 Energy Calculation

The natural potential energy associated to a free plate is expressed in terms
of the bilinear form

$a(u,v)=\nu(\Delta u,\Delta
v)+(1-\nu)\big{[}(u_{xx},v_{xx})+2(u_{xy},v_{xy})+(u_{yy},v_{yy})\big{]},$ (5)

where the expression $(f,g)$ represents the $L^{2}(\Omega)$ inner-product:

$(f,g)=(f,g)_{L^{2}(\Omega)}=\int_{\Omega}f(x,y)g(x,y)d\mathbf{x}.$

We consider the standard kinetic and potential contributions to the total
energy. This is to say, for a state $[w,v]$ we can compute the total energy as

$E(t)=E(w(t),v(t))=U(w(t))+K(v(t)).$

The potential energy is given through the bilinear form as $U(w)=a(w,w)$:

$U(w(t))=\frac{1}{2}D\int_{\Omega}\left[\nu(\Delta
w(t))^{2}+(1-\nu)\left(w^{2}_{xx}(t)+2w(t)^{2}_{xy}+w(t)^{2}_{yy}\right)\right]d\mathbf{x}$

and the kinetic energy

$K(v)=\frac{1}{2}\int_{\Omega}v(t)^{2}d\mathbf{x},$

where we have suppressed the dependence of the states on $\mathbf{x}$.

If there is no damping or flow considered in the model (i.e.,
$k_{0}=k_{1}=a_{1}=a_{2}=0$), then conservation of energy is expected for
solutions—see, for instance, [5, 10]. We demonstrate this here with an initial
state consisting of an initial displacement of $w(x,y;0)=0$ and an initial
velocity of $v(x,y;0)=x$. Figure 3 is produced by the script and demonstrates
the conservation of energy for the numerically computed solution in this
situation.

Figure 3: Conservation of energy demonstrated with a 10x10 mesh grid and
$k_{0}=k_{1}=a_{1}=a_{2}=0$, $v(x,y;0)=x$

Of course, slight variations (errors) can be observed on this small scale due
to numerical error in the finite difference method and time integrator used.

We also demonstrate decay of the energy in the presence of some damping. In a
similar case to the above, but taking a small amount of frictional damping
$k_{0}=0.1$, we observe the energy decreasing. Figure 4 is produced by the
script and demonstrates the dissipation of energy in the system.

Figure 4: Decay of energy demonstrated with a 10x10 mesh grid and
$k_{1}=a_{1}=a_{2}=0$, $k_{0}=0.1$ and $v(x,y;0)=x$

##  6 Table of Variables

In the script that was generated in correspondence to this document, there are
many variables that can be modified in order to fit the desired situation.
Table 1 contains a list of the names and descriptions of each of these
variables.

Variable | Description
---|---
D | Stiffness parameter for the cantilever material
Lx | Length of the cantilever in the x direction
Ly | Length of the cantilever in the y direction
Nx | Number of points in the mesh in the x direction
Ny | Number of points in the mesh in the y direction
N | Total number of node points where w is unknown
$\nu$ | Poisson Ratio
$k_{0}$ | Linear damping parameter
$k_{1}$ | Nonlinear damping parameter
$a_{1}$ | Flow parameter in the x direction
$a_{2}$ | flow parameter in the y direction
anim | Turns the animation on/off
energies | Turns the energy calculation on/off
$t_{0}$ | Initial time of the calculation
$t_{f}$ | Final time of the calculation
ns | Number of time steps taken
Table 1: Descriptions of each of the variables that can be chosen in the code

In addition to these variables, the user can also set the desired initial
conditions using the functions described in Table 2.

Function |  Description
---|---
winit(X, Y) |  Sets the initial position of the cantilever. This is done by defining the internal variable wmat as a function of the X and Y positional arrays.
vinit(X, Y) |  Sets the initial velocity of the cantilever. This is done by defining the internal variable vmat as a function of the X and Y positional arrays.
Table 2: Descriptions of each of the initial condition functions that can be
chosen in the code

##  7 Acknowledgements

The author acknowledges the generous support of the National Science
Foundation in this work through the grant entitled: “Collaborative Research:
Experiment, Theory, and Simulation of Aeroelastic Limit Cycle Oscillations for
Energy Harvesting Applications” (NSF-DMS 1908033).

##  Appendix A Matlab Code

⬇

1close all;

2

3%problem parameters

4D = 1.0e0; %stiffness parameter

5Lx = 1; %length in x direction

6Ly = 1; %length in y direction

7Nx = 10; %number of nodes in x direction

8Ny = 10; %number of nodes in y direction

9N = (Nx-1)*(Ny); %total number of nodes where w is unknown

10nu=0.3;

11k0 = 0.1; %damping parameter

12k1 = 0; %damping parameter

13a1 = 0;

14a2 = 0;

15anim = 1; %see animation

16energies = 1; %compute and plot energies

17

18%timestepping parameters

19t0 = 0; %initial time

20tf = 10;

210e0; %final time

22ns = ceil(50*tf)+1; %number of time steps

23t = linspace(t0, tf, ns); %time vector

24

25

26%initialize spatial domain

27x = linspace(0,Lx,Nx);

28y = linspace(0,Ly,Ny);

29[X,Y] = meshgrid(x,y);

30dx = x(2)-x(1);

31dy = y(2)-y(1);

32

33%initialize solution array W = [wn; vn];

34%get initial diplacement and velocity - form W0

35wn = winit(X,Y);

36vn = vinit(X,Y);

37W0 = [wn; vn]; %initial data vector is of length 2*N

38

39

40% dydt = RHS(t,W0,D,Nx,Ny,N,x,y,dx,dy,nu);

41% dydt = RHS(t,dydt,D,Nx,Ny,N,x,y,dx,dy,nu);

42% dydt = RHS(t,dydt,D,Nx,Ny,N,x,y,dx,dy,nu);

43% dydt = RHS(t,dydt,D,Nx,Ny,N,x,y,dx,dy,nu);

44

45% call the timestepping code

46tic;

47[T,w] = ode15s(@(t,W) RHS(t,W,D,Nx,Ny,N,x,y,dx,dy,nu,k0,k1,a1,a2),t,W0);

48odetime = toc;

49fprintf(’Total run time for ODE integrator: %g\n’,odetime);

50

51%postprocess the results

52%plot the surface

53wmat = [zeros(Ny,1) reshape(w(end,1:N),[Nx-1 Ny])’];

54

55% figure(3); surf(X,Y,wmat); %plot the final surface

56% zlabel({’w(x,y)’});

57% ylabel({’y’});

58% xlabel({’x’});

59% title(’Solution at final time’);

60

61%animate the solution

62wvals=w(:,1:N);

63wmax=max(max(abs(wvals)));

64if anim>0

65 fig = figure(11);%hold on

66 F(ns) = struct(’cdata’,[],’colormap’,[]);

67 grid on;

68 plt=surf(X,Y,[zeros(Ny,1) reshape(w(1,1:N),[Nx-1 Ny])’ ]);

69 colormap(jet(256));

70 colorbar;

71 caxis([-1.1*wmax,1.1*wmax]);

72 zlim([-1.1*wmax,1.1*wmax])

73 view(-40,30);

74 str=sprintf(’t = %3.3f’,T(1));

75 h=text(0.95,0.95,0.9*wmax,str,’FontSize’,12);

76 F(1) = getframe(fig);

77 for j = 2:ns

78 plt.ZData =[zeros(Ny,1) reshape(w(j,1:N),[Nx-1 Ny])’];

79 str=sprintf(’t = %3.3f’,T(j));

80 set(h,’String’,str);

81 drawnow % display updates

82 F(j) = getframe(fig);

83 end

84 movie(fig,F,2,anim*ns/tf);

85

86 vid = VideoWriter(’animation.mp4’,’MPEG-4’);

87 open(vid);

88 writeVideo(vid,F);

89 close(vid);

90end

91

92% compute energies if desired

93if energies==1

94 tic;

95 %compute and plot energies

96 Le = zeros(size(ns));

97 Nle = zeros(size(ns));

98 for j=1:ns

99 Le(j)= computeEnergies(w(j,:),Nx,Ny,x,y,dx,dy,D,nu);

100 end

101 energytime = toc;

102 fprintf(’Total run time for energy calculation: %g\n’,energytime);

103 figure;hold

104 plot(T,Le,’-b’, ’DisplayName’,’E(t)’);

105 legend(’linear E(t)’)

106end

107

108%end

109

110

111

112%function for initial displacement

113function w = winit(X,Y)

114

115Lx = X(1,end); %get Lx

116Ly = Y(end,1); %get Ly

117[Ny,Nx] = size(X); %get Nx and Ny

118%wmat = (-4*X.^5+15*X.^4-20*X.^3+10*X.^2); %define the initial value of w

119wmat = 0.0*X; %define the initial value of w

120% figure; surf(X,Y,wmat); %plot the initial surface

121% zlabel({’w(x,y)’});

122% ylabel({’y’});

123% xlabel({’x’});

124% title(’Initial Displacement’);

125wmat = wmat(:,2:Nx); %trim off left values

126w = reshape(wmat’,[],1); %reshape matrix into the w vector

127end

128

129%function for initial velocity

130function v = vinit(X,Y)

131Lx = X(1,end); %get Lx

132Ly = Y(end,1); %get Ly

133[Ny,Nx] = size(X); %get Nx and Ny

134vmat = X; %define the initial value of w

135% figure; surf(X,Y,vmat); %plot the initial surface

136% % Create zlabel

137% zlabel({’v(x,y)’});

138% % Create ylabel

139% ylabel({’y’});

140% % Create xlabel

141% xlabel({’x’});

142% title(’Initial Velocity’);

143vmat = vmat(:,2:Nx); %trim off left

144v = reshape(vmat’,[],1); %reshape matrix into the w vector

145end

146

147%function for computing energies

148function [le,nle] = computeEnergies(w,Nx,Ny,x,y,dx,dy,D,nu)

149

150%separate the displacements from velocities

151w = reshape(w,[Nx-1 2*(Ny)])’;

152wmat = [zeros(Ny,1),w(1:Ny,:)]; %displacements at interior nodes

153vmat = [zeros(Ny,1),w(Ny+1:2*Ny,:)]; %velocities at interior nodes

154

155%compute the appropriate derivatives at the interior points

156[lx,ly,lxy] = laplace(wmat,Nx,Ny,dx,dy,nu);

157%compute the integral of grad w squared

158integrand = nu * (lx + ly).^2 + (1-nu)*(lx.^2 + 2.*lxy.^2 + ly.^2);

159normpotential2 = trapz(y,trapz(x,integrand,2));

160

161%compute the integral of the velocities squared

162normvel2 = trapz(y,trapz(x,vmat.*vmat,2));

163

164%compute linear energy

165le = 0.5*D*normpotential2 + 0.5*normvel2;

166

167end

168function [lx,ly,lxy] = laplace(w,Nx,Ny,dx,dy,nu)

169

170%initialize the output arrays

171lx = zeros(size(w));

172ly = zeros(size(w));

173lxy = zeros(size(w));

174%apply BCs to wmat so we can compute w_xx and w_yy

175%build the array with ghost nodes

176w = [zeros(1,size(w,2)+2); [zeros(size(w,1),1), w, zeros(size(w,1),1)];
zeros(1,size(w,2)+2)];

177

178%Clamped edge ghost points

179w(:,1) = w(:,3);

180

181%free edge (avoiding free-free corner area)

182for j=2:Nx

183 i = Ny+1;

184 w(i+1,j) = -nu * (w(i,j-1) - 2 * w(i,j) + w(i,j+1)) / dx ^ 2 * dy ^ 2 -
w(i-1,j) + 2 * w(i,j);

185end

186

187for j=2:Nx

188 i = 2;

189 w(i-1,j) = -nu * (w(i,j-1) - 2 * w(i,j) + w(i,j+1)) / dx ^ 2 * dy ^ 2 -
w(i+1,j) + 2 * w(i,j);

190end

191

192for i=3:Ny+1

193 j = Nx+1;

194 w(i,j+1) = -nu * (w(i-1,j) - 2 * w(i,j) + w(i+1,j)) / dy ^ 2 * dx ^ 2 -
w(i,j-1) + 2 * w(i,j);

195end

196

197

198%Now we will apply the conditions for the corners.

199%Labelling convention in this section is relative to the corner

200

201%Top right corner

202i = 2;

203j = Nx + 1;

204%E

205w(i,j+1) = 2 * w(i,j) - w(i,j-1);

206

207%N

208w(i-1,j) = 2 * w(i,j) - w(i+1,j);

209

210%NE

211w(i-1,j+1) = w(i+1,j+1) - w(i+1,j-1) + w(i-1,j-1);

212

213

214%Bottom right corner

215i = Ny + 1;

216j = Nx + 1;

217%E

218w(i,j+1) = 2 * w(i,j) - w(i,j-1);

219

220

221%S

222w(i+1,j) = 2 * w(i,j) - w(i-1,j);

223

224%SE

225w(i+1,j+1) = w(i-1,j+1) - w(i-1,j-1) + w(i+1,j-1);

226

227%loop over interior nodes and compute w_x and w_y using 2nd order centered

228%differences

229for j = 2:Nx+1 %which row of nodes

230 for i=2:Ny+1 %which column of nodes

231 lx(i-1,j-1) = (1/(dx^2))*(w(i,j+1)-2*w(i,j)+w(i,j-1)); %w_xx calculation

232 ly(i-1,j-1) = (1/(dy^2))*(w(i+1,j)-2*w(i,j)+w(i-1,j)); %w_yy calculation

233 lxy(i-1,j-1) = (1/(4*dx*dy)) * ( w(i+1,j+1) - w(i-1,j+1) - w(i+1,j-1) +
w(i-1,j-1));

234 end

235end

236end

237

238function [gx,gy] = grads(w,Nx,Ny,dx,dy)

239

240gx = zeros(size(w));

241gy = zeros(size(w));

242

243%apply BCs to wmat so we can compute w_xx and w_yy

244%build the array with ghost nodes

245w = [zeros(1,size(w,2)+2); [zeros(size(w,1),1), w, zeros(size(w,1),1)];
zeros(1,size(w,2)+2)];

246%Applying hinged-hinged conditions

247w(1,:) = -1.0 * w(3,:);

248w(end,:) = -1.0 * w(end-2,:);

249

250w(:,1) = -1.0 * w(:,3);

251w(:,end) = -1.0 * w(:,end-2);

252%loop over interior nodes and compute w_x and w_y using 2nd order centered

253%differences

254for j = 2:Ny+1 %which row of nodes

255 for i=2:Nx-1 %which column of nodes

256 gx(j-1,i) = (1/(2*dx))*(w(j,i+1)-w(j,i-1)); %w_x calculation

257 gy(j-1,i) = (1/(2*dy))*(w(j+1,i)-w(j-1,i)); %w_y calculation

258 end

259end

260end

261

262%define the forcing function f(x,y,t)

263function z = f(~,~,~)

264z=0;

265end

266

267%define the ODE RHS function

268function dydt=RHS(t,W,D,Nx,Ny,N,x,y,dx,dy,nu,k0,k1,a1,a2)

269fprintf(’Time: %g\n’,t);

270%reshape W into a matrix with 2*Ny rows and Nx-2 columns

271wmat = reshape(W,[Nx-1 2*(Ny)])’;

272w = wmat(1:Ny,:); %displacements at interior nodes

273v = wmat(Ny+1:2*Ny,:); %velocities at interior nodes

274

275%initialize dydt

276dydt = zeros(size([w;v]));

277boundaryWxx = [ones(size(w,1),1) ones(size(w))];

278boundaryWxxx = [ones(size(w,1),1) ones(size(w))];

279

280dydt(1:Ny,1:Nx-1) = v; %place v values in the w_t locations

281

282%pad v vector with zeros on all sides

283v = [zeros(1,size(v,2)+2); [zeros(size(v,1),1), v, zeros(size(v,1),1)];
zeros(1,size(v,2)+2)];

284v = [zeros(1,size(v,2)+2); [zeros(size(v,1),1), v, zeros(size(v,1),1)];
zeros(1,size(v,2)+2)];

285

286

287gradarr=[zeros(size(w,1),1), w];

288[gradx,grady] = grads(gradarr,Nx,Ny,dx,dy);

289normgrad2 = trapz(y,trapz(x,gradx.*gradx + grady.*grady,2));

290

291%pad w vector with zeros on all sides to calculate the grad w squared

292w = [zeros(1,size(w,2)+2); [zeros(size(w,1),1), w, zeros(size(w,1),1)];
zeros(1,size(w,2)+2)];

293

294w = [zeros(1,size(w,2)+2); [zeros(size(w,1),1), w, zeros(size(w,1),1)];
zeros(1,size(w,2)+2)];

295

296%Clamped edge ghost points

297w(:,1) = w(:,3);

298

299%free edge (avoiding free-free corner area)

300for j=2:Nx

301 i = Ny+2;

302 w(i+1,j) = -nu * (w(i,j-1) - 2 * w(i,j) + w(i,j+1)) / dx ^ 2 * dy ^ 2 -
w(i-1,j) + 2 * w(i,j);

303end

304

305for j=2:Nx

306 i = 3;

307 w(i-1,j) = -nu * (w(i,j-1) - 2 * w(i,j) + w(i,j+1)) / dx ^ 2 * dy ^ 2 -
w(i+1,j) + 2 * w(i,j);

308end

309

310for i=4:Ny+1

311 j = Nx+1;

312 w(i,j+1) = -nu * (w(i-1,j) - 2 * w(i,j) + w(i+1,j)) / dy ^ 2 * dx ^ 2 -
w(i,j-1) + 2 * w(i,j);

313end

314

315%free edge, 2nd row (avoiding free-free corner area)

316

317for j=3:Nx-1

318 i = Ny+2;

319 w(i + 2, j) = -(2 - nu) / dx ^ 2 * dy ^ 2 * (w(i + 1, j + 1) - w(i - 1, j
+ 1) - 2 * w(i + 1, j) + 2 * w(i - 1, j) + w(i + 1, j - 1) - w(i - 1, j - 1))
+ 2 * w(i + 1, j) - 2 * w(i - 1, j) + w(i - 2, j);

320end

321

322for j=3:Nx-1

323 i = 3;

324 w(i - 2, j) = (2 - nu) / dx ^ 2 * dy ^ 2 * (w(i + 1, j + 1) - w(i - 1, j +
1) - 2 * w(i + 1, j) + 2 * w(i - 1, j) + w(i + 1, j - 1) - w(i - 1, j - 1)) +
w(i + 2, j) - 2 * w(i + 1, j) + 2 * w(i - 1, j);

325end

326

327for i=5:Ny

328 j = Nx+1;

329 w(i,j+2) = -(2 - nu) * (w(i + 1, j + 1) - w(i + 1, j - 1) - 2 * w(i, j +
1) + 2 * w(i, j - 1) + w(i - 1, j + 1) - w(i - 1, j - 1)) / dy ^ 2 * dx ^ 2 +
2 * w(i, j + 1) - 2 * w(i, j - 1) + w(i, j - 2);

330end

331

332%Now we will apply the conditions for the corners.

333%Labelling convention in this section is relative to the corner

334

335%Top right corner

336i = 3;

337j = Nx + 1;

338%E

339w(i,j+1) = 2 * w(i,j) - w(i,j-1);

340

341%EE

342w(i, j + 2) = (-4 * (nu - 2) * (w(i, j) - w(i, j - 1) - w(i + 1, j + 1) / 2
+ w(i + 1, j - 1) / 2) * dx ^ 2 + 4 * (w(i, j) - w(i, j - 1) + w(i, j - 2) /
4) * dy ^ 2) / dy ^ 2;

343

344%N

345w(i-1,j) = 2 * w(i,j) - w(i+1,j);

346

347%NE

348w(i-1,j+1) = w(i+1,j+1) - w(i+1,j-1) + w(i-1,j-1);

349

350%NN

351w(i-2, j) = (-4 * (nu - 2) * (w(i, j) - w(i + 1, j) - w(i - 1, j - 1) / 2 +
w(i + 1, j - 1) / 2) * dy ^ 2 + 4 * dx ^ 2 * (w(i, j) - w(i + 1, j) + w(i + 2,
j) / 4)) / dx ^ 2;

352

353%NNW

354w(i-2, j-1) = (2 * (w(i, j) - w(i + 1, j - 2) / 2 - w(i + 1, j) - w(i - 1,
j - 1) + w(i + 1, j - 1) + w(i - 1, j - 2) / 2) * (nu - 2) * dy ^ 2 + 2 * dx ^
2 * (w(i - 1, j - 1) - w(i + 1, j - 1) + w(i + 2, j - 1) / 2)) / dx ^ 2;

355

356%SEE

357w(i+1, j+2) = (2 * (w(i, j) - w(i, j - 1) - w(i + 1, j + 1) + w(i + 1, j -
1) + w(i + 2, j + 1) / 2 - w(i + 2, j - 1) / 2) * (nu - 2) * dx ^ 2 + dy ^ 2 *
(w(i + 1, j - 2) + 2 * w(i + 1, j + 1) - 2 * w(i + 1, j - 1))) / dy ^ 2;

358

359

360%Bottom right corner

361i = Ny + 2;

362j = Nx + 1;

363%E

364w(i,j+1) = 2 * w(i,j) - w(i,j-1);

365

366%EE

367w(i, j + 2) = (-4 * (nu - 2) * (w(i, j) - w(i, j - 1) - w(i - 1, j + 1) / 2
+ w(i - 1, j - 1) / 2) * dx ^ 2 + 4 * (w(i, j) - w(i, j - 1) + w(i, j - 2) /
4) * dy ^ 2) / dy ^ 2;

368

369%S

370w(i+1,j) = 2 * w(i,j) - w(i-1,j);

371

372%SE

373w(i+1,j+1) = w(i-1,j+1) - w(i-1,j-1) + w(i+1,j-1);

374

375%SS

376w(i+2, j) = (-4 * (nu - 2) * (w(i, j) - w(i - 1, j) - w(i + 1, j - 1) / 2 +
w(i - 1, j - 1) / 2) * dy ^ 2 + 4 * dx ^ 2 * (w(i, j) - w(i - 1, j) + w(i- 2,
j) / 4)) / dx ^ 2;

377

378%SSW

379w(i+2, j-1) = (2 * (w(i, j) - w(i - 1, j - 2) / 2 - w(i - 1, j) - w(i + 1,
j - 1) + w(i - 1, j - 1) + w(i + 1, j - 2) / 2) * (nu - 2) * dy ^ 2 + 2 * dx ^
2 * (w(i + 1, j - 1) - w(i - 1, j - 1) + w(i - 2, j - 1) / 2)) / dx ^ 2;

380

381%NEE

382w(i-1, j+2) = (2 * (w(i, j) - w(i, j - 1) - w(i - 1, j + 1) + w(i - 1, j -
1) + w(i - 2, j + 1) / 2 - w(i - 2, j - 1) / 2) * (nu - 2) * dx ^ 2 + dy ^ 2 *
(w(i - 1, j - 2) + 2 * w(i - 1, j + 1) - 2 * w(i - 1, j - 1))) / dy ^ 2;

383

384

385%compute the RHS of the velocity equation

386for j = 3:Nx+1 %which column of nodes

387 for i=3:Ny+2 %which row of nodes

388

389 %In vacuo conditions

390 wxxxx = (1/dx^4) * (w(i-2,j) - 4*w(i-1,j) + 6*w(i,j)...

391 - 4 *w(i+1,j) + w(i+2,j));

392

393 wyyyy = (1/dy^4) * (w(i,j-2) - 4*w(i,j-1) + 6*w(i,j)...

394 - 4 *w(i,j+1) + w(i,j+2));

395

396 wxxyy = (1/(dx^2*dy^2)) * ( w(i+1,j+1) - 2* w(i+1,j)...

397 + w(i+1,j-1) - 2* w(i,j+1) + 4 * w(i,j) -2*w(i,j-1)...

398 +w(i-1,j+1) -2 * w(i-1,j)+w(i-1,j-1));

399

400 %damping term

401 wt = v(i,j);

402

403 wtxx = (1/dx^2) * (v(i+1,j) - 2* v(i,j) + v(i-1,j));

404

405 wtyy = (1/dy^2) * (v(i,j+1) - 2* v(i,j) + v(i,j-1));

406

407 damping = - k0*(wt) + k1*(wtxx + wtyy);

408

409 %spatial

410 wx = (1/(2*dx))*(w(i,j+1)-w(i,j-1)); %w_x calculation

411 wy = (1/(2*dy))*(w(i+1,j)-w(i-1,j)); %w_y calculation

412

413 spatial = - a1 * wx - a2 * wy;

414

415 dydt(Ny+(i-2),j-2) = - 1 * D *(wxxxx + 2 * wxxyy + wyyyy)...

416 + f(x(j-2),y(i-2),t)...

417 + spatial + damping;

418 end

419end

420%dydt

421dydt = reshape(dydt’,[2*N 1]);

422

423end

## References

  * [1] Dieter Bahlmann and Ulrich Langer  “A fast solver for the first biharmonic boundary value problem”  In _Numerische Mathematik_ 63.1 Springer, 1992, pp. 297–313
  * [2] MV Barton  “Finite difference equations for the analysis of thin rectangular plates with combinations of fixed and free edges”, 1948
  * [3] Heribert Blum, Rolf Rannacher and Rolf Leis  “On the boundary value problem of the biharmonic operator on domains with angular corners”  In _Mathematical Methods in the Applied Sciences_ 2.4 Wiley Online Library, 1980, pp. 556–581
  * [4] Fioralba Cakoni, George C Hsiao and Wolfgang L Wendland  “On the boundary integral equation method for a mixed boundary value problem of the biharmonic equation”  In _Complex Variables, Theory and Application: An International Journal_ 50.7-11 Taylor & Francis, 2005, pp. 681–696
  * [5] Maria Deliyianni and Justin T Webster  “Theory of Solutions for An Inextensible Cantilever”  In _Applied Mathematics & Optimization_ Springer, 2021, pp. 1–55
  * [6] Maria Deliyianni, Kevin McHugh, Justin T Webster and Earl Dowell  “Dynamic equations of motion for inextensible beams and plates”  In _arXiv preprint arXiv:2107.07454_ , 2021
  * [7] Maria Deliyianni, Varun Gudibanda, Jason Howell and Justin T Webster  “Large deflections of inextensible cantilevers: modeling, theory, and simulation”  In _Mathematical Modelling of Natural Phenomena_ 15 EDP Sciences, 2020, pp. 44
  * [8] Earl Dowell and Kevin McHugh  “Equations of motion for an inextensible beam undergoing large deflections”  In _Journal of Applied Mechanics_ 83.5 American Society of Mechanical Engineers, 2016, pp. 051007
  * [9] Louis W Ehrlich  “Solving the biharmonic equation as coupled finite difference equations”  In _SIAM Journal on Numerical Analysis_ 8.2 SIAM, 1971, pp. 278–287
  * [10] Jason Howell, Katelynn Huneycutt, Justin T Webster and Spencer Wilder  “A thorough look at the (in) stability of piston-theoretic beams”  In _arXiv preprint arXiv:1903.07244_ , 2019
  * [11] Jason S Howell, Daniel Toundykov and Justin T Webster  “A cantilevered extensible beam in axial flow: semigroup well-posedness and postflutter regimes”  In _SIAM Journal on Mathematical Analysis_ 50.2 SIAM, 2018, pp. 2048–2085
  * [12] PCM Lau  “Curvilinear finite difference method for biharmonic equation”  In _International Journal for Numerical Methods in Engineering_ 14.6 Wiley Online Library, 1979, pp. 791–812
  * [13] DM Tang and EH Dowell  “Aeroelastic response and energy harvesting from a cantilevered piezoelectric laminated plate”  In _Journal of Fluids and Structures_ 76 Elsevier, 2018, pp. 14–36
  * [14] DM Tang, Dani Levin and EH Dowell  “Experimental and theoretical correlations for energy harvesting from a large flapping flag response”  In _Journal of Fluids and structures_ 86 Elsevier, 2019, pp. 290–315


The logic $\LL^1_\kappa$ was introduced by Shelah in [3]. In [4], he proved that for a strongly compact cardinal $\kappa$, it admits the following algebraic characterization: two structures are $\LL^1_\kappa$-equivalent if and only if they have isomorphic iterated ultrapowers via $\kappa$-complete ultrafilters. We give  presentation of the logic $\LL^1_\kappa$ and a simplified and slightly modified proof of this result.

§ THE LOGIC $\LL^1_{\KAPPA}$

The logic $\LL^1_{\kappa}$ is defined through a variation of an Ehrenfeucht-Fraïssé game. The states of this game will be triples $(\alpha,f,\pi)$, where $\alpha$ is an ordinal, $\pi$ is a partial isomorphism, and $f$ is a partition function which partitions some subset of the field of $\pi$ into countably many pieces.

[The game $\G^\beta_\theta$] Let $\AA$ and $\BB$ be structures of same signature, let $\beta$ be an ordinal and let $\theta$ be a cardinal. The game
\[
\G^\beta_\theta(\AA,\BB)
\]

is played as follows.

Starting state: The starting state is $(\beta,\emptyset,\emptyset)$.

Further states: Assume that the game is at state $(\alpha, f,\pi)$.

* The player $\1$ chooses some ordinal $\alpha'<\alpha$ and some set $X\in\AA^{\leq\theta}\cup\BB^{\leq\theta}$.

* The player $\2$ chooses a partial partition function $f':\AA\cup\BB\to\omega$ such that ${\dom(f),X\subseteq \dom(f')}$
and such that for all $a\in \dom(f)$,
\[
\]

Then she chooses a partial isomorphism $\pi'\supseteq\pi$ such that
\[{f'^{-1}\{0\}\subseteq \fld(\pi')}.\]

The next state is $(\alpha',f',\pi')$.

The player to first break the rules loses.

Let $\equiv^\beta_\theta$ be the transitive closure of the relation

\[
\text{The player }\2\text{ has a winning strategy in the game }\G^\beta_\theta(\AA,\BB).
\]

A logic is a class function associating to each signature $\tau$ a collection of sentences and a satisfaction relation, satisfying certain regularity properties, see [1].

[The logic $\LL^1_\kappa$] Let $\tau$ be a signature.

* A $\tau$-sentence in $\LL^1_\kappa$ is a class of $\tau_0$-structures which is closed under the relation $\equiv^\beta_\theta$, for some $\tau_0\in[\tau]^{<\kappa}$ and some $\beta,\theta<\kappa$.

* The satisfaction relation of $\LL^1_\kappa$ is defined as
\[
\MM\models\phi\quad:\iff\quad\MM\rest\tau_0\in\phi,
\]

where $\tau_0$ is the signature such that $\phi$ consists of $\tau_0$-structures.

For cardinals of the form $\kappa=\beth_\kappa$, the logic $\LL^1_\kappa$ is a regular logic strictly between the logics $\LL_{\kappa\omega}$ and $\LL_{\kappa\kappa}$.

See [3].

Notice that for any $\tau$-structures $\AA$ and $\BB$,
\[
\AA\equiv_{\LL^1_\kappa}\BB\quad\iff\quad\forall \tau_0\in[\tau]^{<\kappa}\quad\forall\beta,\theta<\kappa\quad \AA\rest\tau_0\equiv^\beta_\theta\BB\rest\tau_0.
\]

Assume that $\kappa=\beth_\kappa$. Assume that $\bar{\AA}=\{\AA_n\}_{n\in \omega}$ is an $\LL_{\kappa\kappa}$-elementary chain of structures. Then, for each $n$,

\[
\AA_n\equiv_{\LL^1_{\kappa}}\bigcup\bar{\AA}.
\]

See [3].

§ ULTRAPOWERS AND $\LL^1_\KAPPA$-THEORIES

If $\UU$ is an ultrafilter on a set $I$ and $\AA$ is a structure, the ultrapower $\AA^{I}/\UU$ will be denoted by $\AA^\UU$. If $\UUU=(\UU_n)_n$ is a sequence of ultrafilters on some sets and $\AA$ is a structure, the iterated ultrapower of $\AA$ along the ultrafilters $(\UU_n)_n$ will be denoted by $\AA^{\UUU}$. In other words, $\AA^{\UUU}$ is the direct limit of the system

\[
\]

\begin{align*}
\end{align*}

and the maps $j_{m,n}:\AA_m\to\AA_n$ are compositions of the ultrapower embeddings. In case the ultrafilters are $\kappa$-complete, we have:

[Łos] If $\UU$ is a $\kappa$-complete ultrafilter on a set $I$ and $\AA$ is a structure, then the ultrapower embedding
\[
\AA\to\AA^\UU,\quad a\mapsto[(a)_{i\in I}]_\UU
\]

is $\LL_{\kappa\kappa}$-elementary.

By the Łos Theorem, thus, if the ultrafilters $\UU_n$ are $\kappa$-complete, then the maps $j_{m,n}:\AA_m\to\AA_n$ are $\LL_{\kappa\kappa}$-elementary.

The direct limit $\AA^{\UUU}$ comes together with embeddings
\[
\]

which are first-order elementary but not more in general. In particular, the limit embeddings might fail to be $\LL_{\kappa\kappa}$-elementary, even if the ultrafilters were $\kappa$-complete. In this case, they might even fail to be $\LL^1_\kappa$-elementary, but by the Proposition <ref>, they still preserve the $\LL^1_\kappa$-theory.

For instance, any ultrapower of a well-founded model by a $\kappa$-complete ultrafilter is again well-founded, since well-foundedness is expressible in the logic $\LL_{\omega_1\omega_1}$, and thus preserved under $\LL_{\kappa\kappa}$-elementary embeddings (in case $\kappa$ is uncountable). However, it is easy to produce an ill-founded model from a well-founded one by iterating the ultrapower construction $\omega$ many times, as will be done in the proof of characterization of $\LL^1_\kappa$.

§.§ Strongly compact cardinals
For cardinals $\lambda\geq\kappa$, denote
\[
\Pow_\kappa(\lambda):=\{x\subseteq\lambda:|x|<\kappa\}.
\]

An ultrafilter $\UU$ on $\Pow_\kappa(\lambda)$ is fine if it is $\kappa$-complete and for each $x\in\Pow_\kappa(\lambda)$, it contains the cone

\[
C_x:=\{y\in\Pow_\kappa(\lambda):x\subseteq y\}.
\]

A cardinal $\kappa$ is $\lambda$-compact if there exists a fine ultrafilter on $\Pow_\kappa(\lambda)$. A cardinal $\kappa$ is strongly compact if it is $\lambda$-compact for every $\lambda\geq\kappa$. The $\lambda$-compact cardinals have the following covering property:

Assume that $\kappa$ is a $\lambda$-compact cardinal and $\UU$ is a fine ultrafilter on $\Pow_\kappa(\lambda)$. Assume that $(H,\in)$ is a transitive model of $\ZFC^-$ closed under $<\kappa$-sequences such that $\kappa,\lambda\in H$. For any set $Y\subseteq H^\UU$ of size at most $\lambda$, there is a set $X\in H^\UU$ such that

\[
Y\subseteq X\quad\text{and}\quad H^\UU\models|X|< j(\kappa),
\]

where $j:H\to H^\UU$ is the ultrapower embedding.

Let $Y\subseteq H^\UU$ be a set of size at most $\lambda$. We find a set $X\in H^\UU$ which covers $Y$ and for which
\[
\]

Say $Y=\{[f_i]_\UU:i<\lambda\}$. Define the function $F:\Pow_\kappa(\lambda)\to H$,
\[
F(x)=\{f_i(x):i\in x\}.
\]

As $H$ is closed under $<\kappa$-sequences, this function $F$ has indeed its range inside $H$, thus $[F]_\UU\in H^\UU$.  Let $X:=[F]_\UU$.

By fineness we have $Y\subseteq X$: for each $i<\lambda$,

\[
{C_{\{i\}}\subseteq\{x:f_i(x)\in F(x)\}\in\UU}.
\]

Also ${H^\UU\models|X|< e(\kappa)}$: simply because

\[
\{x:|F(x)|<\kappa\}=\Pow_\kappa(\lambda)\in\UU.
\]

§ PROOF OF THE CHARACTERIZATION

We now give a proof of the following theorem.

[Shelah, Theorem 1.5 in [4]] Assume that $\kappa$ is a strongly compact cardinal. The following are equivalent:

* $\AA\equiv_{\LL^1_\kappa}\BB$.

* There is a sequence $\bar{\UU}=(\UU_n)_{n<\omega}$ of $\kappa$-complete ultrafilters such that
\[
\AA^{\bar{\UU}}\cong\BB^{\bar{\UU}}.
\]

(1)$\To$(2): Assume that $\AA\equiv_{\LL^1_\kappa}\BB$. For simplicity, assume that the signature $\tau$ of the models $\AA$ and $\BB$ is relational and of size $<\kappa$, and the domains of $\AA$ and $\BB$ are disjoint. For simplicity again, assume that for all $\beta,\theta<\kappa$, the player $\2$ has a winning strategy in the game
\[
\G^\beta_\theta(\AA,\BB).
\]

We will build a countable sequence of ultrafilters $\UUU$ such that the iterated ultrapowers $\AA^{\UUU}$ and $\BB^{\UUU}$ are isomorphic.

Let $\mu$ be a regular cardinal large enough such that the models $\AA$ and $\BB$, $\kappa$, and all the winning strategies are in $H(\mu)$. For all $\beta,\theta<\kappa$, fix some winning strategy $\sigma_{\beta,\theta}$ for the player $\2$ in the game $\G^\beta_\theta(\AA,\BB)$. Choose new unary predicate symbols $A$ and $B$ and a new binary function symbol $\sigma$. Define the structure

\[
\HH:=(H(\mu),\in,A^\HH,B^\HH,\sigma^\HH,R^\HH)_{R\in\tau}
\]

* $A^\HH=\dom(\AA)$

* $B^\HH=\dom(\BB)$

* $\sigma^\HH(\beta,\theta)=\begin{cases}
\sigma_{\beta,\theta},\quad&\text{if }\beta,\theta\in\kappa\\
\emptyset &\text{otherwise.}
\end{cases}$

* For each symbol $R\in\tau$, $R^\HH=R^{\AA}\cup R^\BB$.

We will now build structures $(\HH_n)_n$, $(\AA_n)_n$, $(\BB_n)_n$, ultrafilters $(\UU_n)_n$ and sets $(X_n)_n$, by recursion on $\omega$.

Step $0$: Let $\HH_0:=\HH$, $\AA_0:=\AA$ and $\BB_0:=\BB$.

Step $n+1$:
Assume that $\HH_m$, $\AA_m$ and $\BB_m$ have been defined for all $m\leq n$. For each $m\leq n$, denote
\[
\lambda_m:=|\AA_m|+|\BB_m|+\kappa.
\]

Furthermore, assume that for all $m<n$, we have defined (using the fact that $\kappa$ is strongly compact)

* A fine ultrafilter $\UU_m$ on the set $\Pow_\kappa(\lambda_m)$.

* Its corresponding ultrapower embedding
\[
\]

* A set $X_m\in\HH_{m+1}$ such that the pointwise images $e_m[\AA_m]$ and $e_m[\BB_m]$ are subsets of $X_m$ and
\[
\HH_{m+1}\models|X_m|<e_m(\kappa),
\]
using the covering property of compact cardinals as in Lemma <ref>.

We now define the ultrafilter $\UU_n$, the model $\HH_{n+1}$, an embedding $e_n$, the set $ X_n $, and the models $ \AA_{n+1}$ and $ \BB_{n+1}$.

* Let $\UU_{n}$ be any fine ultrafilter on $\Pow_\kappa(\lambda_n)$. This is possible because $\kappa$ is strongly compact.

* Let $\HH_{n+1}:=\HH_n^{\UU_n}.$

* Let $e_n:\HH_n\to\HH_{n+1}$
be the ultrapower embedding. Notice that this embedding is $\LL_{\kappa\kappa}$-elementary and its critical point is $\kappa$.

* Let $X_n\in\HH_{n+1}$ be a set such that

\[
e_n[\AA_n],e_n[\BB_n]\subseteq X_n\quad\text{and}\quad \HH_{n+1}\models |X_n|<e_n(\kappa).
\]

This is possible by the covering properties of $\lambda_n$-compact cardinals, by Lemma <ref>.

* Finally, let

\begin{align*}
& \AA_{n+1}:=\AA_n^{\UU_n}\\
\end{align*}

We have the directed system

\[
\left(\HH_n,e_{m,n}\right)_{m<n<\omega},
\]

where each $e_{m,n}:\HH_m\to\HH_n$ is an $\LL_{\kappa\kappa}$-elementary embedding, obtained by composing the ultrapower embeddings. Let $\HH^{\UUU}$ be the direct limit of this system.

The restricted maps

\begin{align*}
\end{align*}

are also $\LL_{\kappa\kappa}$-elementary.
We get the directed systems

\[
\left(\AA_n,e^\AA_{m,n}\right)_{m<n<\omega}\quad\text{and}\quad\left(\BB_n,e^\BB_{m,n}\right)_{m<n<\omega},
\]

and we can take the direct limits of these systems, denote them by $\AA^{\UUU}$ and $\BB^{\UUU}$, respectively.

We have the first-order elementary limit embeddings:

\begin{align*}
\end{align*}

The models $\AA^{\UUU}$ and $\BB^{\UUU}$ are isomorphic.

Notice first that for each $n$, the $n$th iterates $\AA_n$ and $\BB_n$ are isomorphic to the structures $A^{\HH_n}$ and $B^{\HH_n}$,
respectively. Thus also
\[
\AA^{\UUU}\cong A^{\HH^{\UUU}}\quad\text{and}\quad \BB^{\UUU}\cong B^{\HH^{\UUU}}.
\]

It is thus enough to show that $A^{\HH^{\UUU}}$ and $B^{\HH^{\UUU}}$ are isomorphic.

By the first-order elementarity of the map $e_{0,\omega}$,

\begin{align*}
\HH^{\UUU}\models\quad \q\forall\beta,\theta<e_{0,\omega}(\kappa)\quad&\sigma^{\HH^{\UUU}}(\beta,\theta)\textit{ is a winning strategy for the player }\2\textit{ in }\\
&\textit{the game }\G^\beta_\theta(A^{\HH^{\UUU}},B^{\HH^{\UUU}})\q.
\end{align*}

We now fix some parameters $\beta$ and $\theta$ below $e_{0,\omega}(\kappa)$ in order to consider the game $\G^\beta_\theta(A^{\HH^{\UUU}},B^{\HH^{\UUU}})$, computed in $\HH^{\UUU}$.

Let $\beta:=e_{1,\omega}(\kappa)$ and for each $n$, denote
\[
\beta_n:=e_{n+2,\omega}(\kappa).
\]

The sequence $(\beta_n)_n$ is a descending sequence of ordinals of $\HH^{\UUU}$ below $\beta$.

For each $n$, denote

\[
\bar{X}_n:=e_{n+1,\omega}(X_n).
\]

By construction, the sets $(\bar{X}_n)_n$ cover the domains of the models $\AA^{\UUU}$ and $\BB^{\UUU}$.

\[
\theta:=\max\{|\bar{X}_n|^{\HH^{\UUU}},\beta\}.
\]

Both $\beta$ and $\theta$ are below $e_{0,\omega}(\kappa)$, and each $\bar{X}_n$ has size $\leq\theta$ in $\HH^{\UUU}$.

Then we describe a play of the player $\1$ in the game $\G^\beta_\theta(\AA^{\UUU},\BB^{\UUU})$:

* At the $(2n+1)$th step, he plays the ordinal $\beta_{2n+1}$ and the set $\AA^{\UUU}\cap \bar{X}_{2n+1}$.

* At the $(2n+2)$th step, he plays the ordinal $\beta_{2n+2}$ and the set $\BB^{\UUU}\cap \bar{X}_{2n+2}$.

Every finite initial segment of this play is as an element in the model $\HH^{\UUU}$.
Hence, the player $\2$ must be able to win against this play; otherwise, there would be some finite play of the player $\1$ which the player $\2$ loses and this would contradict the fact that in the model $\HH^{\UUU}$, the player $\2$ has a winning strategy in the game $\G^\beta_\theta(A^{\HH^{\UUU}},B^{\HH^{\UUU}})$.

She can thus win, and eventually, after $\omega$ many steps, she will have produced a chain of partial isomorphisms $(\pi_n)_n$ such that

\[
\bigcup_n\pi_n:A^{\HH^{\UUU}}\cong B^{\HH^{\UUU}}.
\]

This ends the proof of the Claim.

(2)$\To$(1): Assume that $\bar{\UU}=(\UU_n)_n$ are $\kappa$-complete ultrafilters, each $\UU_n$ on some set $I_n$, and
We show that $\AA\equiv_{\LL^1_\kappa}\BB$.

\[
\begin{cases} \AA_0:=\AA\\
\AA_{n+1}:=\AA_n^{\UU_n}
\end{cases}
\]

\[
\begin{cases} \BB_0:=\BB\\
\BB_{n+1}:=\BB_n^{\UU_n}.
\end{cases}
\]

Without loss of generality we may identify each $\AA_n$ with its image under the embedding into the direct limit and get that for each $n$,

\[
\AA_n\elem_{\LL_{\kappa\kappa}}\AA_{n+1}\quad\text{and}\quad\AA^{\UUU}\cong\bigcup_n\AA_n.
\]

and similarly for the models $\BB_n$.
The chains $(\AA_n)_n$ and $(\BB_n)_n$ are thus $\LL_{\kappa\kappa}$-elementary, and by the Union Lemma <ref>,

\[
\AA\equiv_{\LL_\kappa^1}\AA^{\UUU}\cong\BB^{\UUU}\equiv_{\LL^1_\kappa}\BB.
\]

This shows that, indeed, $\AA\equiv_{\LL_\kappa^1}\BB$, as wanted.

[1]
J. Barwise and S. Feferman.
Model-Theoretic Logics.
Cambridge University Press, 1985.

[2]
Akihiro Kanamori.
The higher infinite: large cardinals in set theory from their
Springer Science & Business Media, 2008.

[3]
Saharon Shelah.
Nice infinitary logics.
Journal of the American Mathematical Society, 25(2):395–427,

[4]
Saharon Shelah.
Isomorphic limit ultrapowers for infinitary logic.
Israel Journal of Mathematics, 246(1):21–46, 2021.


Keehang KwonmlabelA Daeseong KangnlabelB [labelA]The author is a professor of
Computer Eng., DongA University<EMAIL_ADDRESS>[labelB]The author is
a professor of Electronics Eng., DongA University. 11 11 200311

# Removing Qualified Names in Modular Languages

(2003; 2003)

###### keywords:

modules, qualified names.

Although the notion of qualified names is popular in module systems, it causes
severe complications. In this paper, we propose an alternative to qualified
names. The key idea is to import the declarations in other modules to the
current module before they are used. In this way, all the declarations can be
accessed locally. However, this approach is not efficient in memory usage. Our
contribution is the module weakening scheme which allows us to import the
minimal parts. As an example of this approach, we propose a module system for
functional languages.

##  1 Introduction

Modularity is the key technique for dealing with large programs. Most modern
languages (including object-oriented ones) employ qualified names of the form
$m.f$ to access a method $f$ in a module $m$. Although the notion of qualified
names is easy to implement, it leads to unnecessarily long names and it runs
counter to the core of knowledgebase representation, i.e., the conciseness.
Thisl problem should be eliminated to preserve clean and concise codes. In
order to do so, we consider here an alternative to qualified names.

The key idea is to import the declarations in $m$ to the current module before
they are used. In this way, all the declarations can be accessed locally.

To be specific, we propose to add the following:

  * •

import declarations (MI declarations) of the form

$/m$

where $/m$ is a module name.

  * •

interactive/querying declarations (MQ expressions) of the form (inspired by
the work in [4])

$(f(t_{1},\ldots,t_{n})=V)^{/m}$

where $f(t_{1},\ldots,t_{n})=V$ is a query to a module $m$ and $V$ is a free
variable not appearing in $x_{1},\ldots,x_{n}$.

The former one has the following semantics: the declarations in $/m$ are to be
added to the current program. This expression thus supports the idea of
importing all the declarations of a module.

The latter one has the following semantics: first evaluate
$f(t_{1},\ldots,t_{n})$ w.r.t. $/m$ and set $V$ to the resulting value $w$.
This expression thus supports the idea of importing some ( logical consequence
of) declarations of a module. This is called module weakening/module
customization. Note that the notion of MQ declarations is a novel feature
which is not present in traditional languages. For example, $(fib(3)=V)^{mf}$
is a querying declaration where the value of $v$ is not known. Its value is
later determined by evaluating $fib(3)$ w.r.t. the “fibonacci” module $mf$.

Although our approach can be applied to other programming paradigms, this
paper focus on functional languages. That is, we extend a functional language
with MI/MQ declarations.

##  2 A Case Study: Functional Languages

The theory of recursive functions, which we call $REC$, provides a basis for
functional programming. It includes operations of composition, recursion,
$etc$. Although $REC$ is quite expressive, it does not contain local module
mechanisms.

To fix this problem, a modern way to add local declarations is to introduce
declarations-implication expressions (DI expressions) of the form:

$D\hskip 2.0pt\mbox{\small$\rightarrow$}\hskip 2.0ptE$

where $D$ is a set of function declarations. This expression is adapted from
the work in [8].

The above has the following intended semantics: the MI/MQ declarations in $D$
are first processed and then are to be added to the current program in the
course of evaluating $E$.

This paper proposes an extension of the core functional languages with
DI/MI/MQ expressions.

##  3 Examples

We assume a fibonacci module $mf$ and a prime module $mp$ which respectively
contains the definitions of $fib(n)$ and $prime(n)$.

% fibonacci module.
---
% $fib(n)$ returns $n$th Fibonacci number.
/mf = % fibonacci module
$fib(1)=1$.
$fib(2)=1$.
$fib(n+2)=fib(n)+fib(n+1)$

% module prime
---
% prime(n) returns true if $n$ is prime
/mp = % prime module
$prime(n)=prime\\_aux(n,n-1)$.
$prime\\_aux(X)=\ldots$.

An illustration of DI expressions is provided by the following definition of
the function $primefib(n)$ which returns true if $n$th Fibonacci number is
prime:

/mw =
---
$primefib(X)$ ${\rm=}$ |
|  $/mf\hskip 2.0pt\mbox{\small$\rightarrow$}\hskip 2.0pt/mp\hskip
2.0pt\mbox{\small$\rightarrow$}\hskip 2.0ptprime(fib(n))$

The body of the definition above contains DI expressions. As an example,
evaluating $primefib(3)$ would result in adding all the declarations in
$mf,mp$, and then evaluating $prime(fib(3))$. The machine returns true, as 2
is prime.

MQ declarations are intended to add only some (logical consequences) of a
module to the current module. An illustration of this aspect is shown here.

$primefib(n)$ ${\rm=}$ |  $(fib(n)=v)^{/mf}$
---|---
|  $\hskip 2.0pt\mbox{\small$\rightarrow$}\hskip 2.0pt(prime(v)=w)^{/mp}\hskip
2.0pt\mbox{\small$\rightarrow$}\hskip 2.0ptprime(fib(n))$

The body of the definition above contains MQ declarations. As a particular
example, evaluating $primefib(3)$ would result in adding only the two
declarations $fib(3)=2$ and $prime(2)=true$ to the program, and then
evaluating $prime(fib(3))$. That is, $v$ is set to 2, $w$ is set to $true$.

##  4 The Language

The language is a version of the core functional languages with DI/MI/MQ
expressions. It is described by $E$\- and $D$-rules given by the abstract
syntax as follows:

$E::=$ |  $c\;|\;x\;|\;h(E,\ldots,E)\;|\;D\hskip 2.0pt\mbox{\small$\rightarrow$}\hskip 2.0ptE\;|\;T$
---|---
$D::=$ | $/m\;|\;f(t_{1},\ldots,t_{n})=E\;|\;(f(t_{1},\ldots,t_{n})=v)^{/m}\;|\;$
| $D\land D$

In the abstract syntax, $E$ and $D$ denote the expressions and the
definitions, respectively. In the rules above, $c$ is a constant, $x,v$ are
variables, $t$ is a term which is either a variable or a constant, and $m$ is
a module name. A set of function definitions $D$ is called a program in this
language.

We will present the semantics of this language in the style of [7]. It
consists of two steps. The first step is to preprocess and instantiate all the
querying declarations in $D$ by invoking queries to their corresponding
module. The second step is described as a set of rules in Definition 1. The
evaluation strategy assumed by these rules is an eager evaluation. Note that
execution alternates between two phases: the evaluation phase defined by eval
and the backchaining phase by bc.

In the evaluation phase, denoted by $eval(D,E,K)$, the machine tries to
evaluate an expression $E$ from the program $D$, a set of definitions, to get
a value $K$. Note that these rules written in logic-programming style, i.e.,
$eval(D,E,K)$ is true if the evaluation result of $E$ in $D$ is $K$. For
instance, if $E$ is a function call $h$, the machine first evaluates all of
its arguments and then looks for a definition of $h$ in the program in the
backchaining mode.

The rules (1) – (4) describe the backchaining mode, denoted by
$bc(D_{1},D,h,K)$. In the backchaining mode, the machine tries to evaluate a
function call $h$ by using the function definition in the program $D_{1}$.

######  Definition 1.

Let $E$ be an expression and let $D$ be a program. Then the notion of
evaluating $\langle D,E\rangle$ to a value $K$ — $eval(D,E,K)$ — is defined as
follows:

  * (1)

$bc(h(c_{1},\ldots,c_{n})=E,D,h(c_{1},\ldots,c_{n}),K)$
if $eval(D,E,K)$. % switch to evaluation mode.

  * (2)

$bc(D_{1}\land D_{2},D,h(c_{1},\ldots,c_{n}),K)$
if $bc(D_{1},D,h(c_{1},\ldots,c_{n}),K)$. % look for $h$ in $D_{1}$

  * (3)

$bc(D_{1}\land D_{2},D,h(c_{1},\ldots,c_{n}),K)$
if $bc(D_{2},D,h(c_{1},\ldots,c_{n}),K)$. % look for $h$ in $D_{2}$

  * (4)

$bc(h(x_{1},\ldots,x_{n})=E,D,h(c_{1},\ldots,c_{n}),K)$
if
$bc(h(c_{1}/x_{1},\ldots,c_{n}/x_{n})=E^{\prime},D,h(c_{1},\ldots,c_{n}),K)$
where $E^{\prime}=[c_{1}/x_{1},\ldots,c_{n}/x_{n}]E$. % argument passing to
$h$ and $E$.

  * (5)

$eval(D,\top,\top)$. % $\top$ is always a success.

  * (6)

$eval(D,c,c)$. % success if $c$ is a constant.

  * (7)

$eval(D,h(c_{1},\ldots,c_{n}),K)$
if $bc(D,D,h(c_{1},\ldots,c_{n}),K)$. % switch to backchaining by making a
copy of $D$ for a function call.

  * (8)

$eval(D,h(E_{1},\ldots,E_{n}),K)$
if $eval(D,E_{i},c_{i})$ and $eval(D,h(c_{1},\ldots,c_{n}),K)$. % evaluate the
arguments first.

  * (9)

$eval(D,D_{1}\hskip 2.0pt\mbox{\small$\rightarrow$}\hskip 2.0ptE,K)$
if $eval(D\hskip 2.0pt\mbox{\small$\wedge$}\hskip 2.0ptD_{1},E,K)$ % DI
expressions.

  * (10)

$eval(D,D_{1}\hskip 2.0pt\mbox{\small$\wedge$}\hskip 2.0pt/m\hskip
2.0pt\mbox{\small$\wedge$}\hskip 2.0ptD_{2}\hskip
2.0pt\mbox{\small$\rightarrow$}\hskip 2.0ptE,K)$
if $eval(D,D_{1}\hskip 2.0pt\mbox{\small$\wedge$}\hskip 2.0ptD_{3}\hskip
2.0pt\mbox{\small$\wedge$}\hskip 2.0ptD_{2}\hskip
2.0pt\mbox{\small$\rightarrow$}\hskip 2.0ptE,K)$, provided that $D_{3}$ is the
declarations contained in module $m$. % MI expressions.

  * (11)

$eval(D,D_{1}\hskip 2.0pt\mbox{\small$\wedge$}\hskip
2.0pt(f(t_{1},\ldots,t_{n})=v)^{m}\hskip 2.0pt\mbox{\small$\wedge$}\hskip
2.0ptD_{2}\ \hskip 2.0pt\mbox{\small$\rightarrow$}\hskip 2.0pt\ E,K)$
if $eval(D,D_{1}\hskip 2.0pt\mbox{\small$\wedge$}\hskip
2.0pt(f(t_{1},\ldots,t_{n})=w)\hskip 2.0pt\mbox{\small$\wedge$}\hskip
2.0ptD_{2}\hskip 2.0pt\mbox{\small$\rightarrow$}\hskip 2.0ptE[w/v],K)$,
provided that $w$ is the value obtained by evaluating $f(t_{1},\ldots,t_{n})$
w.r.t. the module $m$. % MQ expressions.

Note that, in rule (11), a module that is queried could itself query other
modules. For simplicity, other popular constructs such as if-then-else and
pattern matching are not shown above. If $eval(D,E,K)$ has no derivation, it
returns a failure.

##  5 Conclusion

In this paper, we proposed an extension to functional languages with DI/MI/MQ
expressions. These expressions are particularly useful for avoiding qualified
names in functional languages.

The MQ expressions can be implemented by treating them like exceptions. That
is, when an MQ expression is encountered, suspend the current execution,
switch to another execution ( That is, evaluating a query with respect to
another module) and then resume the suspended execution.

Our ultimate interest is to design a module system for Computability Logic
[2]–[6].

## References

  * [1] G. Japaridze, “Sequential operators in computability logic”, Information and Computation, vol.206, No.12, pp.1443-1475, 2008.
  * [2] G. Japaridze. Introduction to computability logic. Annals of Pure and Applied Logic, 2003, 123(1/3): 1-99.
  * [3] G. Japaridze. Computability logic: a formal theory of interaction. In Interactive Computation: The New Paradigm, Goldin D, Smolka S A, Wegner P (eds.), Springer, 2006, pp. 183-223.
  * [4] G. Japaridze. In the beginning was game semantics. In Games: Unifying Logic, Language and Philosophy, Majer O, Pietarinen A -V, Tulenheimo T (eds.), Springer, 2009, pp. 249-350.
  * [5] G. Japaridze. Towards applied theories based on computability logic. Journal of Symbolic Logic, 2010, 75(2): 565-601.
  * [6] G. Japaridze. On the system CL12 of computability logic. http://arxiv.org/abs/1203.0103, June 2013.
  * [7] E. Komendantskaya and V. Komendantsky, “On uniform proof-theoretical operational semantics for logic programming”, In J.-Y. Beziau and A.Costa-Leite, editors, Perspectives on Universal Logic, pages 379–394. Polimetrica Publisher, 2007.
  * [8] D. Miller, “A logical analysis of modules in logic programming”, Journal of Logic Programming, vol.6, pp.79–108, 1989.


exponential coordinates of the first kind.

First notice that $X_{3}=[X_{1},X_{2}]=X_{1}X_{2}-X_{2}X_{1}$ and
$X_{4}=[X_{1},X_{3}]=X_{1}^{2}X_{2}-2X_{1}X_{2}X_{1}+X_{2}X_{1}^{2}$. Hence,
since $f$ is horizontally affine, $X_{3}f=2X_{1}X_{2}f$, by exploiting that
$(X_{1}X_{2}+X_{2}X_{1})f=0$. Moreover, exploiting
$(X_{1}X_{2}+X_{2}X_{1})f=X_{1}^{2}f=0$ we get $X_{4}f=3X_{1}^{2}X_{2}f$.

From $[X_{1},X_{4}]=0$ we deduce
$0=[X_{1},X_{4}]f=(3X_{1}^{3}X_{2}-X_{1}^{2}X_{2}X_{1})f=4X_{1}^{3}X_{2}f$,
where in the second equality we are using that $X_{1}^{2}f=0$ and in the third
one we are using that $(X_{1}X_{2}+X_{2}X_{1})f=0$. Thus $X_{1}^{3}X_{2}f=0$.
From $[X_{2},X_{3}]=0$ we deduce
$0=[X_{2},X_{3}]f=(2X_{2}X_{1}X_{2}+X_{2}X_{1}X_{2})f=3X_{2}X_{1}X_{2}f$,
where in the second equality we are using that $X_{2}^{2}f=0$. Hence
$X_{2}X_{1}X_{2}f=0$. From $[X_{2},X_{4}]=0$ we deduce
$0=[X_{2},X_{4}]f=(3X_{2}X_{1}^{2}X_{2}+2X_{1}X_{2}X_{1}X_{2}-X_{2}X_{1}^{2}X_{2})f=2X_{2}X_{1}^{2}X_{2}f$,
where in the second equality we are using $X_{2}^{2}f=0$ and in the third one
we are using $X_{2}X_{1}X_{2}f=0$, which we obtained before. Then
$X_{2}X_{1}^{2}X_{2}f=0$.

Since $X_{1}(X_{1}^{2}X_{2}f)=X_{2}(X_{1}^{2}X_{2}f)=0$, we get that
$X_{1}^{2}X_{2}f$ is constant. Thus there exists $k\in\mathbb{R}$ such that
$X_{1}^{2}X_{2}f\equiv k$. Then $X_{1}(X_{1}X_{2}f)\equiv k$ and
$X_{2}(X_{1}X_{2}f)=0$ readily imply that that there exists $h\in\mathbb{R}$
such that $X_{1}X_{2}f=kx_{1}+h$ in exponential coordinates of the second kind
described above. Thus $X_{1}(X_{2}f)=kx_{1}+h$ and $X_{2}(X_{2}f)=0$ readily
imply that there exists $v\in\mathbb{R}$ such that
$X_{2}f=kx_{1}^{2}/2+hx_{1}+v$. Moreover,
$X_{2}(X_{1}f)=-X_{1}X_{2}f=-kx_{1}-h$ and $X_{1}(X_{1}f)=0$ yield
$X_{3}(X_{1}f)\equiv-k$ and $X_{4}(X_{1}f)\equiv 0$ so that by integrating the
system of PDEs we obtain that there exists $m\in\mathbb{R}$ such that
$X_{1}f=-kx_{3}-hx_{2}+m$.

Thus if $f$ is horizontally affine, there exist $k,h,v,m\in\mathbb{R}$ such
that $X_{2}f=kx_{1}^{2}/2+hx_{1}+v$ and $X_{1}f=-kx_{3}-hx_{2}+m$. Thus we
obtain that $X_{3}f=(X_{1}X_{2}-X_{2}X_{1})f=2kx_{1}+2h$ and
$X_{4}f=(X_{1}X_{3}-X_{3}X_{1})f=3k$. Thus integrating the system of PDEs one
obtains that there exists $n\in\mathbb{R}$ such that
$f=k(3x_{4}-x_{1}x_{3})+h(2x_{3}-x_{1}x_{2})+vx_{2}+mx_{1}+n$. Thus if $f$ is
horizontally affine,
$f\in\mathrm{span}\\{1,x_{1},x_{2},2x_{3}-x_{1}x_{2},3x_{4}-x_{1}x_{3}\\}$ and
it is readily verified that each element of the previous vector space is
actually horizontally affine, and thus this is a characterization of the
horizontally affine maps in $\mathbb{E}^{1}$.

We can check, through simple computations involving BCH formula, that

$\exp(x_{4}X_{4})\exp(x_{3}X_{3})\exp(x_{2}X_{2})\exp(x_{1}X_{1})=\exp(a_{1}X_{1}+a_{2}X_{2}+a_{3}X_{3}+a_{4}X_{4}),$

implies that $x_{1}=a_{1}$, $x_{2}=a_{2}$, $x_{3}=a_{3}+a_{1}a_{2}/2$, and
$x_{4}=a_{4}+a_{1}a_{3}/2+a_{1}^{2}a_{2}/6$. Thus one obtains, by using
exponential coordinates of the first kind, that $f$ is horizontally affine if
and only if $f\in\mathrm{span}\\{1,a_{1},a_{2},a_{3},6a_{4}+a_{1}a_{3}\\}$. As
a consequence, already in the easiest step-3 Carnot group, one has a
horizontally affine function that is not affine in exponential coordinates of
the first kind, namely
$\widetilde{f}(a_{1},a_{2},a_{3},a_{4}):=6a_{4}+a_{1}a_{3}$. As a consequence
the sublevel sets of $\widetilde{f}$ are precisely monotone sets that are not
half-spaces.

One can also write down the explicit expression of a family of
$\\{X_{1},X_{2}\\}$-polynomial distributions $f$. It can be proved through
some computations involving the explicit expressions of $X_{1},X_{2}$ above
that the vector space of the distributions $f$ on $\mathbb{E}^{1}$ such that
$X_{1}f=X_{2}^{2}f=0$ is ${\rm
span}\\{1,x_{2},x_{3},x_{4},x_{2}x_{4}-x_{3}^{2}/2\\}$, where the functions
are written in exponential coordinates of the second kind associated to
$(X_{1},X_{2},X_{3},X_{4})$.

Free group of step 3 and rank 2. Let $\mathbb{F}_{23}$ be the free Carnot
group of step 3 and rank 2, with Lie algebra $\mathfrak{f}_{23}$ equipped with
the stratification

$\mathfrak{f}_{23}:=\mathrm{span}\\{X_{1},X_{2}\\}\oplus\mathrm{span}\\{X_{3}\\}\oplus\mathrm{span}\\{X_{4},X_{5}\\},$

with nontrivial bracket relations being $[X_{2},X_{1}]=X_{3}$,
$[X_{3},X_{1}]=X_{4}$, $[X_{3},X_{2}]=X_{5}$. In exponential coordinates of
the second kind associated to the adapted basis
$(X_{1},X_{2},X_{3},X_{4},X_{5})$ we can write

$X_{1}=\partial_{1},\quad
X_{2}=\partial_{2}-x_{1}\partial_{3}+(x_{1}^{2}/2)\partial_{4}+x_{1}x_{2}\partial_{5},\quad
X_{3}=\partial_{3}-x_{1}\partial_{4}-x_{2}\partial_{5},\quad
X_{4}=\partial_{4},\,X_{5}=\partial_{5},$

see [BLD19, pages 21-22].

It can be shown through tedious computations involving the explicit
expressions of $X_{1},X_{2}$ above that if a distribution $f$ on
$\mathbb{F}_{23}$ is such that $X_{1}f=0$, and $X_{2}^{2}f=0$, then $f$ is
represented by a polynomial in exponential chart (this comes from Theorem 1.3)
and the vector space of such $f$’s is
$\mathrm{span}\\{1,x_{2},x_{3},x_{4},x_{2}x_{4}-x_{3}^{2}/2,x_{5}+x_{2}x_{3}/2\\}$.

SL(2,$\mathbb{R}$). Let $\mathrm{SL}(2,\mathbb{R})$ be the group of $2\times
2$ real matrices with determinant equal to one. Every element of
$\mathrm{SL}(2,\mathbb{R})$ in a neighbourhood of the identity can be written
as

$\begin{bmatrix}x_{1}&x_{2}\\\ x_{3}&\frac{1+x_{2}x_{3}}{x_{1}}\\\
\end{bmatrix},$

for some $(x_{1},x_{2},x_{3})$ in a neighbourhood of $(1,0,0)$. Thus we can
use $(x_{1},x_{2},x_{3})$ as coordinates from an open neighbourhood of
$(1,0,0)$ in $\mathbb{R}^{3}$ to an open neighbourhood of the identity matrix
in $\mathrm{SL}(2,\mathbb{R})$. It can be computed, see [Tor12, Example 7.16],
that, in such a neighbourhood of $(1,0,0)$, the left-invariant vector fields
$X_{1},X_{2},X_{3}$ such that
$(X_{i})_{|_{(1,0,0)}}=(\partial_{i})_{|_{(1,0,0)}}$ for all $1\leq i\leq 3$,
are

$X_{1}=x_{1}\partial_{1}-x_{2}\partial_{2}+x_{3}\partial_{3},\quad
X_{2}=x_{1}\partial_{2},\quad
X_{3}=x_{2}\partial_{1}+\frac{1+x_{2}x_{3}}{x_{1}}\partial_{3}.$

We have that $[X_{2},X_{3}]=X_{1}$, and then $\\{X_{2},X_{3}\\}$ Lie generates
the lie algebra of $\mathrm{SL}(2,\mathbb{R})$. The coordinate function
$x_{3}$ satisfies $X_{2}x_{3}=0$ and
$X_{3}^{2}x_{3}=X_{3}((1+x_{2}x_{3})/x_{1})=-x_{2}(1+x_{2}x_{3})/x_{1}^{2}+(1+x_{2}x_{3})/x_{1}\cdot
x_{2}/x_{1}=0$. Thus the coordinate function $x_{3}$ is
$\\{X_{2},X_{3}\\}$-polynomial but $X_{1}^{k}x_{3}=x_{3}$ for every $k\geq 0$,
so that in the previous coordinates $x_{3}$ is not a polynomial à la Leibman
on $\mathrm{SL}(2,\mathbb{R})$ even if it is $\\{X_{1},X_{2}\\}$-polynomial.

Orientation-preserving affine functions on $\mathbb{R}$. This example shows
that a $k$-polynomial distribution with respect to a subset $S$ that Lie
generates $\mathfrak{g}$ may not be a polynomial distribution according to
Section 3.1, and thus it may not be a polynomial à la Leibman, see Theorem
1.2. Let us consider the Lie group of orientation-preserving affinity of the
real line

$\mathrm{Aff}^{+}(\mathbb{R}):=\\{\psi:t\in\mathbb{R}\mapsto
yt+x:x,y\in\mathbb{R},y>0\\},$ (A.1)

endowed with the product

$(x,y)\cdot(\overline{x},\overline{y})=(y\overline{x}+x,y\overline{y}),$

that comes from the composition of maps. We identify
$\mathrm{Aff}^{+}(\mathbb{R})$ with $\mathbb{R}\times(0,+\infty)$ by means of
the choice of coordinates $(x,y)$. The identity element of the group is
$(0,1)$ and the left-invariant vector fields $X,Y$ such that

$X_{|_{(0,1)}}=(\partial_{x})_{|_{(0,1)}},\qquad
Y_{|_{(0,1)}}=(\partial_{y})_{|_{(0,1)}},$

are

$X_{|_{(x_{0},y_{0})}}=y_{0}(\partial_{x})_{|_{(x_{0},y_{0})}},\qquad
Y_{|_{(x_{0},y_{0})}}=y_{0}(\partial_{y})_{|_{(x_{0},y_{0})}},\qquad\text{for
all $(x_{0},y_{0})\in\mathbb{R}\times(0,+\infty)$}.$

We claim that the analytic function $f(x,y):=(x+1)\log y$ on
$\mathrm{Aff}^{+}(\mathbb{R})$ is $2$-polynomial with respect to $\\{X,Y\\}$
but it is not polynomial according to Section 3.1, and thus it is not a
polynomial à la Leibman, see Theorem 1.2. Indeed, first
$X^{2}f=(y\partial_{x})^{2}((x+1)\log y)=0$, and
$Y^{2}f=(y\partial_{y})^{2}((x+1)\log y)=0$, and then $f$ is 2-polynomial with
respect to $\\{X,Y\\}$. Second, notice that for every $\alpha\in\mathbb{N}$ we
have, by induction, that
$Y^{\alpha}Xf=(y\partial_{y})^{\alpha}(y\partial_{x})((x+1)\log
y)=(y\partial_{y})^{\alpha}(y\log y)=y\log y+\alpha y$: thus $f$ cannot be a
polynomial according to Section 3.1, and then it is not a polynomial à la
Leibman, see Theorem 1.2.

Let us claim moreover that $f$ is not $\mathfrak{g}$-polynomial, even if it is
$S$-polynomial, thus showing that there is no propagation of the property of
being $S$-polynomial with a Lie generating $S$ in the non-nilpotent case.
Indeed, it can be proved by induction that $(X+Y)^{n}f=y\log y+ny$ for every
$n\geq 2$, and $(X+Y)f=y\log y+x+1$. Thus there does not exist any $n\geq 0$
such that $(X+Y)^{n}f\equiv 0$.

Let us further notice that the map $(\alpha,\beta)\to f\circ\exp(\alpha
X+\beta Y)$ is not polynomial. Indeed, simple computations lead to show that
$\exp(\alpha X+\beta Y)=(\alpha/\beta(e^{\beta}-1),e^{\beta})$ for every
$(\alpha,\beta)\in\mathbb{R}\times(\mathbb{R}\setminus\\{0\\})$, while
$\exp(\alpha X)=(\alpha,1)$, for every $\alpha\in\mathbb{R}$. Then
$f\circ\exp(\alpha X+\beta Y)=\alpha(e^{\beta}-1)+\beta$, for every
$(\alpha,\beta)\in\mathbb{R}^{2}$, which is not a polynomial in
$(\alpha,\beta)\in\mathbb{R}^{2}$.

## References

  * [ABB19] Agrachev, A.; Barilari, D.; Boscain, U.: A Comprehensive Introduction to Sub-Riemannian Geometry. Cambridge Studies in Advanced Mathematics. Cambridge University Press, 2019.
  * [BLD19] Bellettini, C.; Le Donne, E.: Sets with constant normal in Carnot groups: properties and examples. Preprint on arXiv, arXiv:1910.12117 (2019).
  * [BLU07] Bonfiglioli, A.; Lanconelli, E.; Uguzzoni, F.: Stratified Lie Groups and Potential Theory for their Sub-Laplacians. Springer Monographs in Mathematics, Springer-Verlag Berlin Heidelberg, New York, (2007).
  * [Buc70] Buckley, J. T.: Polynomial functions and wreath products. Illinois J. Math. 14 (1970) 274–282.
  * [CK10] Cheeger, J.; Kleiner, B.: Metric differentiation, monotonicity and maps to $L^{1}$. Invent. Math. 182 (2010), no. 2, 335–370.
  * [CKN11] Cheeger, J.; Kleiner, B.; Naor, A.: Compression bounds for Lipschitz maps from the Heisenberg group to $L^{1}$. Acta Math. 207 (2011), no. 2, 291–373.
  * [CG90] Corwin, L. G.; Greenleaf, F. P.: Representations of nilpotent Lie groups and their applications. Part I: Basic theory and examples. Volume 18 of Cambridge Studies in Advanced Mathematics. Cambridge University Press, Cambridge, 1990.
  * [Ehr56] Ehrenpreis, L.: Some properties of distributions on Lie groups. Pacific J. Math. 6 (1956), 591–605.
  * [Fol95] Folland, G. B.: A course in abstract harmonic analysis. Studies in Advanced Mathematics. CRC Press, Boca Raton, FL (1995).
  * [FS82] Folland, G. B.; Stein, E. M.: Hardy spaces on homogeneous groups. Volume 28, Princeton University Press, (1982).
  * [FG02] Fritzsche, K.; Grauert, K.: From holomorphic functions to complex manifolds. Graduate Text in Mathematics, Springer (2002).
  * [GR15] Garofalo, N.; Rotz, K.: Properties of a frequency of Almgren type for harmonic functions in Carnot groups. Calc. Variations and Partial Differential Equations, 54 (2015), no. 2, 2197-2238.
  * [GT06] Green, B. J.; Tao, T.C.: The quantitative behaviour of polynomial orbits on nilmanifolds, Annals of Math. 175 (2012), no. 2, 465–540.
  * [Hor67] Hörmander, L.: Hypoelliptic second order differential equations. Acta Math., 119 (1967), 147-171.
  * [KP20] Kyed, D.; Densing Petersen, H.: Polynomial cohomology and polynomial maps on nilpotent groups. Glasgow Mathematical Journal, 62(3) (2020), 706-736. doi:10.1017/S0017089519000429
  * [LD17] Le Donne, E.: Primer on Carnot Groups: Homogeneous groups, Carnot-Carathéodory Spaces, and Regularity of Their Isometries. Anal. Geom. Metr. Spaces 5 (2017), 116-137.
  * [LDMR20] Le Donne, E.; Morbidelli, D.; Rigot, S.: Horizontally affine maps on step-two Carnot groups. Preprint on arXiv, arXiv:2004.08129 (2020).
  * [LDNG19] Le Donne, E.; Nicolussi Golo, S.: Metric Lie groups admitting dilations. Accepted in Arkiv för Matematik. Preprint on arXiv, arXiv:1901.02559 (2019).
  * [LDT20] Le Donne, E.; Tripaldi, F.: A cornucopia of Carnot groups in low dimensions. Preprint on arXiv, arXiv:2008.12356.
  * [Lei02] Leibman, A.: Polynomial mappings of groups. Israel J. Math. 129 (2002), 29-60. MR 1910931. Zbl 1007.20035. http://dx.doi.org/10.100.
  * [Mon02] Montgomery, R.: A tour of subriemannian geometries, their geodesics and applications. Mathematical Surveys and Monographs, vol. 91, American Mathematical Society, Prov- idence, RI, (2002).
  * [Mor18] Morbidelli, D.: On the inner cone property for convex sets in two-step carnot groups, with applications to monotone sets. Publ. Mat. (to appear), arXiv e-prints, arXiv:1808.06513 (2018).
  * [NY18] Naor, A.; Young, R.: Vertical perimeter versus horizontal perimeter. Ann. of Math. (2) 188 (2018), no. 1, 171–279.
  * [Pan89] Pansu, P.: Métriques de Carnot-Carathéodory et quasiisométries des espaces symétriques de rang un. Ann. of Math., 129 (1989), 1-60.
  * [Pas68] Passi, I. B. S.: Polynomial maps on groups. Journal of Algebra 9 (1968), 121–151.
  * [Ric06] Rickly, M.: First-order regularity of convex functions on Carnot groups. J. Geom. Anal. 16 (2006), no. 4, 679–702.
  * [Rud91] Rudin, W.: Functional Analysis. McGraw-Hill (1991).
  * [Sze85] Székelyhidi, L.: Regularity properties of polynomials on groups. Acta Math. Hungar. 45 (1985) 15–19.
  * [Tor12] Torres del Castillo, G. F.: Differentiable manifolds. A theoretical physics approach. Birkhäuser/Springer, New York, (2012), viii+275.
  * [Var84] Varadarajan, V. S.: Lie Groups, Lie Algebras, and Their Representations. Volume 102 of Graduate Texts in Mathematics (1984).