# Geometric mechanics of random kirigami

Gaurav Chaudhary  School of Engineering and Applied Sciences, Harvard
University, Cambridge, MA 02138.  Lauren Niu  Department of Physics, Harvard
University, Cambridge, MA 02143.  Marta Lewicka  Department of Mathematics,
University of Pittsburgh, Pittsburgh, PA 15260.  Qing Han  Department of
Mathematics, University of Notre Dame, South Bend, IN 46556.  L Mahadevan
<EMAIL_ADDRESS>School of Engineering and Applied Sciences, Harvard
University, Cambridge, MA 02138.  Department of Physics, Harvard University,
Cambridge, MA 02143.

###### Abstract

The presence of cuts in a thin planar sheet can dramatically alter its
mechanical and geometrical response to loading, as the cuts allow the sheet to
deform strongly in the third dimension. We use numerical experiments to
characterize the geometric mechanics of kirigamized sheets as a function of
the number, size and orientation of cuts. We show that the geometry of
mechanically loaded sheets can be approximated as a composition of simple
developable units: flats, cylinders, cones and compressed Elasticae. This
geometric construction yields simple scaling laws for the mechanical response
of the sheet in both the weak and strongly deformed limit. In the ultimately
stretched limit, this further leads to a theorem on the nature and form of
geodesics in an arbitrary kirigami pattern, consistent with observations and
simulations. By varying the shape and size of the geodesic in a kirigamized
sheet, we show that we can control the deployment trajectory of the sheet, and
thence its functional properties as a robotic gripper or a soft light window.
Overall our study of random kirigami sets the stage for controlling the shape
and shielding the stresses in thin sheets using cuts.

###### pacs:

Valid PACS appear here

Kirigami, the art of paper cutting, is now increasingly being seen as a
paradigm for the design of mechanical metamaterials that exhibit exceptional
geometric and structural properties Zhang  _et al._ (2015); Blees  _et al._
(2015); Bertoldi  _et al._ (2017). The basis for kirigami is the well known
observation that the mechanical response of thin sheets is exclusively due to
the geometrical scale separation induced by slenderness which makes bending
deformations inexpensive compared to stretching. In kirigami, the presence of
cuts provides for an extra degree of control via the internal localization of
large bending deformations at cuts; these allow for the nature and scale of
internal large-scale bending modes by varying the number, size and location of
the cuts. This raises a number of questions associated with both the forward
problem of understanding the mechanics of these topologically and
geometrically complex materials as well as the inverse problem of designing
the cuts to obtain different types of articulated deformations for shape
optimization. Recent work in the context of the forward problem has focused
primarily on the mechanics of kirigami with simple distributions of periodic
cuts, aimed at characterizing the response using a combination of theory,
experiment and computation Rafsanjani and Bertoldi (2017); Moshe  _et al._
(2019); Sadik and Dias (2021). In contrast, the inverse problem of designing
cuts that allow for articulated shape transformations has been limited
primarily to geometric optimization Choi  _et al._ (2019, 2021), without much
in the nature of the mechanical response of the resulting structures. To
design cut patterns to control the shape and response of kirigamized sheets,
we need to combine aspects of both these classes of problems by understanding
the geometric mechanics of sheets with multiple, aperiodic cuts. Here, we take
a step in this direction by describing the geometry and mechanics of kirigami
sheets with aperiodic, randomly located cuts in the dilute limit, so that the
cuts do not themselves intersect.

## Geometry and mechanics of a sheet with a single cut

To get a sense of the geometry of a kirigamized sheet, in Fig. 1 a, we show
the shape of a thin circular sheet of radius $R$, thickness $h$ ($R/h\gg 1$)
and a single horizontal cut of length $l$. When the nominal strain $\gamma>0$
induced by the applied vertical force crosses a threshold, the initially
planar sheet buckles out of the plane into a complex geometrical shape. If the
sheet has multiple cuts, as in Fig. 1 b, the deformed geometry is even more
complex. However, casual observations of the sheet show that the underlying
constituents of the deformed sheet are actually simple conical and cylindrical
domains connected by transition layers; increasing the topological complexity
of the sheet increases the number, size, shape and orientational order of such
domains.

Figure 1: (a) When a thin circular paper sheet ($R=8.9$ cm, $h=0.01$ cm) with
a cut of length $l=6.4$ cm) is loaded along diametrically opposite ends, it
buckles into a complex 3D surface, but can be approximated as a composition of
two cones (black dashed lines), and two cylindrical cores (in dotted lines).
(b) A perforated sheet with the same dimensions as in (a) with several cuts
deforms into an even more beautiful 3D structure when pulled as shown. Figure
2: (a) The mean curvature of a loaded sheet with two point loads at $O_{1}$
and $O_{2}$ for various applied strains $\gamma$. The grey lines are the
geodesics connecting the loading points without intersecting the cuts, and the
blue lines are tracking the cut shape. Highly localized bending deformation,
as indicated by the mean curvature, can be observed. (b) Geodesic connecting
the loading points $O_{1}$, $O_{2}$ (extracted from Fig. 2 a) better align
with the axis $O_{1}O_{2}$ with increasing strain. (c) The evolving shape of
the cut shows resemblance to Euler’s elastica. (d) Force displacement curves
for the case shown in (a) with $l/R=0.86$ and $\theta=\pi/2$ subjected to
point load at the opposite ends as shown in the inset. The curves show three
distinct regimes corresponding to the cases shown in Fig. 2 a. Green curve
corresponds to the case shown in Fig. 2 a and red curve corresponds to the
case with cut edges buckling out-of-plane in opposite direction. The rescaled
strain in a form similar to “finite extensibiltiy” models shows a power-law
scaling with an exponent $-1$. (e) A zoomed in view of the cut corner for the
case Fig. 2 a (iii), and a schematic showing the relevant length scales that
govern the mechanics in this regime. (f) Scaling force with
($1/(l^{3}\sin^{3}{\theta})$) collapses the plateau mechanical response for
the cases with varying cut length $l$ and orientation $\theta$. Unscaled data
shown in the inset. Red curves correspond to the cases with $\theta=\pi/2$,
and $l/R=0.86,1.29,1.71$. Green curves are the cases with
$\theta=\pi/4,\pi/3,5\pi/12$ and $l/R=0.86$. (d) Strongly strained kirigami
sheets qualitatively show a “freely-jointed” polymer-chain like divergence as
the length $O_{1}O_{2}\to l_{g}$. Shown are the cases with $\theta=\pi/2$ and
$l/R-[0.57-1.71]$.

We use numerical simulations to study the geometrical mechanics of kirigamized
sheets (see SI S2), starting from the Föppl-von Kármán plate energy of a
triangulated plate that is minimized using a conjugate gradient method. Cuts
are defined as thin rectangular slits of length $l$, and small width $w$, with
a semicircular tip of diameter $w$ added to the tip of the cuts of width $w$.

Figure 2 a shows the mean curvature of the deformed sheet for different values
of the strain $\gamma$ when the cut is initially orthogonal ($\theta=\pi/2$)
to the loading axis (the line connecting the loading points) (Fig. 2 a (i)).
It is evident that for very low strains $\gamma\ll 1$, the sheet just
stretches, but remains planar. When $\gamma>\gamma_{c}$, the sheet buckles
with bending deformations becoming localized along two conical domains
centered near the ends of the cut, and a large cylindrical domain of nearly
uniform mean curvature $\kappa$ appears to connect the cut edges and the
loading points on both sides. Further stretching of the sheet causes an
increase in the curvature of these cylindrical domains, and reduces the
Euclidean distance between the end-points ($A$ and $B$) of the cuts to $S$
from its initial length $l$ (see SI Movie S1). This deformation localization
enables the straightening (in three dimensions) of the lines connecting the
ends of the cuts to the points of force application $O_{1},O_{2}$ shown in
Fig. 2 a, as shown in Fig. 2 b. Simultaneously, the free edges associated with
the cut deform into a shape that resembles Euler’s elastica Love (2013) as
shown in the Fig. 2 c. When the applied strain becomes very large, the end
points of these elasticae, corresponding to the ends of the cut, come together
so that the polygons connecting the ends of the cut to the points of force
application straighten out and merge. If the sheet thickness is very small
($h/R\rightarrow 0$), we thus expect that the sheet will be locally flat-
folded onto itself, a limit that we will address later.

Moving from this geometric description of the sheet to its mechanical
response, the force-displacement response of the kirigami sheet described in
Fig. 2 a-c is shown in Fig. 2 d. At very low values of strain $\gamma$, the
force is linearly proportional to the strain, but once the sheet buckles, the
force flattens out to a plateau-like regime as the sheet stretches by bending
out the plane. Eventually, as the ends of the cut come together, the sheet
stiffens as it cannot deform any further without significant stretching and
the force increases showing a near-divergence. When the sheet buckles out of
the plane, there are two almost equivalent modes of deformation: a symmetric
mode when both cylindrical domains on either side of the cut are in-phase, and
an antisymmetric mode both cylindrical domains on either side of the cut are
out-of-phase. The plateau force for these two cases differ marginally, but the
linear and the divergent response away from the plateau is indistinguishable
as can be seen in Fig. 2 d.

To understand the origin of this divergence, we use a scaling approach. The
torque due to the applied force $f$ acting over a length $S$ corresponding to
the (small) distance between the ends of the cut is balanced by the internal
elastic torque $EI\kappa_{c}$ where $\kappa_{c}$ is the characteristic mean
curvature in the neighborhood of the end of the cut, with $EI$ being the
bending stiffness of the sheet ($E$ is the Youngs’s modulus and $I$ is the
second moment of area around the axis normal to the deformed area near the cut
corner), so that $fS\sim EI\kappa_{c}$. As shown in Fig. 2 e, at the ends of
the cut of small width $w$, the characteristic curvature
$\kappa_{c}\sim\theta_{c}/w$, where $\theta_{c}$ is the angle at the slit
corner. As the sheet is pulled apart by the forces so that they are $2R\gamma$
apart, geometry implies that $\theta_{c}S/2\sim 2R\gamma$ and $S\approx
2\sqrt{(l/2)^{2}-2R^{2}\gamma}$. Substituting these geometric relations into
the overall torque balance then yields the relation

$f\sim\frac{EI}{wR}\frac{\gamma}{(\bar{l})^{2}-\gamma}$ (1)

where $\bar{l}=l/(2\sqrt{2}R)$. Writing the stretch ratio of deformation as
the end-to-end displacement $R(1+\gamma)$ normalized by the length of the
shortest segment connecting the force application points to the ends of the
cut, i.e. the piecewise linear geodesic length $l_{g}=\sqrt{(l/2)^{2}+R^{2}}$,
eq. 1 can be expressed in a more familiar form $\hat{f}\sim
1/(1-R(1+\gamma)/l_{g})$ where $\hat{f}=f/(EIwR)$ which we see is similar to
the divergent response of a freely-jointed polymer chain Grosberg and Khokhlov
(1994). The inset to Fig. 2 d shows that the mechanical response with the
rescaled definition of the stretch agrees remarkably well with this simple
scaling estimate. We pause to note that the divergent mechanical response of
freely-jointed chain is intimately linked to the balance between entropic
effects and a finite chain length, quite unlike the divergent response of the
athermal kirigami sheet which is due to the localization of curvature of the
sheet at the ends of the cut.

Having understood the geometry and mechanical response of a sheet with a
single symmetrically placed cut, we ask what would happen when the length $l$
and/or its orientation $\theta$ is varied. Fig. 2 f inset shows the mechanical
response for various cases where $l$ is varied keeping the cut orientation
orthogonal to the clamped axis ($\theta=\pi/2$) in red curves. The mechanical
response in all cases is qualitatively similar to Fig. 2 d, with a shift in
the applied strain $\gamma$ at the onset of plateau response, the magnitude of
the plateau force, and the strain at the onset of divergent response. Larger
$l$ result in a lower value of the threshold in $\gamma$ and a lower force
plateau persisting for longer, before the force diverges. Similarly changing
the orientation $\theta$ of the cut changes the mechanical response and the
plateau value increases as the initial cut direction is more aligned with the
direction of the clamping axis. The deformed geometric configurations for all
the cases are shown in Fig. S5-6.

To quantify these observations, we note that the linear response at very small
strain corresponds to the planar stretching of sheet. The presence of the cut
leads to a stress intensification near the ends of the cut. Following Inglis’
seminal work Inglis (1913), the stress near the tip of a cut of length $2l$
and a very small radius of curvature $w$ can be approximated as
$\sigma_{0}\sqrt{l/w}$, where $\sigma_{0}$ is the far field applied stress on
a large plate. Assuming the planar elastic deformation is confined near the
crack tip, we get $\sigma_{0}\sqrt{l/w}\sim E\gamma$. Hence we expect
$\sigma_{0}\sim f\sim 1/\sqrt{l}$ (see SI Fig. S4b). Following the initial
linear increase in the force with applied strain, the sheet buckles when the
compressive load on a plate whose size scales as the size of the cut reaches
the buckling threshold. Since the buckling stress $\sigma_{b}\sim B/hl^{2}$,
with $B$ being the flexural rigidity of the sheet, thus the buckling force
$f_{b}\sim(B/h^{2}R)(1/l^{2})$ (see SI Fig. S4d). The strain at which sheet
buckles ($\gamma_{b}$) can be estimated from a balance of the linearly
increasing planar stress ($\sim E\gamma_{b}/\sqrt{l/w}$) with the buckling
stress, so that $\gamma_{b}\sim 1/l^{3/2}$ (see also SI Fig. S5-6).

Following the onset of buckling, for small out-of-plane displacements
$\delta\ll l$ of the cut boundary, the curvature of the deformed cylindrical
core of the sheet scales as $\kappa\sim\delta/l^{2}$. Thus the total bending
energy of the sheet scales as $U\sim B\kappa^{2}A\sim B\delta^{2}R/l^{3}$,
with $A\sim Rl$ being the area of the localization size approximated by area
under quadrilateral $O_{1}AO_{2}B$ shown in Fig. 2 b. Thus the force scales as
$f=\partial U/\partial(R\gamma)\sim BR\delta/l^{3}$. For a cut oriented at an
arbitrary angle $\theta$ to the clamped axis, we replace $l$ by its orthogonal
projection to the loading axis $l\sin{\theta}$. Fig. 2 f shows the collapsed
data from Fig. 2 f inset, indicating that the bending energy localization in
the elasticae corresponding to the two edges of cut determines the magnitude
of plateau response in the force-displacement curves. We note that we have
ignored the contribution of energy from the conical domains since the size of
conical domains is much smaller than the cylinder-like domains just above the
onset of out-of-plane deformation, and further mean curvature decays away from
the cone tip as $\kappa\sim 1/r$. At large applied strain, when the geodesics
are relatively better aligned with the loading axis, application of further
strain induces a strong bending deformation at the cut corners. Sheets with a
finite tearing threshold stress generally tear as a result of this
deformation, and this raises a class of different questions about the nature
and shape of the curve of tearing Zehnder and Potdar (1998). Fig. 2 g depicts
the rescaled force-displacement curves showing the divergence response of
sheets with varying cut length which agree well with the scaling arguments
presented above in eq. 1.

## Geometry and mechanics of a sheet with multiple cuts

Figure 3: (a) Localized bending deformation in the sheets with multiple cuts
aligned orthogonally to the loading axis show similar geometric structures
regardless the number of cuts, shown are front and side views of the cases
with 2 and 12 cuts. (b) Mechanical response for kirigami sheets with different
number parallel of cuts ($\theta=\pi/2,l=0.57R$). (inset) Rescaled force with
the length $d$ (labelled in (a)) collapsed the plateau regime. (c) Localized
deformation in a random kirigami sheet with 21 randomly distributed cuts of
same length (dashed lines) with a total length of all cuts $8R$. The mean
curvature of the deformed sheet is projected on the initial flat
configuration. (d) Mechanical response of a randomly kirigamized sheet
perforated with 21 randomly distributed uniform cuts, with varying total cut
length has similar behavior as sheets with single and multiple structured
cuts. (inset) Rescaled force data collapsed the plateau force reasonably well.

We now turn to understand how the geometric mechanics of sheets kirigamized
with a single cut translates into our understanding of sheets with multiple
cuts, as shown in Fig. 3. As examples, the mean curvature maps in Fig. 3 a
show the similarity in the localized deformation for the cases with 2 and 12
cuts that are perpendicular to the loading axis (see SI Movie S2, other cases
are shown in SI Fig. S9). It is clear that the deformed geometry in such cases
consists of four conical domains, and $N+1$ connected elasticae where $N$ is
the number of cuts. Indeed, irrespective of the number of cuts, the conical
domains localize near the ends of the cuts that are nearest to the loading
points and are connected by a flat sheet for the cases with larger $N$. Fig. 3
b shows the mechanical response of sheets kirigamized with varying number of
cuts and shows that increasing the number of cuts softens the system, reducing
the plateau force and delaying the transition to the ultimately divergent
force-displacement response. To characterize the mechanical response in the
case with multiple cuts, we approximate the bending energy localized in $N+1$
elasticae close to the onset of out-of-plane deformation. This results in the
form the energy $U\sim BR^{2}d\gamma/l^{3}$, and hence $f\sim\partial
U/\partial(R\gamma)\sim BRd/l^{3}$ where $d$ is the distance between the
loading point and the nearest cut (see SI-S3). Consistent with this, rescaling
the force with $d$ collapses the data over the scale of intermediate
deformations as shown in the inset of Fig. 3 b.

At large strains the strong bending deformations near the cone tips similar to
the case with single cut discussed earlier. Predicting the location of conical
domains will enable predicting the stress concentrations, and the potential
sites of structural failure in practical applications. Our observations
suggest that the conical domains appear at the end of a cut if the addition of
that cut increases the geodesic length. To formalize this, we define a binary
participation ratio (_PR_) for each cut as

$PR_{i}=\begin{cases}1&\text{if $l_{g}^{0}-l_{g}^{i}>0$,}\\\ 0&\text{if
$l_{g}^{0}-l_{g}^{i}=0$}\end{cases}$

where $l_{g}^{0}$ is the geodesic length evaluated for a given cut
arrangement, and $l_{g}^{i}$ is the geodesic length with the $i^{\text{th}}$
cut removed from the arrangement. For the cases shown in Fig. 3 a-b, two cuts
that are nearest to the loading points have a $PR=1$ and all other cuts have
$PR=0$. For such cases, any cut with a projected length $l\sin\theta>d$ the
distance of the force application point does not influence $l_{g}$, and thence
the mechanical response. Similarly, for the cases with two cuts of varying
projected lengths (see SI Fig. S7, S8), conical domains disappear at the
corners of cuts with $PR=0$. In Fig. 3 a-b, we see that an increase in the
number of cuts increases the length $l_{g}$. Since the divergent force-
displacement response emerges as $2R(1+\gamma)\to l_{g}$, the plateau response
in the force-displacement curves is observed at larger $\gamma$ for larger
$N$. In cases with $N>2$, fixing the location of the cuts closest to the
points of force application sets the trajectory of geodesics as well as
$l_{g}$, regardless the presence of the inner cuts; hence the divergence
transition occurs at same strain regardless of the presence of the inner cuts.
In fact the force-displacement curves overlap at all applied strains for such
cases, and the geometry of the deformed sheets is identical (see SI Fig. S10).

The geometry of the geodesics controls the mechanical response of the overall
system when the cuts are randomly distributed. In Fig. 3 c, we show an example
with 21 cuts of the same length $8R$, with the location of the cut midpoint
and its orientation randomly chosen such that a minimum separation exists
between the cuts, and from the clamped points (varying the cut length does not
change any of our results qualitatively). The localization of deformation in
multiple elasticae and conical domains is evident while the sheet remains flat
and undeformed near some cuts.

To obtain the average response for the random kirigami cases, we repeat the
simulations keeping $\sum l$, sum of all cut lengths, a constant. Fig. 3 d
shows the mechanical response for random kirigami. The results represent mean
statistics of the response for 10 samples per case. The overall nature of the
force-displacement curve is similar to that of a single cut, and as expected
the plateau force decreases for the cases with longer cuts. The cases with
smaller $\sum l$ show a very weak deviation from the initial linear response,
but with increasing $\sum l$ a clear plateau is observed. The force-
displacement divergence is observed when the geodesic connecting the points of
force application straightens out under applied strain. Relatively large
variance in force beyond the initial linear response exists. This is due to
the strong dependence of mechanics on the cut length and location beyond the
initial linear regime. Similarly, increasing the number of cuts while keeping
$\sum l$ constant, results in similar observations (see SI Fig. S12) with a
lower number of longer length cuts resulting in a lower plateau force response
and an increased variance.

It is evident from the geometry of random kirigami that the deformation gets
localized near a few cuts (see SI Fig. S11). Just as for the case with
structured cuts, where the mechanical response depends on the distance of cut
from the point of application of force, its projected length and its
participation ratio $PR$, for random cuts that are not very close to one
another, a similar scenario arises. At the onset of the buckling transition
from the initial planar stretching response at very low strains, each cut with
$PR>0$ introduces a soft bending deformation mode in the sheet with a
characteristic bending force $f^{*}$, given by the smallest buckling load,
i.e.

$f^{*}\sim Eh^{3}\leavevmode\nobreak\
\min_{i}\left[\frac{d_{i}}{l_{i}^{3}\sin^{3}{\theta_{d_{i}}}}+\frac{2R-d_{i}}{l_{i}^{3}\sin^{3}{\theta_{2R-d_{i}}}}\right]2R$
(2)

where $d_{i}$ is the minimum distance of the midpoint of the cut of length
$l_{i}$ from the loading points, and $2R-d_{i}$ is the distance of the cut
midpoint from the farthest loading point, and $\theta_{d_{i}}$ and
$\theta_{2R-d_{i}}$ are the angles that the cut makes with the line joining
the midpoint of the cut to the points of force application. The above result
follows from the assumption that at the onset of plateau regime, the
characteristic mean curvature for bending localization is set by the cut that
corresponds to $f^{*}$, so that the two terms in eq. 2 follow from the energy
of two elasticae that exist on both sides of the cut. Here, we note that the
cuts with $PR=0$ do not alter the plateau response near its onset as seen with
the cases in Fig. 3 b (see SI Fig. S7-8). These observations allow us to
determine the rescaled mechanical response shown in the inset of Fig. 3 d,
providing a reasonable collapse in the plateau region of force-displacement
data. The spread in the scaled data is likely due to the simplification that
the area of the sheet where energy is localized is assumed to span the sheet
(hence the factor $2R$ in eq. 2), and that additional cuts which buckle
following the onset of first buckling also contribute to localize the bending
deformation. All together, this allows us to reduce a given random kirigamized
sheet to a simpler, “mechanical equivalent” (see SI Fig. S13). Since the cuts
with $PR=0$ do not effect the geometric mechanics of the sheet, simply
removing them from the given cut arrangement results in a sheets with
identically geometric mechanics with reduced number of cuts. Thus our simple
scaling approach reduces the complexity of random kirigami using elementary
geometric mechanics.

In the ultimate deformation limit of sheets with multiple cuts, the shape of
the sheet is determined by the (3-dimensional) straightening of the shortest
path (corresponding to the geodesic) connecting the points of force
application. Geodesics for all the cases have three components: two line
segments connecting the points of application of force to the ends of the cuts
nearest to these points, and a piecewise polygonal geodesic passing through
the ends of all cuts.

## Zero-thickness, flat-folded kirigami

For an initially flat sheet, our observations suggest that the shortest path
between the points of force application for simple cut patterns is just a
polygonal curve that connects these points, i.e. all geodesics in a planar
sheet with random cuts are polygonals. When a very thin sheet is deformed by
boundary forces, its ultimate shape is characterized by the formation of sharp
creases as the sheet folds on itself, as shown in Fig. 5. These observations
of the geometry of strongly deformed kirigamized sheets show that the
polygonal geodesic connecting the points of force application in the plane,
becomes approximately straight $\mathbb{R}^{3}$. When the sheet is flat-
folded, the geodesic is rectified leading to a configuration that is a
piecewise affine isometric immersion of the plane. We leave precise theorems
and proofs of these statements for a separate study Han  _et al._ (2021), but
provide an intuitive argument for them here.

Figure 4: The path-shortening algorithm yields a polygonal competing to be a
geodesic. Figure 5: Strongly stretched inextensible sheet of negligible
bending rigidity, $h\to 0$, can be folded to 2d sheets. This is illustrated
through three different kirigamized sheets. The two projections are
experimentally realized by creasing the sheets: (i) 1 crease along each
cylindrical core, and (ii) 3 creases to flatten the conical domains near the
cut corners. The flat foldability of a strongly stretched sheet leads to
straightening of geodesics (shown in red curve) in 3d space.

We represent the given set of cuts $L$ contained in an open, bounded, convex
domain $\Omega\subset\mathbb{R}^{2}$, as the set of edges of a graph $G$.
Without loss of generality $G$ is planar, i.e. each pair of its edges
intersects at most at a single common vertex.

The first result, i.e. the polygonal structure of geodesics, is shown via a
path-shortening algorithm. Given $p,q\in\bar{\Omega}\setminus L$ and a
piecewise $C^{1}$ curve $\tau:[0,1]\to\bar{\Omega}\setminus L$ with
$\tau(0)=p$ and $\tau(1)=q$, one successively replaces its portions by
segments, as follows. Let $t_{1}\in(0,1]$ be the first time that the segment
$\overline{p\tau(t_{1})}$ intersects $L$. If $\tau(t_{1})=q$, then
$\overline{pq}$ is the desired geodesic connecting $p$ and $q$. Otherwise,
$\overline{p\tau(t_{1})}$ must contain some of the vertices of the cuts. Call
$p_{1}$ the closest one of these vertices to $\tau(t_{1})$ and concatenate the
segment $\overline{p_{1}\tau(t_{1})}$ with the curve $\tau$ restricted to
$[t_{1},1]$. The process is now repeated from $p_{1}$. After finitely many
such steps one obtains a polygonal connecting $p$ and $q$, with length shorter
length than that of $\tau$.

The second result, i.e. existence of a geodesic-rectifying piecewise affine
isometry, follows via the folding algorithm below. We take
$p,q\in\partial\Omega$ and denote by $d$ the length of the geodesics between
$p$ and $q$.

Step 1. Sealing portions of inessential cuts that do not affect
$\mbox{dist}(p,q)$. To this end, label cuts (the edges of $G$) by
$l_{1},\ldots,l_{n}$. Move the first endpoint vertex of $l_{1}$ toward its
second vertex, and start “sealing” the portion of the cut $l_{1}$ left behind.
The length $d$ of the geodesics connecting $p$ and $q$ is nonincreasing: it
may drop initially, decrease continuously, or it may initially remain
constant. The sealing process is stopped when $d$ becomes strictly less than
the original geodesic distance, and the new position point is labeled as the
new vertex endpoint of $l_{1}$. In the next step, the second endpoint is moved
along $l_{1}$ toward the (new) first endpoint and the process is repeated,
thus possibly sealing the cut $l_{1}$ further. The same procedure is carried
out for each $l_{i}$ in the given order $i=1,\ldots,n$. It follows that upon
repeating the same process for the newly created configuration, labeled the
minimal configuration, it will not be further altered.

While different ordering of cuts and vertices may yield different minimal
configurations and new geodesics may be created in the cut-sealing process,
all original geodesics are preserved. Also, since the updated set $L$ is a
subset of the original $L$, finding an isometry relative to the new $L$ yields
an isometry for the original cut set $L$.

Step 2: Ordering the geodesics and the connected components of
$\Omega\setminus L$. There are two important properties of any minimal
configuration: the graph $G$ has no loops (i.e. it is a collection of its
connected components that are trees), and each vertex that is a leaf is a
vertex of some geodesic.

With these properties, one proceeds to label all geodesics in a consecutive
order, with $\tau_{1}\preceq\ldots\preceq\tau_{N}$. Here,
$\tau_{r}\preceq\tau_{r+1}$ means that the concatenated polygonal from $p$ to
$q$ via $\tau_{r}$ and then back to $p$ via $\tau_{r+1}$ encloses a region
$D_{r}$ and it is oriented counterclockwise with respect to $D_{r}$. Next, one
labels and orders the trees $\\{T_{m}\\}_{m=1}^{s}$ in $\bar{D}_{r}$ so that
$D_{r}\setminus L$ is partitioned into subregions $\\{P_{m}\\}_{m=0}^{s}$ and
$\\{Q_{m}\\}_{m=1}^{s}$ in the following way: each $P_{m}$ is a polygon
bounded by the “right most” path from the tree $T_{m}$, the “left most” path
from $T_{m+1}$, and the intermediate portions of $\tau_{r}$ and $\tau_{r+1}$
which are concave with respect to $P_{m}$. Each $Q_{m}$ is a finite union of
polygons enclosed within the single tree $T_{m}$, again bounded by portions of
geodesics $\tau_{r}$ and $\tau_{r+1}$. Note that $\tau_{r}$ and $\tau_{r+1}$
may have nontrivial overlaps and some of $\\{Q_{m}\\}_{m=1}^{s}$ may be empty.

Step 3: Constructing a desired isometry. We fix the segment $I=[0,de_{1}]$
along the $x_{1}$-axis in $\mathbb{R}^{3}$ and construct an isometric
immersion $u$ of $\Omega\setminus L$ into $R^{3}$, such that $u(p)=0$,
$u(q)=de_{1}$, $u(\tau_{r})=I$ for $r=1,\ldots,N$, and where each segment on
$\tau_{r}$ is mapped onto a designated subsegment of $I$.

The resulting $u$ consists exclusively of planar folds and returns the image
that is a subset of $\mathbb{R}^{2}$.By Step 2, for each $r=1,\cdots,N-1$ we
have $D_{r}=\bigcup_{m=0}^{s}P_{m}\cup\bigcup_{m=1}^{s}Q_{m}$. We construct
$u$ separately on $P_{0},Q_{1},Q_{2},\ldots,Q_{m},P_{m}$, where the step to
construct $u$ on $P_{1},\ldots,P_{m-1}$ is highly technical Han  _et al._
(2021). Since the exterior region
$D_{0}=\Omega\setminus\bigcup_{r=0}^{N-1}D_{r}$ does not contain trees, the
two outermost geodesics $\sigma_{1}$ and $\sigma_{N}$ are convex, and so the
definition of $u$ on $D_{0}$ consists of several simple folds.

It also turns out that the condition $p,q\in\partial\Omega$ is essential:
there exist minimal configurations for $p,q\notin\partial\Omega$, that do not
admit any isometry $u$ with the property that the Euclidean distance from
$u(p)$ to $u(q)$ equals the geodesic distance from $p$ to $q$ in
$\Omega\setminus L$ (for further details, we refer to Han  _et al._ (2021)).

## Functional kirigami structures

Figure 6: (a) The kirigami grasper can be employed in active and passive modes
requiring different energy input. A simple framework of grasper function on a
generic mechanical displacement of the sheet. (b) A kirigami grasper design
for grasping objects. The deformed sheet shows the concave regions (for active
grasping), and open holes (for passive grasping). (c) Passsive grasping of an
racquet ball. (d) Active grasping of a tennis ball.

The geometric mechanics of ordered and disordered kirigami leads us naturally
to questions of design for function, which we demonstrate with two examples.
The first is the use of kirigamized sheets to robotic grasping. To be
effective, a grasper must enable controlled gripping, lifting and relocating
objects of varying scales and shapes with minimal external energy input. We
define an active grasper as one that requires continuous application of
tensile force (external energy) to grasp and relocate, while a passive grasper
only requires external work to be done in order to grasp and release the
object, with _energy-free_ relocation as shown in the schematic in Fig. 6 a.
Kirigami enables both designs; while the active gripper accommodates the
target object in the curved features of the deformed sheet, the passive
gripper utilizes the holes/cuts in the structure to accommodate objects. The
passive mode thus requires a prestretch to deform the 2d cut into a 3d
slot/hole of the size similar to the target object. On releasing the
prestretch the cut boundary forms contact with the object, grasping it. The
arrested object is released by applying an extensional strain to the sheet.
The sequence is demonstrated in Fig. 6 c and SI Movie S3. In the passive mode,
the conical tips and the cylindrical core of the kirigami sheet enables
confinement of the target object as shown in Fig. 6 c and SI Movie S3.

The kirigamized grasper has multiple smaller cuts, in addition to the larger
cut in the middle as shown in Fig. 6 b. The additional cuts are prescribed in
a way that all cuts have $PR>0$, and buckle under the applied strain. Further,
the largest cut is split in the middle by a small cut, oriented along the
pulling axis. Together these features enable symmetric deformation of the flat
kirigami sheet, and improve the stability in handling due to additional points
of contact with the grasped object as shown in Fig. 6 b. We note that the
simplest kirigamized gripper with a single cut (Fig. 1 a, 2 a) deforms
asymmetrically, and hence cannot be effectively used as a grasper in both
active and passive modes. In practical scenarios however symmetric deformation
can be realized if the sheet is stretched significantly Yang  _et al._ (2021).

The range of applicability of a kirigami gripper can be understood from a
balance of the forces due to the bending of kirigami sheet, and the weight of
the object to grasp. For a cut of length $l$, the characteristic grasping
force of the deformed sheet can be written as $f_{g}\sim\mu(B/l)c(\gamma)$,
with $\mu$ being the coefficient of friction and $c(\gamma)$ is a strain
dependent geometric factor. For an object with effective density $\rho$ and
size $l$, force balance yields $\mu(B/l)c(\gamma)\sim\rho gl^{3}$. This
provides a non-dimensional kirigami grasper parameter $\rho gl^{4}/(\mu
Bc(\gamma))$, which has values ranging between $1-100$ for the successful
grasping demonstrations using a plastic gripper.

Although our discussion so far was restricted to cuts that are straight
rectangular slits, the ideas established are applicable to general shaped
cuts. Inspired by the recent efforts towards design of kirigami-inspired
mechanically deployable structures Zhang  _et al._ (2015), we demonstrate a
simple deployable kirigami structure here, whose force-dependent shielding and
transmission can be tuned with simple geometric parameters.

Figure 7: (a) Kirigami sheets with curved cuts (dashed lines) show very
different response to an applied strain. Red line is the geodesic connecting
points of application of force $O_{1},O_{2}$. (b) A phase diagram quantifying
the mechanical shielding and transmittance response of the deployed kirigami
structures.

Fig. 7 a shows a case of physical kirigamized sheets perforated with five
concentric circular arcs. The arcs extend between $\pm d\phi$ and
$\pm(\pi-d\phi)$ w.r.t. the vertical axis, and the radius is linearly
increased between the five arcs. The planar kirigami sheets (with marginally
different $d\phi$) show very different geometric mechanics under a small
applied strain. For the case when $d\phi=-\pi/30$, the geodesic connecting
$O_{1}$ and $O_{2}$ skirts the outer cut without intersecting any of the inner
cuts. The inner cuts have $PR=0$, and hence do not introduce any soft
deformation modes. Under an applied deformation, the outer frame localizes the
bending while the inner structure stays nearly planar without any deformation
or stays “mechanically shielded”. A small change in $d\phi(-\pi/180)$ results
in the cut arrangement such that the geodesic connecting $O_{1}$ and $O_{2}$
meanders through the corners of all cuts. Under an applied load, this
structure shows a large relative out-of-plane displacement of different
domains of the sheet (see SI Movie S4). Further changing $d\phi=\pi/60$,
results in a straight geodesic connecting $O_{1}$ and $O_{2}$, and hence
tensile loading of this structure results in the in-plane deformation.

We quantify the functional response of such kirigamized structures using
$d\phi$ as a tunable parameter. A straightforward observation is that force
required to deform such class of structures increases monotonically with
$d\phi$ (see SI Fig. S14a). The geometric consequence of different geodesic
paths can be quantified in terms of two functional features of this family of
kirigami structures: transmittance and shielding (Quantification of functional
responsec:materials). In practical scenarios this corresponds to the light
transmitted though an optical window when illuminated with a light rays
perpendicular to the rest plane. Shielding is proposed to be linked its
ability to restrict the mechanical deformation to the boundary, and
effectively protecting the interior. We restrict the deformation to small
strains $0.01-1\%$, a practically relevant regime. And since the geometric
mechanics is strain-dependent, we take a representative value for comparison
corresponding to the highest applied strain.

Fig. 7 b shows a phase diagram showing shielding and transmittance as a
function of $d\phi$. The properties are evaluated at a strain of 1%. Both
flatness and transmittance curves display a non-monotonic trends. The cases
with large $l_{g}$ show high shielding ability since their interior remains
relatively flat, while a high transmittance is achievable in the cases with
$l_{g}\to 2R$. Further the strain sensitivity of transmittance is enhanced
with increasing $l_{g}/2R$ (see SI Fig. S14b). For the cases with $l_{g}=2R$ a
weak dependence on the functional properties can be seen since the structures
deform primarily by planar stretching.

## Discussion

Our study of random kirigami has shown how elementary geometric and energetic
concepts allow us to understand the three-dimensional structure and mechanical
response of kirigamized sheets. This leads to a geometric view of how cuts
respond or not, along with scaling arguments for all regimes of deformation.
The resulting simplicity of the framework reduces a complex nonlinear problem
to geometrical constructions, and thus eases the search for novel engineering
solutions using kirigami in such instances as grasping and windowing, which we
hope are just the beginning of a different way of thinking about using
topological and geometrical mechanical metamaterials.

## Materials and Methods

### Numerical Simulations

For our numerical experiments, use the finite difference scheme outlined in
Weischedel  _et al._ (2012) and van Rees  _et al._ (2017) to represent the
sheets and minimize the Föppl-von Kármán plate energy on triangular meshes of
the geometry, fixing the two clamped regions at various distances apart, with
otherwise free boundary conditions. The Young (elastic) modulus, Poisson
ratio, and thickness used to calculate the elastic energy are $1$ GPa, $0.4$,
and $0.01$ cm respectively. The geodesics of the mesh were computed using an
open source code Peyré (2011) based on fast marching approach Kimmel and
Sethian (1998).

### Quantification of functional response

Transmission of the kirgami window is quantified as the difference between the
area of undeformed sheet and the projected area of deformed sheet on the rest
plane, normalized by the undeformed sheet area. The shielding effect is
quantified by the average flatness of the sheet. We define flatness as the
variance in the radial correlation distribution of the average face normal of
an area element with the orientation of the element at sheet center,
$(\sum_{i}\hat{n}_{0}.\hat{n}(r_{i})a(r_{i}))/A(r_{i}+\Delta r)$. Here,
$\hat{n}_{0}$ is the unit normal to a infinitesimal element at $r_{i}\to 0$
and $\hat{n}(r_{i})$ is the unit normal to an element at a distance $r_{i}$
from the center with an area $a(r_{i})$. This quantity represents the average
orientation of all elements located between $r_{i}$ and $r_{i}+\Delta r$ with
the orientation of the center.

###### Acknowledgements.

We thank the Bertoldi lab, Jeremy A. Guillette and FAS Academic Technology at
Harvard University for sharing resources that helped our experiments. The work
was supported partially by NSF grants DMS-2006439 (ML), BioMatter DMR 1922321
and MRSEC DMR 2011754 and EFRI 1830901, the Simons Foundation and the Seydoux
Fund (LM).

## References

  * Zhang  _et al._ (2015) Y. Zhang, Z. Yan, K. Nan, D. Xiao, Y. Liu, H. Luan, H. Fu, X. Wang, Q. Yang, J. Wang, _et al._ , Proceedings of the National Academy of Sciences 112, 11757 (2015).
  * Blees  _et al._ (2015) M. K. Blees, A. W. Barnard, P. A. Rose, S. P. Roberts, K. L. McGill, P. Y. Huang, A. R. Ruyack, J. W. Kevek, B. Kobrin, D. A. Muller, _et al._ , Nature 524, 204 (2015).
  * Bertoldi  _et al._ (2017) K. Bertoldi, V. Vitelli, J. Christensen, and M. Van Hecke, Nature Reviews Materials 2, 1 (2017).
  * Rafsanjani and Bertoldi (2017) A. Rafsanjani and K. Bertoldi, Physical Review Letters 118, 084301 (2017).
  * Moshe  _et al._ (2019) M. Moshe, E. Esposito, S. Shankar, B. Bircan, I. Cohen, D. R. Nelson, and M. J. Bowick, Physical Review Letters 122, 048001 (2019).
  * Sadik and Dias (2021) S. Sadik and M. A. Dias, Journal of the Mechanics and Physics of Solids 151, 104370 (2021).
  * Choi  _et al._ (2019) G. P. Choi, L. H. Dudte, and L. Mahadevan, Nature Materials 18, 999 (2019).
  * Choi  _et al._ (2021) G. P. Choi, L. H. Dudte, and L. Mahadevan, Physical Review Research 3, 043030 (2021).
  * Love (2013) A. E. H. Love,  _A treatise on the mathematical theory of elasticity_ (Cambridge university press, 2013).
  * Grosberg and Khokhlov (1994) A. Y. Grosberg and A. R. Khokhlov,  _Statistical physics of macromolecules_ (Amer Inst of Physics, 1994).
  * Inglis (1913) C. E. Inglis, Trans Inst Naval Archit 55, 219 (1913).
  * Zehnder and Potdar (1998) A. T. Zehnder and Y. K. Potdar, International Journal of Fracture 93, 409 (1998).
  * Han  _et al._ (2021) Q. Han, M. Lewicka, and L. Mahadevan, arXiv preprint (2021).
  * Yang  _et al._ (2021) Y. Yang, K. Vella, and D. P. Holmes, Science Robotics 6 (2021).
  * Weischedel  _et al._ (2012) C. Weischedel, A. Tuganov, T. Hermansson, J. Linn, and M. Wardetzky, (2012).
  * van Rees  _et al._ (2017) W. M. van Rees, E. Vouga, and L. Mahadevan, Proceedings of the National Academy of Sciences 114, 11597 (2017), https://www.pnas.org/content/114/44/11597.full.pdf .
  * Peyré (2011) G. Peyré, IEEE Computing in Science and Engineering 13, 94 (2011).
  * Kimmel and Sethian (1998) R. Kimmel and J. A. Sethian, Proceedings of the National Academy of Sciences 95, 8431 (1998).


* [14] Craig S Hamilton, Regina Kruse, Linda Sansoni, Sonja Barkhofen, Christine Silberhorn, and Igor Jex.  Gaussian boson sampling.  Physical review letters, 119(17):170501, 2017.
  * [15] Daniel J Brod, Ernesto F Galvão, Andrea Crespi, Roberto Osellame, Nicolò Spagnolo, and Fabio Sciarrino.  Photonic implementation of boson sampling: a review.  Advanced Photonics, 1(3):034001, 2019.
  * [16] Han-Sen Zhong, Hui Wang, Yu-Hao Deng, Ming-Cheng Chen, Li-Chao Peng, Yi-Han Luo, Jian Qin, Dian Wu, Xing Ding, Yi Hu, et al.  Quantum computational advantage using photons.  Science, 370(6523):1460–1463, 2020.
  * [17] Marlan O Scully and M Suhail Zubairy.  Quantum optics, 1999.
  * [18] Max Hofheinz, EM Weig, M Ansmann, Radoslaw C Bialczak, Erik Lucero, M Neeley, AD O’connell, H Wang, John M Martinis, and AN Cleland.  Generation of fock states in a superconducting quantum circuit.  Nature, 454(7202):310–314, 2008.
  * [19] H Wang, M Hofheinz, M Ansmann, RC Bialczak, E Lucero, M Neeley, AD O’connell, D Sank, J Wenner, AN Cleland, et al.  Measurement of the decay of fock states in a superconducting quantum circuit.  Physical Review Letters, 101(24):240401, 2008.
  * [20] Christopher S Wang, Jacob C Curtis, Brian J Lester, Yaxing Zhang, Yvonne Y Gao, Jessica Freeze, Victor S Batista, Patrick H Vaccaro, Isaac L Chuang, Luigi Frunzio, et al.  Efficient multiphoton sampling of molecular vibronic spectra on a superconducting bosonic processor.  Physical Review X, 10(2):021060, 2020.
  * [21] Gerhard Rempe, F Schmidt-Kaler, and Herbert Walther.  Observation of sub-poissonian photon statistics in a micromaser.  Physical review letters, 64(23):2783, 1990.
  * [22] Benjamin TH Varcoe, Simon Brattke, Matthias Weidinger, and Herbert Walther.  Preparing pure photon number states of the radiation field.  Nature, 403(6771):743–746, 2000.
  * [23] Clément Sayrin, Igor Dotsenko, Xingxing Zhou, Bruno Peaudecerf, Théo Rybarczyk, Sébastien Gleyzes, Pierre Rouchon, Mazyar Mirrahimi, Hadis Amini, Michel Brune, et al.  Real-time quantum feedback prepares and stabilizes photon number states.  Nature, 477(7362):73–77, 2011.
  * [24] Edo Waks, Eleni Diamanti, Barry C Sanders, Stephen D Bartlett, and Yoshihisa Yamamoto.  Direct observation of nonclassical photon statistics in parametric down-conversion.  Physical review letters, 92(11):113602, 2004.
  * [25] Merlin Cooper, Laura J Wright, Christoph Söller, and Brian J Smith.  Experimental generation of multi-photon fock states.  Optics express, 21(5):5309–5317, 2013.
  * [26] M Kitagawa and Y Yamamoto.  Number-phase minimum-uncertainty state with reduced number uncertainty in a kerr nonlinear interferometer.  Physical Review A, 34(5):3974, 1986.
  * [27] Daniel F Walls and Gerard J Milburn.  Quantum optics.  Springer Science & Business Media, 2007.
  * [28] PD Drummond and DF Walls.  Quantum theory of optical bistability. i. nonlinear polarisability model.  Journal of Physics A: Mathematical and General, 13(2):725, 1980\.
  * [29] Robert W Boyd.  Nonlinear optics.  Academic press, 2020.
  * [30] Hermann Haus.  Waves and fields in optoelectronics.  Prentice-Hall, Inc., 1984.
  * [31] Mehmet Fatih Yanik, Shanhui Fan, and Marin Soljačić.  High-contrast all-optical bistable switching in photonic crystal microcavities.  Applied Physics Letters, 83(14):2739–2741, 2003.
  * [32] Thomas Fink, Anne Schade, Sven Höfling, Christian Schneider, and Ataç Imamoglu.  Signatures of a dissipative phase transition in photon correlation measurements.  Nature Physics, 14(4):365–369, 2018.
  * [33] A Imamoḡlu, Helmut Schmidt, Gareth Woods, and Moshe Deutsch.  Strongly interacting photons in a nonlinear cavity.  Physical Review Letters, 79(8):1467, 1997.
  * [34] Aymeric Delteil, Thomas Fink, Anne Schade, Sven Höfling, Christian Schneider, and Ataç İmamoğlu.  Towards polariton blockade of confined exciton–polaritons.  Nature materials, 18(3):219–222, 2019.
  * [35] Nicholas Rivera and Ido Kaminer.  Light–matter interactions with photonic quasiparticles.  Nature Reviews Physics, 2(10):538–561, 2020.
  * [36] François Bondu, P Fritschel, CN Man, and A Brillet.  Ultrahigh-spectral-purity laser for the virgo experiment.  Optics letters, 21(8):582–584, 1996.
  * [37] Marlan O Scully and Willis E Lamb Jr.  Quantum theory of an optical maser. i. general theory.  Physical Review, 159(2):208, 1967.
  * [38] Shanhui Fan and John D Joannopoulos.  Analysis of guided resonances in photonic crystal slabs.  Physical Review B, 65(23):235112, 2002.
  * [39] Shanhui Fan, Wonjoo Suh, and John D Joannopoulos.  Temporal coupled-mode theory for the fano resonance in optical resonators.  JOSA A, 20(3):569–572, 2003.
  * [40] Melvin Lax.  Quantum noise vii: The rate equations and amplitude noise in lasers.  IEEE Journal of Quantum Electronics, 3(2):37–46, 1967.
  * [41] Y Yamamoto, S Machida, and O Nilsson.  Amplitude squeezing in a pump-noise-suppressed laser oscillator.  Physical Review A, 34(5):4025, 1986.
  * [42] S Machida and Y Yamamoto.  Observation of sub-poissonian photoelectron statistics in a negative feedback semiconductor laser.  Optics communications, 57(4):290–296, 1986.
  * [43] WH Richardson, S Machida, and Y Yamamoto.  Squeezed photon-number noise and sub-poissonian electrical partition noise in a semiconductor laser.  Physical review letters, 66(22):2867, 1991.
  * [44] Yoshihisa Yamamoto, Susumu Machida, and Wayne H Richardson.  Photon number squeezed states in semiconductor lasers.  Science, 255(5049):1219–1224, 1992.
  * [45] H Ritsch.  Quantum noise reduction in lasers with nonlinear absorbers.  Quantum Optics: Journal of the European Optical Society Part B, 2(2):189, 1990.
  * [46] DANIEL F Walls, MJ Collett, and AS Lane.  Amplitude-noise reduction in lasers with intracavity nonlinear elements.  Physical Review A, 42(7):4366, 1990.
  * [47] Seema Lathi and Yoshihisa Yamamoto.  Influence of nonlinear gain and loss on the intensity noise of a multimode semiconductor laser.  Physical Review A, 59(1):819, 1999.
  * [48] VSC Canela and HJ Carmichael.  Bright sub-poissonian light through intrinsic feedback and external control.  Physical review letters, 124(6):063604, 2020.
  * [49] HM Wiseman and GJ Milburn.  Noise reduction in a laser by nonlinear damping.  Physical Review A, 44(11):7815, 1991.
  * [50] H Ritsch, MAM Marte, and P Zoller.  Quantum noise reduction in raman lasers.  EPL (Europhysics Letters), 19(1):7, 1992.
  * [51] Amnon Yariv, Rashit Nabiev, and Kerry Vahala.  Self-quenching of fundamental phase and amplitude noise in semiconductor lasers with dispersive loss.  Optics letters, 15(23):1359–1361, 1990.
  * [52] J Kitching, R Boyd, A Yariv, and Y Shevy.  Amplitude noise reduction in semiconductor lasers with weak, dispersive optical feedback.  Optics letters, 19(17):1331–1333, 1994.
  * [53] Marin Soljačić and John D Joannopoulos.  Enhancement of nonlinear effects using photonic crystals.  Nature materials, 3(4):211–219, 2004.
  * [54] Kerry J Vahala.  Optical microcavities.  nature, 424(6950):839–846, 2003.
  * [55] Chia Wei Hsu, Bo Zhen, Jeongwon Lee, Song-Liang Chua, Steven G Johnson, John D Joannopoulos, and Marin Soljačić.  Observation of trapped light within the radiation continuum.  Nature, 499(7457):188–191, 2013.
  * [56] Chia Wei Hsu, Bo Zhen, A Douglas Stone, John D Joannopoulos, and Marin Soljačić.  Bound states in the continuum.  Nature Reviews Materials, 1(9):1–13, 2016.
  * [57] Mikhail F Limonov, Mikhail V Rybin, Alexander N Poddubny, and Yuri S Kivshar.  Fano resonances in photonics.  Nature Photonics, 11(9):543–554, 2017.
  * [58] Eric D Black.  An introduction to pound–drever–hall laser frequency stabilization.  American journal of physics, 69(1):79–87, 2001.
  * [59] Anton Frisk Kockum, Adam Miranowicz, Simone De Liberato, Salvatore Savasta, and Franco Nori.  Ultrastrong coupling between light and matter.  Nature Reviews Physics, 1(1):19–40, 2019.
  * [60] P Forn-Díaz, L Lamata, E Rico, J Kono, and E Solano.  Ultrastrong coupling regimes of light-matter interaction.  Reviews of Modern Physics, 91(2):025005, 2019.
  * [61] Philip Krantz, Morten Kjaergaard, Fei Yan, Terry P Orlando, Simon Gustavsson, and William D Oliver.  A quantum engineer’s guide to superconducting qubits.  Applied Physics Reviews, 6(2):021318, 2019.
  * [62] Ondřej Černotík, Aurélien Dantan, and Claudiu Genes.  Cavity quantum electrodynamics with frequency-dependent reflectors.  Physical review letters, 122(24):243601, 2019.
  * [63] Howard J Carmichael.  Statistical methods in quantum optics 2: Non-classical fields.  Springer Science & Business Media, 2009.
  * [64] Yoshihisa Yamamoto and Atac Imamoglu.  Mesoscopic quantum optics.  Mesoscopic Quantum Optics, 1999.
  * [65] Weng W Chow, Stephan W Koch, and Murray III Sargent.  Semiconductor-laser physics.  Springer Science & Business Media, 2012.
  * [66] Hermann Haken.  Waves, photons, atoms, volume 1.  North Holland, 1981.


# Solution of Definite Integrals using Functional Link Artificial Neural
Networks

Satyasaran Changdar  Department of Information Technology, Institute of
Engineering and Management, Salt Lake, Kolkata  Snehangshu Bhattacharjee
Department of Information Technology, Institute of Engineering and Management,
Salt Lake, Kolkata

###### Abstract

This paper discusses a new method to solve definite integrals using
feedforward artificial neural networks. The objective is to build a neural
network that would be a novel alternative to pre-established numerical methods
and with the help of a learning algorithm, be able to solve definite
integrals, by minimising a well constructed error function. The proposed
algorithm, with respect to existing numerical methods, is effective and
precise and well-suited for purposes which require integration of higher order
polynomials. The observations have been recorded and illustrated in tabular
and graphical form.

Keywords— Artificial neural networks, definite integral, Functional Link
Artificial Neural Network (FLANN), numerical methods

##  1 Introduction

Through centuries, men have been captivated by the power of the human brain.
Our very ability to capture, comprehend, analyse and draw conclusions from
seemingly innocent pieces of information, is a marvel which deserves
admiration. Only recently has science been able to gauge intelligence and its
parameters, and yet it is far from fathoming the intricate machinery of sense
and logic completely.
With several scientific attempts being made to study activities and patterns
of the brain since the nineteenth century, it was only in 1943 that McCulloch
and Pitts introduced a simple model of a neural network [16], which had
weighted inputs summed up to give the output. With the advent of computers in
the 1950s, it was finally possible to implement the theories. The concept of a
”thinking machine” has been much argued over, and the Dartmouth Project in
1956 followed by Frank Rosenblatt’s work in development of a Perceptron [20]
was a major impetus to the idea. Mark I Perceptron was built in the year 1958,
at Cornell University, with its main function being character recognition.
Soon after, Bernard Widrow and Marcian Hoff of Stanford successfully launched
neural network models ADALINE and MADALINE [23] for real-world purposes.
In 1982, John Hopfield of Caltech, with his charismatic paper presented to the
National Academy of Sciences [6], started rekindling interest in what was, by
far, an overlooked field in science. The US-Japan Joint Conference on
Cooperative/Competitive Neural Networks propelled efforts considerably and by
1985, the Neural Networks for Computing conference was held, which was to
become an annual affair. Also, IEEE’s first International Conference on Neural
Networks became a success in 1987. The last decade of the twentieth century
was crucial in bringing the study of neural networks to the fore with many
sectors in innovation and research recognising and targeting it.
ANN is being extensively used as a tool in various disciplines for
applications which include modelling complex mathematical structures and
obtaining solutions to mathematical problems. Mathematical models, with a
framework of ANN, can provide sustainable and effective solutions in
commercially applicable fields [7] [2] [3] [4]. Functions such as those
discussed below are frequently encountered in real-life problems of
engineering. Calculus is extensively used in manufacturing and other core
sectors for applications like improving architecture and determining the
infrastructure of bridges. It plays an important role in applications that
deal with massive data and require precision and calculus equations to make
predictions like in meteorology and in the finance sector. Even in the case of
search engines which demand a number of variables to be taken into account,
calculus comes into play. The next logical step would be incorporation of
neural networks in practical applications in the near future.
In the present paper, we attempt to discuss a way of solving Definite
Integrals by applying a direct, clinical and robust approach using the
effectiveness of neural networks instead of more cumbersome numerical methods
involving analytical calculations.

##  2 Motivation

Determining an Indefinite Integral refers to the process of determining a
function from its derivative _f(x)_ [22]. The functions that could possibly
have _f(x)_ as a derivative are the so-called antiderivatives of _f(x)_ ,and
the formula that leads to them is called the indefinite integral of _f(x)_.
Mathematically:

$\int{f(x)dx}=g(x)+C$ (1)

In this equation, _g(x)_ is one of the antiderivatives of _f(x)_ , and as the
constant _C_ changes, the equation gives a family of antiderivatives.
Let _g(x)_ be the antiderivative of the continuous function _f(x)_ defined on
_[a,b]_. Then the definite integral of _f(x)_ over _[a,b]_ is denoted by

$\int_{a}^{b}{f(x)dx}=g(b)-g(a)$ (2)

In the process of solving the above equation, we come across a number of
problems:
1\. The integrand _f(x)_ may be known only at certain points.
2\. A formula for the integrand _f(x)_ may be known, but it may be difficult
or impossible to find an antiderivative, that is, an elementary function.
3\. It may be possible to find an antiderivative symbolically, but it may be
easier to compute a numerical approximation than to compute the
antiderivative.
In the proposed method, we have discussed a solution in accordance to the
above points _(2)_ and _(3)_.
Usually, in cases where analytical solutions may not be arrived at easily, we
rely on numerical methods for solving the equation. This again leads to
specialised methods, specific to the purpose. A more appropriate approach
would be to design a solution that aptly fits in the places of these varied
methods and gives a decently precise and accurate solution. This leads us to
consider Artificial Neural Networks as an alternative [12] [13] [14]. The
efficiency of Artificial Neural Networks is due to its nonlinear learning
capability [8]. This powerful tool finds application in functional
approximation and optimisation. Due to application of FLANN (Functional Link
Artificial Neural Network) [15] [21], we no longer need hidden layers. This
makes the learning algorithm simpler. The FLANN architecture effectively
expands the input vector giving more flexibility and decision-making
capability to the network, due to a more expansive input space.

##  3 Mathematical modelling

We have already established that the definite integral of a function can be
calculated as shown in equation 2

$\int_{a}^{b}{f(x)dx}=g(b)-g(a)$ (3)

In the proposed model we would replace the _g(b)_ and _g(a)_ by _N(b)_ and
_N(a)_ respectively. Here, _N(x)_ represents a function of _x_ which is
obtained as the output of the neural network. Therefore our former equation
reduces to

$\int_{a}^{b}{f(x)dx}=N(b)-N(a)$ (4)

###  3.1 Model of the Neural Network

Output of the Neural Network

$N(x)=w_{1}\Phi_{1}(x)+w_{2}\Phi_{2}(x)+w_{3}\Phi_{3}(x)+......+w_{n}\Phi_{n}(x)$
(5)

In general, a model of a neural network comprises simple computational units
(nodes) which may be input nodes, output nodes or hidden nodes. The units,
with appropriate interconnections, is represented by a graph with the nodes as
the vertices and each edge between two neurons is labelled with a synaptic
weight. The network may or may not have multiple layers, depending on the
complexity of problems. As layers and number of neurons increase, the training
process becomes more complicated and requires intense computation. An
alternative to this cumbersome approach is the FLANN architecture
[15][18][17]. In our present study, we have employed the classification method
of Functional Link Artificial Neural Network (FLANN). FLANN architecture uses
a single layer feed-forward network, which can be considered as an alternative
to multi-layer neural network for solving problems. The drawback with multi-
layer neural network is its linear nature which becomes a hindrance for
mapping complex non-linear problems. Since the practical realm of data deals
mostly with non-linearity, therefore, choosing multi-layer networks for this
is not wise. Yet again, employing a multi-layer neural network for this
purpose further complicates the task. But in FLANN architecture, the input
vector is functionally expanded to overcome such problems, such that even
though it is a single layer neural network, it is well-suited to the task. The
model which has been used is shown in Fig.1.

Figure 1: Schematic diagram of the neural network model

If the input vector is _x_ , then after expansion, the input pattern is a set
of elements :

$\Phi_{1}(x),\Phi_{2}(x),.....\Phi_{n}(x)$ (6)

which are appropriate functions that will serve in the approximation process.
The set of functions for expansion may be set in accordance to the complexity
of the problem and the restraints posed on it.

###  3.2 Modelling the Error Function

As discussed above in equation 4

$\int_{a}^{b}{f(x)dx}=N(b)-N(a)$ (7)

Our objective is to build a function _N(x)_ using the neural network [10] [11]
[24] [9] which will approximate the function _g(x)_ upto a desired precision,
so that this function _N(x)_ can be later on employed to compute the values of
the integral at _a_ and _b_. Therefore,

$g(x)\approx N(x)$ (8)

In the ideal scenario, _g(x)_ equals _N(x)_

$g(x)-N(x)=0$ (9)

Differentiating equation 9 with respect to _x_ gives,

$g^{\prime}(x)-N^{\prime}(x)=0$ (10)

or,

$f(x)-N^{\prime}(x)=0$ (11)

where,

$N(x)=w_{1}\Phi_{1}(x)+w_{2}\Phi_{2}(x)+w_{3}\Phi_{3}(x)+......+w_{n}\Phi_{n}(x)$
(12)
$N^{\prime}(x)=w_{1}\Phi_{1}^{\prime}(x)+w_{2}\Phi_{2}^{\prime}(x)+w_{3}\Phi_{3}^{\prime}(x)+......+w_{n}\Phi_{n}^{\prime}(x)$
(13)

Therefore equation 11 becomes

$f(x)-(w_{1}\Phi_{1}^{\prime}(x)+w_{2}\Phi_{2}^{\prime}(x)+w_{3}\Phi_{3}^{\prime}(x)+......+w_{n}\Phi_{n}^{\prime}(x))=0$
(14)

The error function can now be summarised as:

$E(x)=\frac{1}{2}{\sum_{i=1}^{k}{{(f(x_{i})-N^{\prime}(x_{i}))}}^{2}}$ (15)

where $x_{i}$ refers to each training sample point; _k_ points belonging to
the interval _(a,b)_ have been considered.

###  3.3 Minimisation of Error

In this study, we have already constructed an error function, minimisation
thereof would lead to an appropriate solution to the given definite integral.
The error function uses a parameter which corresponds to the result derived
from the neural network. This parameter is modified with every modification of
the network (modification here corresponds to updation of the weights in the
network). The error optimisation algorithm used is gradient descent algorithm
[1] [19]. This algorithm uses first-order iterative optimisation to find the
minimum of a function. The algorithm is designed to proceed in steps
proportional to the gradient of the function, in the negative direction of the
current point. In the gradient descent algorithm, initially, the gradient of
the function is to be calculated. This is done by differentiating the error
function with respect to each of the weights one by one. These are then used
to update the weights. The process leads to reduction of the error gradually
until we reach a desired precision.

Differentiating 15 with respect to each weight:

$\frac{dE}{dw_{1}}=-\sum_{i=1}^{k}{e_{i}\Phi_{1}^{\prime}(x_{i})}$ (16)
$\frac{dE}{dw_{2}}=-\sum_{i=1}^{k}{e_{i}\Phi_{2}^{\prime}(x_{i})}$ (17)
$\frac{dE}{dw_{3}}=-\sum_{i=1}^{k}{e_{i}\Phi_{3}^{\prime}(x_{i})}$ (18)

and likewise, till

$\frac{dE}{dw_{n}}=-\sum_{i=1}^{k}{e_{i}\Phi_{n}^{\prime}(x_{i})}$ (19)

where,

$\emph{$e_{i}$}=f(x_{i})-N^{\prime}(x_{i})$ (20)

The weights are then updated as:

$w_{1}^{+}=w_{1}-\eta\frac{dE}{dw_{1}}$ (21)
$w_{2}^{+}=w_{2}-\eta\frac{dE}{dw_{2}}$ (22)
$w_{3}^{+}=w_{3}-\eta\frac{dE}{dw_{3}}$ (23)

and likewise, till

$w_{n}^{+}=w_{n}-\eta\frac{dE}{dw_{n}}$ (24)

where, $\eta$ is the _learning rate_.

##  4 Results and observations

Let us work out an example using the method discussed above. Consider the
function,

$f(x)=\sqrt{1+(x^{2})}$ (25)

We would like to integrate the function from 0 to 2, i.e

$\int_{0}^{2}{\sqrt{1+(x^{2})}dx}$ (26)

Using Maclaurin Series, a function _g(x)_ can be expanded as,

$g(x)=\sum_{n=0}^{\infty}{g^{(n)}(0)\frac{x^{n}}{n!}}$ (27)

where $g^{(n)}$ is the $n^{\text{th}}$ derivative of _g(x)_.
Expanding the concerned function using Maclaurin’s Series,

$\sqrt{(1+(x^{2}))}=1+\frac{x^{2}}{2}-\frac{x^{4}}{8}+\frac{x^{6}}{16}-.....$
(28)

Approximating equation 28 using the neural network would require a good choice
of the functional links $\phi_{1},\phi_{2},\phi_{3},....\phi_{n}$
In accordance to the Maclaurin Series expansion of _f(x)_ , which is a
polynomial function, we choose

$\phi_{1}(x)=x$ (29) $\phi_{2}(x)=x^{2}$ (30) $\phi_{3}(x)=x^{3}$ (31)

and likewise, till

$\phi_{n}(x)=x^{n}$ (32)

Therefore, equation 5 can be reconstructed as

$N(x)=w_{1}x+w_{2}x^{2}+w_{3}x^{3}+......+w_{n}x^{n}$ (33)

Similarly, equation 13 becomes

$N^{\prime}(x)=w_{1}+2w_{2}x+3w_{3}x^{2}+......+nw_{n}x^{n-1}$ (34)

The error function for the problem under consideration, becomes

$E(x)=\frac{1}{2}{\sum_{i=1}^{k}{(f(x)-(w_{1}+2w_{2}x+3w_{3}x^{2}+..+nw_{n}x^{n-1}))^{2}}}$
(35)

The objective is to optimise this error function for which we initialise the
network with random weights, a suitable learning rate $\eta$ is chosen, and
with every iteration, the weights are updated using the gradient descent
algorithm, as follows:

$w_{1}^{+}=w_{1}-\eta\frac{dE}{dw_{1}}$ (36)
$w_{1}^{+}=w_{1}-\eta\sum_{i=1}^{k}{-e_{i}}$ (37)
$w_{2}^{+}=w_{2}-\eta\frac{dE}{dw_{2}}$ (38)
$w_{2}^{+}=w_{2}-\eta\sum_{i=1}^{k}{-e_{i}(2x_{i})}$ (39)
$w_{3}^{+}=w_{3}-\eta\frac{dE}{dw_{3}}$ (40)
$w_{3}^{+}=w_{3}-\eta\sum_{i=1}^{k}{-e_{i}(3{x_{i}}^{2})}$ (41)

and likewise, till

$w_{n}^{+}=w_{n}-\eta\frac{dE}{dw_{n}}$ (42)
$w_{n}^{+}=w_{n}-\eta\sum_{i=1}^{k}{-e_{i}(n{x_{i}}^{n-1})}$ (43)

The above procedure is repeated till the error is minimised to a desired
precision. Finally, the weights are locked at the most favourable magnitudes
(where the error function is minimised). These weights, when now incorporated
in the network (i.e. $N(x)$ which has been approximated to $g(x)$, the assumed
anti-derivative of $f(x)$), provides us the result of the integral in any
interval within _[0,2]_ , i.e, the domain in which the network has been
trained. The results obtained from this neural network versus the exact values
calculated analytically, are illustrated in Fig.2

Figure 2: Graph comparing exact value with observed value using FLANN for
$I=\int_{0}^{b}{f(x)dx}$ where $f(x)=\sqrt{1+x^{2}}$, b$\in$(0,2]

On analysing the change in error with the number of iterations, we obtain
Fig.3

Figure 3: Error vs. Iteration while integrating $\sqrt{1+x^{2}}$

Similarly, the proposed algorithm is applied to integrate another function:

$\int_{0}^{2}{{2^{x}}dx}$ (44)

The Fig.4 is obtained in this case.

Figure 4: Graph comparing exact value with observed value using FLANN for
$I=\int_{0}^{b}{f(x)dx}$ where $f(x)=2^{x}$, b$\in$(0,2]

In order to illustrate the feasibility and efficiency of the proposed method,
a few examples are selected, the method is applied to them, and the results
are recorded in Fig.5. The given functions are integrated over [0,2].

Figure 5: Comparison of a few functions solved analytically and using the
discussed model

##  5 Examples for comparison with other methods

Conventional numerical methods, when implemented in calculating integrals, are
not always sufficient or accurate in approximation[5]. For instance, the
Trapezoidal rule gives an error:

$E^{T}_{n}(f)\equiv\int^{b}_{a}{f(x)dx}-T_{n}(f)\equiv-\frac{h^{2}(b-a)}{12}f^{2}(c_{n})$

for some $c_{n}$ in the interval (a,b).

As is evident, the error in Trapezoidal rule increases in the following
scenarios:
1\. If $h$ increases, where $h$ is the width of each subinterval between
$[a,b]$
2\. If $(b-a)$ is a large value
3\. If the function $f(x)$ is a polynomial of a high degree; which causes
$f^{2}(c_{n})$ to yield a large value.

Similarly, Simpson’s One-third rule gives an error:

$E^{S}(f)\equiv\int^{b}_{a}{f(x)dx}-S(f)\equiv-\frac{1}{180}(b-a)h^{4}f^{4}(c_{n})$

where $c_{n}$ lies in the interval (a,b).

The error, in this case increases for the following:
1\. If $h$ increases, where $h$ is the width of each subinterval between
$[a,b]$
2\. If $(b-a)$ is a large value
3\. If the function $f(x)$ is a polynomial of a very high degree, which causes
$f^{4}(c_{n})$ to yield a large value.

In the following example, we have considered the function

$\int^{b}_{a}{f(x)dx}=\int^{6}_{0}{x^{6}dx}$ (45)

Since the polynomial considered is of a high degree, therefore $f^{2}(c_{n})$
and $f^{4}(c_{n})$ yields a considerably larger value where $c_{n}$ lies in
the interval $(0,6)$. In cases where $f(x)$ is a higher degree polynomial and
$(b-a)$ is also large, both the Trapezoidal and Simpson’s Rule would fail to
come up with precise results.

To show the promising results of the proposed FLANN architecture, we have
compared the results of the equation 45 with those obtained from FLANN,
Trapezoidal method, Simpson’s rule and the same calculated analytically. While
training the neural network for the concerned problem, we have considered 10
training points _(i.e. k=10)_ belonging to the interval (a,b). Similarly,
using _k_ points in Trapezoidal and Simpson’s Rule implies that the interval
$[a,b]$ gets divided into _k+1_ sub-intervals, thus the width _h_ of each
subinterval becomes $\frac{b-a}{k+1}$. These results have been shown in Fig.6
and Fig.7.

Figure 6: Graph comparing the results of the integral $I=\int_{0}^{b}{f(x)dx}$
where $f(x)=x^{6}$, b$\in$(0,6] calculated using FLANN, Trapezoidal method,
Simpson’s rule and the same calculated analytically Figure 7: This graph is a
magnified version of Fig.6, which clearly shows the precision attained by the
proposed FLANN algorithm when compared to other techniques. Note that the
exact value has completely coincided with the calculated value using the FLANN
technique.

##  6 Practical Applications

While solving real life problems we come across various situations where
solving a definite integral analytically becomes difficult or seems
impossible. To tackle situations like these we have been using numerical
integration techniques. One such definite integral which is difficult to solve
analytically is the elliptic integral. Mathematics defines an elliptic
integral as any function f(x) which can be expressed in the form

$f(x)=\int_{c}^{x}{R(t,\sqrt{P(t)})dt}$ (46)

where R is a rational function of its two arguments, P is a polynomial of
degree 3 or 4 with no repeated roots, and c is a constant. In general,
integrals in this form cannot be expressed in terms of elementary functions
easily. This is where numerical integration techniques proves to be of immense
help and can often give precise results. Integrals like elliptic integrals
which are often difficult to solve have huge applications in astrophysics and
various engineering domains.

##  7 Conclusion

The main objective of the paper has been to establish a new method using
neural networks, which is applicable in solving definite integrals in a more
efficient and novel approach as compared to existing methods. We construct a
neural network model based on FLANN architecture, which is used to approximate
the integrand by training the neural network weights. The solved examples,
compared with other existing methods, show that the algorithm is effective and
feasible in all cases.
The proposed algorithm has the following advantages:
1\. The FLANN architecture used successfully incorporates non-linearity in the
model, which is a requirement in most practical scenarios.
2\. The neural network, when trained over the interval [a,b], can easily
estimate the values in another interval [$a_{1}$,$b_{1}$] where a $\leq$
$a_{1}$,$b_{1}$ $\leq$ b.
3\. Even when computing integrals of higher order polynomials, the method
evaluates to highly precise levels.

The successful implementation of the algorithm can depend on certain
improvisations, as follows:
1\. Change of optimisation method.
2\. Better choice of the functional links.
3\. Use of high-speed and more reliable processors.

## References

  * [1] Pierre Baldi.  Gradient descent learning algorithm overview: A general dynamical systems perspective.  IEEE Transactions on neural networks, 6(1):182–195, 1995.
  * [2] Satyasaran Changdar and Soumen De.  Analysis of non-linear pulsatile blood flow in artery through a generalized multiple stenosis.  Arabian Journal of Mathematics, 5(1):51–61, 2016.
  * [3] Satyasaran Changdar and Soumen De.  Analytical solution of mathematical model of magnetohydrodynamic blood nanofluid flowing through an inclined multiple stenosed artery.  Journal of Nanofluids, 6(6):1198–1205, 2017.
  * [4] Satyasaran Changdar and Soumen De.  Investigation of nanoparticle as a drug carrier suspended in a blood flowing through an inclined multiple stenosed artery.  Bionanoscience, 8(1):166–178, 2018.
  * [5] Joe D Hoffman and Steven Frankel.  Numerical methods for engineers and scientists.  CRC press, 2001.
  * [6] John J Hopfield.  Neural networks and physical systems with emergent collective computational abilities.  Proceedings of the national academy of sciences, 79(8):2554–2558, 1982.
  * [7] ALAN Horwitz and ARYA Ebrahimpour.  Engineering applications in differential and integral calculus.  International Journal of Engineering Education, 18(1):78–88, 2002\.
  * [8] Nicolaos Karayiannis and Anastasios N Venetsanopoulos.  Artificial neural networks: learning algorithms, performance evaluation, and applications, volume 209.  Springer Science & Business Media, 2013.
  * [9] Manoj Kumar and Neha Yadav.  Multilayer perceptrons and radial basis function neural network methods for the solution of differential equations: a survey.  Computers & Mathematics with Applications, 62(10):3796–3811, 2011\.
  * [10] Isaac E Lagaris, Aristidis Likas, and Dimitrios I Fotiadis.  Artificial neural networks for solving ordinary and partial differential equations.  IEEE transactions on neural networks, 9(5):987–1000, 1998.
  * [11] Isaac E Lagaris, Aristidis C Likas, and Dimitris G Papageorgiou.  Neural-network methods for boundary value problems with irregular boundaries.  IEEE Transactions on Neural Networks, 11(5):1041–1049, 2000.
  * [12] Susmita Mall and Snehashish Chakraverty.  Regression-based neural network training for the solution of ordinary differential equations.  International Journal of Mathematical Modelling and Numerical Optimisation, 4(2):136–149, 2013.
  * [13] Susmita Mall and Snehashish Chakraverty.  Chebyshev neural network based model for solving lane–emden type equations.  Applied Mathematics and Computation, 247:100–114, 2014.
  * [14] Susmita Mall and Snehashish Chakraverty.  Application of legendre neural network for solving ordinary differential equations.  Applied Soft Computing, 43:347–356, 2016.
  * [15] Susmita Mall and Snehashish Chakraverty.  Single layer chebyshev neural network model for solving elliptic partial differential equations.  Neural Processing Letters, 45(3):825–840, 2017.
  * [16] Warren S McCulloch and Walter Pitts.  A logical calculus of the ideas immanent in nervous activity.  The bulletin of mathematical biophysics, 5(4):115–133, 1943.
  * [17] Santosh Kumar Nanda and Debi Prasad Tripathy.  Application of functional link artificial neural network for prediction of machinery noise in opencast mines.  Advances in Fuzzy Systems, 2011:4, 2011.
  * [18] Santosh Kumar Nanda, Debi Prasad Tripathy, Simanta Kumar Nayak, and Subhasis Mohapatra.  Prediction of rainfall in india using artificial neural network (ann) models.  International Journal of Intelligent Systems and Applications, 5(12):1, 2013.
  * [19] Arunselvan Ramaswamy and Shalabh Bhatnagar.  Analysis of gradient descent methods with nondiminishing bounded errors.  IEEE Transactions on Automatic Control, 63(5):1465–1471, 2018.
  * [20] Frank Rosenblatt.  The perceptron: a probabilistic model for information storage and organization in the brain.  Psychological review, 65(6):386, 1958.
  * [21] Deepti Moyi Sahoo and Snehashish Chakraverty.  Functional link neural network learning for response prediction of tall shear buildings with respect to earthquake data.  IEEE Transactions on Systems, Man, and Cybernetics: Systems, 48(1):1–10, 2018.
  * [22] James Stewart.  Essential calculus: Early transcendentals.  Cengage Learning, 2012.
  * [23] Bernard Widrow and Michael A Lehr.  Artificial neural networks of the perceptron, madaline, and backpropagation family.  In Neurobionics, pages 133–205. Elsevier, 1993.
  * [24] Neha Yadav, Anupam Yadav, and Manoj Kumar.  Neural network methods for solving differential equations.  In An Introduction to Neural Network Methods for Differential Equations, pages 43–100. Springer, 2015.