# Scaling the tail beat frequency and swimming speed in underwater undulatory
swimming

Jesús Sánchez-Rodríguez1,2, Christophe Raufaste1,3 and Médéric Argentina1

( 1Université Côte d’Azur, CNRS, INPHYNI, France
2 Laboratory of Fluid Mechanics and Instabilities, École Polytechnique
Fédérale de Lausanne, CH-1015 Lausanne, Switzerland
3Institut Universitaire de France (IUF), France)

###### Abstract

Due to its great efficiency and maneuverability, undulatory swimming is the
predominant form of locomotion in aquatic vertebrates. A myriad of animals of
different species and sizes oscillate their bodies to propel themselves in
aquatic environments with swimming speed scaling as the product of the animal
length by the oscillation frequency. Although frequency tuning is the primary
means by which a swimmer selects its speed, there is no consensus on the
mechanisms involved. In this article, we propose scaling laws for undulatory
swimmers that relate oscillation frequency to length by taking into account
both the biological characteristics of the muscles and the interaction of the
moving swimmer with its environment. Results are supported by an extensive
literature review including approximately 1200 individuals of different
species, sizes and swimming environments. We highlight a crossover in length
around 0.5-1 m. Below this value, the frequency can be tuned between 2-20 Hz
due to biological constraints and the interplay between slow and fast muscles.
Above this value, the fluid-swimmer interaction must be taken into account and
the frequency is inversely proportional to the length of the animal. This
approach predicts a maximum swimming speed around 5-10 m.s-1 for large
swimmers, consistent with the threshold to prevent bubble cavitation.

##  1 Introduction

Beyond a few centimeters in length, most aquatic vertebrates propel themselves
through the water by deforming their spines and propagating deformation waves
through the body [1]. Fish, cetaceans, reptiles, amphibians, and birds
oscillate the head, body, tail, and/or fins, as appropriate, with a variety of
gaits described by a specific classification [2]. Despite the complexity of
treating each case separately, the kinematics of underwater undulatory
swimmers can be captured to first order with a few parameters such as the
wavelength of the deformation $\lambda$, the tail beat amplitude $A$ and the
tail beat frequency $f$. There is now considerable evidence that $\lambda$ and
$A$ are strongly related to animal length $L$, regardless of the size and
shape of the animal, or the swimming conditions. In the example of fish, the
wavelength scales as the animal length, with a factor of the order of unity
[3, 4]. The same applies to the tail beat amplitude that follows $A\simeq
0.2L$ [5, 6, 7, 8] from tadpoles of a few centimeters to whales of 20 meters
in length (Methods 4.1). These simple allometric scaling relations reveal
general physical laws that are valid over several orders of magnitude in size.
Momentum balance and minimal energy expenditure associated with the
hydrodynamic interaction between the moving body and the surrounding water
appear to drive the selection of the amplitude, as well as the determination
of the swimming speed, $U\simeq 3.3Af$ [5, 9, 10, 8]. Therefore, the swimming
speed is proportional to the oscillation frequency for a given swimmer and
scales as $Lf$, with a proportionality factor between 0.4 and 1 for fish and
cetaceans [11, 3, 12, 13, 14]. Taking the average factor, the relationship
$U\simeq 0.7Lf$ is therefore a very good approximation of the swimming speed
to within a factor of 2 at most. The swimming speed is thus intrinsically
linked to the oscillation frequency, but unlike the aforementioned scaling
laws that follow clear and widely documented trends over several orders of
magnitude in length, no consensus has been reached on the law that sets the
tail beat frequency. Most studies agree that the frequency decreases with the
length [5] and scaling laws $f\sim L^{-n}$ with an exponent $n$ ranging
between 0.5 and 1 are often reported together with models referring to
biological constraints, to the hydrodynamic interactions of the swimmer with
its environment or even to the effect of gravity [15, 16, 17, 18, 19]. The
difficulty of establishing a clear law and identifying the mechanisms at play
is related to the fact there is only a factor of 100 between the highest
frequency recorded in the smallest fish and the lowest frequency recorded in
the largest cetaceans, typically 20 and 0.2 Hz respectively (Fig. 1), while
measurements show a large dispersion for a given length. Firstly, several
parameters such as the previous training or the condition of the animal, its
age or sex, and the water or body temperature influence the final performance
of the swimmer [20, 21, 22]. Secondly, the frequency is not fixed for a given
swimmer but is the parameter that is adjusted to determine the swimming speed,
and a factor of 10 in frequency can easily be observed for a given specimen
[5].

Obtaining experimental results that maintain homogeneity in all these features
for a wide range of aquatic animal sizes is simply impossible and it is thus
not surprising to find diverse experimental laws in the literature. In
particular, comparing animals with the same level of activity, such as
sustained, prolonged or burst, would be necessary [23] but difficult to
implement in experiments [24, 25]. As a consequence, instead of focusing on a
specific activity gait, we propose to gather all data available in the
literature, regardless of the level of activity or any other specific
characteristic, to build a database of more than one thousand entries. This
approach provides a complete picture of the dependency of the frequency with
the length and captures both the main trend and the dispersion associated with
a given length. We thus propose a frequency selection mechanism that, on the
one hand, balances the swimmer’s muscle force and the reactive forces
generated by the fluid when the animal is in motion and, on the other hand,
considers the type of muscles, slow and fast. This works uncovers allometric
relations for the swimming frequency and speed that take the form of scaling
laws in the limits of very small or very large swimmers.

##  2 Results

###  2.1 Tail beat frequency measurements

We collected about 1200 data points from references listed in Methods 4.1,
with no discrimination on the basis of activity level or any of the other
parameters mentioned above, to avoid any possible bias in the length-frequency
relationship. In Fig. 1, it appears that the tail beat frequency is correlated
to the length, with all measurements located within a band in the $L-f$ plane.
For most of the lengths, the upper bound is well identified as the burst
activity level [26, 27, 28] and the frequency varies approximately by a factor
of 10 for a given $L$. For the longest animals, typically cetaceans with $L>5$
m, the magnitude of the band decreases. As we will discuss later, this
decrease is most likely associated with a lack of burst frequency
measurements: unlike the smaller animals, these animals were only observed in
their natural environment and were not forced to swim at peak activity. If we
define the upper and lower bounds of the band as fast and slow, we observe
that both follow the same behavior: the frequency is constant and maximum at
small lengths, typically $[f]_{\mathrm{fast}}\simeq 20$ Hz and
$[f]_{\mathrm{slow}}\simeq 2$ Hz respectively, before decreasing at larger
$L$. The change in tendency occurs around $L=0.5-1$ m. Since the span in $f$
observed for a given $L$ is of the same order of magnitude as the span in $f$
of a given specimen to adjust its speed, we conclude that the frequency
intervals are primarily associated with variations in the level of activity.

Figure 1: Tail beat frequency $f$ as a function of length $L$ for amphibians
(yellow), fish (blue), reptiles (green), birds (red) and mammals (purple).
Thick black and grey lines represent the burst and sustained activity levels,
respectively, fitted with the model. Thin lines are the scaling laws in the
limits of very small and very large swimmers. Values of the parameters for the
fast bound are $\left[L_{c}\right]_{\mathrm{fast}}=(0.5\pm 0.2)$ m,
$\left[f_{0}\right]_{\mathrm{fast}}=(21\pm 2)$ Hz,
$\left[\kappa\right]_{\mathrm{fast}}=(4\pm 3)$, while for the slow bound we
find $[L_{c}]_{\mathrm{slow}}=(1.0\pm 0.5)$ m,
$[f_{0}]_{\mathrm{slow}}=(2.0\pm 0.3)$ Hz, $[\kappa]_{\mathrm{slow}}=(9\pm
11)$.

In the next section, we interpret the data set by taking into account both the
muscle activity and the mechanical interaction of the swimmer with his
environment. Building on previous models, we account for limitations of muscle
in terms of both maximal force and frequency to explain the transition around
$L=0.5-1$ m between a regime characterized by a frequency that is constant and
a regime where the frequency decreases with length. Finally, we consider the
nature of the muscles, slow versus fast, to predict the slow and fast bounds
of the band in the $L-f$ plane.

###  2.2 Model and scaling laws

Locomotion occurs as a result of undulating movements produced by the
contraction of blocks of muscle segments [29]. While the muscles contract on
one side of the swimmer, those on the opposite side relax, alternately flexing
the entire body from side to side. Vertebrate muscles share a large number of
structural and functional features [3] that can be reasonably well described
by a limited number of parameters. This is the case for the relationship
between the force in the muscle, $F$, and the muscle contraction velocity,
$v$, as described by Hill’s muscle model [30]:

$F=F_{0}\frac{1-\displaystyle\frac{v}{v_{0}}}{1+\displaystyle\kappa\displaystyle\frac{v}{v_{0}}},$
(1)

where $F_{0}$ is the maximum isometric force generated in the muscle, obtained
as $v$ tends to zero, and $v_{0}$ is the maximum contraction velocity over
which no force can be produced. In its dimensionless form, the force $F/F_{0}$
is a decreasing and convex function of the velocity $v/v_{0}$, whose degree of
curvature is quantified by the parameter $\kappa$ (Methods 4.2). Since the
muscle fibers work in parallel, the maximum force scales as the cross area of
the muscle: $F_{0}\sim L^{2}\sigma_{0}$ where $\sigma_{0}$ is the maximum
isometric force per unit cross-sectional area [31]. The contraction velocity
$v$ drives the tail velocity and we can expect the scaling $v\sim Af$, which
leads to $v/v_{0}=f/f_{0}$, where $f_{0}$ is the maximum tail beat frequency
expected from the physiological limit of the muscle. $f_{0}$ can be inferred
by measuring the twitch contraction time $T$ of the muscle [32], that is, the
period of a single contraction and relaxation cycle produced by an action
potential within the muscle fibers, in order to deduce the maximum frequency
in the form $f_{0}=\frac{1}{2T}$ since one period consists of two antagonistic
contractions.

The swimmer’s movements set the surrounding fluid in motion to propel it
along. Given that the only force exerted by the animals on the fluid is
provided by the muscles, there must exist a balance between the forces
produced by the muscles and the reaction force of the fluid. Above a few
centimeters in length, aquatic organisms have a mode of locomotion based on
inertia [10]: while the body is oscillating, boluses of water of mass $\rho
L^{3}$ are set in motion with an acceleration $Af^{2}$ normal to the tail,
resulting in a lateral force that scales as $\rho L^{3}Af^{2}$, where $\rho$
is the density of water. As discussed in the Introduction and in the Method
4.1, $A$ is proportional to $L$ and thus the lateral force that pushes the
fluid scales as $\rho L^{4}f^{2}$. By balancing the latter with the force
exerted by the muscle (Eq. (1)), we obtain:

$\left(\frac{L}{{L_{c}}}\right)^{2}=\left(\frac{{f_{0}}}{f}\right)^{2}\frac{1-\displaystyle\frac{f}{f_{0}}}{1+\displaystyle\kappa\displaystyle\frac{f}{f_{0}}},$
(2)

where we have introduced the length $L_{c}$ that marks the crossover between
two regimes with different scaling laws:

$\displaystyle L_{c}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{\sigma_{0}}{\rho}}\frac{1}{f_{0}}$ (3)
$\displaystyle f$ $\displaystyle\simeq$ $\displaystyle
f_{0},\quad\mathrm{if}\quad L\ll L_{c}$ (4) $\displaystyle f$
$\displaystyle\simeq$ $\displaystyle
f_{0}\frac{L_{c}}{L},\quad\mathrm{if}\quad L\gg L_{c}.$ (5)

Eq. (2) predicts nicely the two bounds of the frequency band. In Fig. 1, we
have drawn the best fits (Methods 4.3) with the set of parameters
$[L_{c}]_{\mathrm{fast}}=0.5\pm 0.2$ m, $[f_{0}]_{\mathrm{fast}}=21\pm 2$ Hz
and $[\kappa]_{\mathrm{fast}}=4\pm 3$ for the fast bound and
$[L_{c}]_{\mathrm{slow}}=1.0\pm 0.5$ m, $[f_{0}]_{\mathrm{slow}}=2.0\pm 0.3$
Hz and $[\kappa]_{\mathrm{slow}}=9\pm 11$ for the slow bound. Unlike $f_{0}$
and $L_{c}$, $\kappa$ has large standard deviations from the fitted values; in
fact $\kappa$ weakly shapes the fit since it only plays a role in the
transition between the two limit regimes. These fits are also predictive: on
the example of humans and underwater undulatory swimming, also called dolphin
kick, the fit of the fast bound predicts a maximum kicking frequency around
3-4 Hz for a swimmer of about 2 meters, in good agreement with data recorded
in elite swimmers [33].

In addition, this model reconciles the two approaches and validates each in
its own length range.

For $L\ll L_{c}$, the frequency is fixed by a biological constraint, $f_{0}$,
in the spirit of the approach initiated by Wardle [32] who proposed that the
maximum tail beat frequency corresponds to the maximum frequency expected from
the muscles.

For $L\gg L_{c}$, the frequency decreases with length as a consequence of the
interplay between another biological constraint, here $\sigma_{0}$, and the
interaction of the swimmer with its environment [31, 16, 18]. In this limit,
the model predicts that the frequency is given by $f=cL^{-1}$, where we have
introduced the speed $c=L_{c}f_{0}=\sqrt{\sigma_{0}/\rho}$, that we estimate
to be $[c]_{\mathrm{slow}}\simeq 2.0\ \mathrm{m.s^{-1}}$ and
$[c]_{\mathrm{fast}}\simeq 10\ \mathrm{m.s^{-1}}$ for the slow and fast
bounds, respectively. The ratio of the wavelength of deformation $\lambda$ to
the body length ranges from 0.74 to 1.14 for anguillorm and thunniform
swimmers, respectively [3, 4], and $c\simeq\lambda f$ can also be interpreted
as the speed of the wave propagating along the swimmer’s body. Still in the
perspective of a mechanical approach, we remark here that the scaling $f\sim
L^{-1}$ from Eq. (5) is also consistent with an approximation of the swimmer
by an elastic beam of length $L$, radius $R$ and Young Modulus $E$. In such a
case, the natural frequency of bending waves scales as
$f_{e}\sim\frac{R}{L^{2}}\sqrt{\frac{E}{\rho}}$ [34]. By assuming geometrical
similarity, $R\sim L$, and body elasticity compatible with stresses generated
by muscle fibers $E\sim\sigma_{0}$, the natural bending frequency scales as
the tail beat frequency $f$ of large animals, defined in Eq. (5). This
suggests that large animals undulate their bodies to resonate with their
natural bending modes, as suggested by various studies [35, 36, 37].

Finally, we now understand why there are so many different exponents in the
literature resulting from attempts to describe the frequency-length
relationship as a scaling law (e.g., in [15, 16, 17, 18, 19]). Special
attention must be paid to the analysis of a range of lengths, because 1) the
scaling laws are only valid in either of the two limiting regimes and 2) the
data must be measured at the same activity level for the comparison to be
meaningful.

###  2.3 Considering the type of muscle fiber and comparison with biological
data

In order to connect the biology of swimmers to the fast and slow bounds in
Fig. 1, we take into account the different types of muscle fibers. Muscle
fibers can be roughly characterized as fast or slow, in part because the
latter type has a much lower level of ATP activity and a smaller contraction
speed but increased activity of oxidative enzymes. Therefore, slow fibers are
intended to produce forces over a prolonged period of activity [38], while
fast fibers use anaerobic chemical reactions and are adapted to rapid
movements, while producing higher forces. As a result, fast and slow muscles
are primarily solicited in burst and sustained activity levels, respectively
[39]. In what follows, we hypothesize that the presence of the lower frequency
boundary of the $f-L$ graph reflects the use of slow fibers, while to reach
high frequencies, fast fibers are exploited. The fast and slow frequency
bounds in Fig. 1 will therefore be modeled with the same Eq. (2), but with
parameters values that account for fast and slow muscles, respectively.

The parameters obtained from the fits match the biological data very well.
Wardle et al. suggests inferring $f_{0}$ from measurements of the twitch
contraction time [32], resulting in values between 5 and 25 Hz for 4 cm to 2.3
m fish [11]. These measurements are in agreement with the value
$[f_{0}]_{\mathrm{fast}}=21\ \mathrm{Hz}$ that fits the fast bound. Given that
$\sigma_{0}\simeq 200$ kN.m-2 for fast muscles is rather constant among
species [40], the estimate $L_{c}\sim 0.7$ m is in excellent agreement with
the fit $[L_{c}]_{\mathrm{fast}}=0.5\ \mathrm{m}$. For slow muscles, the
maximum frequency is smaller, and we employ the same approach to study the
slow bound, although measurements are rarer in this case. In the example of 10
cm salmon and 1 m sharks, $f_{0}$ ranges from 0.5–2 Hz, again in good
agreement with the $[f_{0}]_{\mathrm{slow}}=2.0\ \mathrm{Hz}$ obtained from
the fit of the slow bound. Measurements of $\sigma_{0}$ for slow muscles are
found between 20 and 80 kN.m-2 [41, 40, 42, 33]. If we take 50 kN.m-2 as the
typical value, we find $L_{c}\sim 3.5$ m, whose order of magnitude is coherent
with the value $[L_{c}]_{\mathrm{slow}}=1\ \mathrm{m}$ found from the fit.
Finally, the values of $\kappa$ adjusted for the fast and slow bounds,
$[\kappa]_{\mathrm{fast}}=4$ and $[\kappa]_{\mathrm{slow}}=9$, are in the
range of values recorded in vertebrates, typically between 2.5 and 10 [43,
44].

This framework with six parameters, three for each bound, appears coherent
because the filament arrangement in striated muscles is very similar along all
vertebrates [45].

In addition to the activity level, temperature has an effect on the twitch
contraction time, and therefore on $f_{0}$ [32, 46]. Cold water generally
depresses locomotory muscle function: the lower the temperature, the lower
$f_{0}$ is and a factor of 5 difference can be easily found between 2 and
30∘C. This can lead to deviations from the main trend at extreme temperatures.
This is probably the main explanation of the surprising low tail beat
frequency of Greenland sharks in Arctic waters ($L\simeq 3$ m and $f\simeq
0.15$ Hz in Fig. 1), a factor of 2–3 below the fit of the slow bound [17].
Conversely, there are examples of thermal acclimation of fish in warm waters
that could lead to $f_{0}$ values as high as 50 Hz [47].

###  2.4 Scaling the swimming speed

The relationship $U\simeq 0.7Lf$ intrinsically relates the tail beat frequency
to the swimming speed to a very good approximation, with a factor of 2 at most
for fish and cetaceans for $U$. This allows us to infer scaling laws for a
given activity level:

$\displaystyle U$ $\displaystyle\simeq$ $\displaystyle
0.7f_{0}L,\quad\mathrm{if}\quad L\ll L_{c}$ (6) $\displaystyle U$
$\displaystyle\simeq$ $\displaystyle 0.7f_{0}L_{c},\quad\mathrm{if}\quad L\gg
L_{c}.$ (7)

For small animals, we therefore expect swimming speed to increase with length,
whereas it should saturate at a constant value, $0.7\sqrt{\sigma_{0}/\rho}$
from Eqs. (3) and (6), for large animals. This is consistent with studies that
found a tail beat frequency scaling as $L^{-1}$ for large swimmers, and also
found a nearly constant swimming speed in this case [48, 18]. Still for large
swimmers, the swimming speed should range approximately between 1 and 8 m.s-1
depending on the activity level, given the values of $f_{0}$ and $L_{c}$ found
from the fits of the slow and fast bounds, respectively.

Figure 2:  Swimming speed $U$ as a function of length $L$. Following the law
$U=0.7Lf$, we show our estimates of swimming speed from the tail beat
frequency measurements displayed in Fig. 1 (closed circles): amphibians
(yellow), fish (blue), reptiles (green), birds (red), and mammals (purple).
Brown squares correspond to the data gathered by Hirt et al. for maximum
swimming speeds [49] using the mass-length relationship (Methods 4.1). Open
translucent squares represent either non-peer reviewed papers or data coming
from estimations and not actual measurements. Open opaque squares represent
data obtained using rod-mounted devices. The other data are represented by
closed opaque squares. The black and gray thick lines represent the fast and
slow boundaries predicted by the model and fitted parameters in Fig. 1,
respectively. Thin lines are the scaling laws in the limit of very small and
very large swimmers.

In order to properly describe the mechanisms at play, it is essential to
compare data at the same activity level. We therefore rely on the study by
Hirt et al., who reported data on maximum swimming speed [49]. In Fig. 2, we
superimpose these data on the expected swimming speeds inferred from all
frequency measurements. Various points need to be discussed.

First for $L\lesssim 0.5-1$ m, we observe that both data sets exhibit the same
upper limit. For most lengths, the match is perfect and only small differences
are observed, but they are at most a factor of 2. This regime is consistent
with the scaling law expected in the limit $L<L_{c}$: quantitatively, Hirt et
al. found that the maximum swimming speed scales as $U\propto M^{0.36}$, on
average, or equivalently $U\propto L$, in agreement with the scaling law
proposed in Eq. (4).

Second, for $L\gtrsim 0.5-1$ m, we observe a significant difference. While our
scaling law inferred from frequency measurements predicts a constant speed,
around 5–10 m.s-1, data gathered by Hirt et al. suggest a humped shape with a
maximum around 30–40 m.s-1 obtained for $L\simeq 1-3$ m followed by smaller
speeds for larger animals. In fact, we propose that the two data sets differ
in the two regions for two different reasons.

  1. 1.

While Hirt et al. did an enormous amount of work gathering data from the
literature, we suggest that the maximum is artificial if we apply relevant
filters. Most of the highest maximum speed data collected by Hirt et al., for
fish ranging from about 1 to 3 m in length, are estimates or predictions based
on in vitro physiological measurements. These measurements have been shown to
significantly overestimate the expected maximum speed of what were thought to
be the fastest swimmers, like billfish. First, these fish have lengths
$L>L_{c}$ and consequently their tail beat frequencies are significantly
smaller than the maximum frequency $f_{0}$ expected from the muscles: we
predict that a 2 m-long swimmer with $f_{0}=20\ \mathrm{Hz}$ would swim with a
maximum tail beat frequency four times smaller, below 5 Hz (Fig. 1), and thus
would have a maximum speed reduced by a factor of 4 in comparison to
predictions based on the twitch contraction time only. Second, recent
estimates based on measurements of twitch contraction times of anaerobic
muscles also provide upper bounds lower than $10\ \mathrm{m.s^{-1}}$ for
billfish and other large marine predatory fish [14]. Actually, reported values
of maximum speeds agree with this argument and refute some incredibly huge
values that had been estimated for animals of this size. For marlins, direct
measurements using speedometers showed that the observed maximum swimming
speed was around $2.25\ \mathrm{m.s^{-1}}$ [50], considerably lower than the
estimates of $30\ \mathrm{m.s^{-1}}$[51]. Burst speeds of sailfish also show
values around 8 m.s-1 measured with high-speed video and accelerometry [28], a
value much smaller than 30 m.s-1 [51]. In addition, there is some theoretical
evidence that the maximum speed should be smaller than 15 m.s-1, because
cavitation should appear at greater speeds, which should damage the flesh of
the swimmer [52]. In Fig. 2, we have used open translucent squares to
represent data that were not actual measurements but estimates, or that were
taken from non-peer-reviewed studies (Methods 4.4). If we remove these points
from the analysis, we find that the two upper bounds of the data sets match
very well, with the exception of four data points obtained for tuna and
barracuda that are still significantly faster than our fit of the fast bound
(open opaque squares in Fig. 2). Note that all of these points were measured
using rod-mounted devices that measure the speed at which the line is pulled
from the reel when a fish is hooked and pulling on the line ([5, 53, 54] and
Methods 4.4). Given the large fluctuations in the measurements made with this
method [53], it is likely that it overestimates the maximum speed, which is
also supported by the fact that barracuda and tuna do not show particularly
high maximum frequencies [14].

  2. 2.

For $L\gtrsim 5$ m, we attribute the discrepancy to a lack of tail beat
frequency measurements for very long swimmers, typically cetaceans, at a burst
level of activity. This would explain the jump in frequencies around $L\sim 5$
m in Fig. 1 for the fast bound. In this figure, the measured frequencies above
$L=5$ m correspond to swimming speeds between 1 and 4 m.s-1 [16, 19, 55].
Swimming speeds up to 10 m.s-1 were recorded in sperm whales [56], for which
tail beat frequencies were unfortunately not measured, but which should be
consistently higher than those plotted in Fig. 1 for the same length. Unlike
sperm whales, killer whales and some other large marine predators, most
cetaceans are filter feeders and do not have predators due to their size.
Therefore, they do not often need to move at maximum speeds, which would favor
data closer to the slow bound than the fast bound.

Following these considerations, it is reasonable to consider that the maximum
speed is constant for $L\gtrsim 0.5-1$ m, with a typical value around 5–10
m.s-1, in very good agreement with the prediction $0.7f_{0}L_{c}$ inferred
from the fit of the fast bound. From the definition of $L_{c}$, it means that
$U\sim\sqrt{\sigma_{0}/\rho}$ for $L>L_{c}$ and that swimming speed is
directly associated with the maximum stress generated by the muscles to push
the surrounding water. In their model, Hirt et al. state that heavier (and
consequently longer) animals need more time to accelerate to achieve maximum
speed and this fact would prevent the heaviest animals from being the fastest.
Here we suggest that the effect of a finite acceleration time would be a
second-order effect, unlike the other locomotion modes running and flying
[49].

##  3 Discussion

From the results, we conclude that length and activity level determine
swimming frequency and speed at the leading order. The data collected on
natural swimmers are explained using a simple model that accounts for
biological characteristics, through Hill’s muscle model, as well as the
interaction of the undulating swimmer with its environment. This model
requires only a few parameters that could be refined in the future to account
for specific characteristics of each swimmer (body temperature, swimming gait,
etc.). This work broadens our understanding of animal locomotion but should
also help in designing biomimetic and autonomous swimming robots [57, 58, 59,
60, 61], with the constraint that artificial swimmers have their own internal
characteristics that replace the biological ones discussed in the present
study.

Our study highlights a crossover at a length $L_{c}\sim 0.5-1$ m that
separates two limits: while small swimmers are constrained by biology only,
large swimmers are constrained by their environment as well. For a given
activity level, different scaling laws are found for swimming frequency and
speed in the two limits. This should also be the case for other quantities,
such as the muscle power for locomotion. For this quantity, our model predicts
scaling laws in $L^{5}$ and $L^{2}$ for very small and very large swimmers,
respectively (Methods 4.2). Measurements of oxygen consumption [24, 25] over a
wide range of lengths and activity levels might be a way to test these
predictions. In the framework of the model, $L_{c}$ also marks a significant
change in the way muscles are used. Small swimmers use muscles at their
maximum speed but negligible force in comparison to their maximal capabilities
(Methods 4.2). Large swimmers exhibit the opposite behavior. Given that muscle
power for locomotion scales as the product of force, frequency and length, we
expect muscle power for locomotion to be negligible in these two limits in
comparison to the maximum power available. Remarkably, only intermediate fish,
with lengths around $L_{c}$, would use the full capacity of muscle power to
undulate and move through water. In light of this comment, we can question in
the future whether intermediate fish are more likely to use economic
locomotion strategies (e.g., intermittent swimming [62], schooling [63], etc.)
compared to very small and large swimmers.

##  4 Methods

###  4.1 Data and allometry plots

We retrieved the length-frequency data available in the literature for a total
of 1202 animals, with a range of different species, morphologies and sizes. We
regrouped the data according to the classical division of vertebrates:
amphibians, fish, birds, reptiles and mammals. In figures, these data points
are displayed in yellow, blue, red, green and purple, respectively. In the
cases where the length data were not reported [16] but instead the mass of the
animals, the length was deduced through the allometric relation deduced in
Fig. 3b $L=0.44M^{0.33}$ ($L$ in meters and $M$ in kilograms), assuming
geometric similarity. Additionally, we recovered length-amplitude and length-
mass data. The references we used are provided in Tab. 1.

| References
---|---
| Frequency vs length or mass | Amplitude vs length or mass | Mass vs length
Amphibians | [64] | [64] | -
Fish | [5, 7, 65, 66, 26, 67, 68, 69, 70, 71, 72, 28, 17] | [5, 7, 65, 66, 26, 67, 72] | [5, 73, 74, 75, 76, 77]
Reptiles | [78] | [78] | -
Birds | [48, 79] | - | [48]
Mammals | [27, 80, 16, 81, 6, 55, 82] | [27, 80, 6, 82] | [83, 80, 84, 85, 27, 82]
Table 1: References reporting relations between length, frequency, amplitude
or mass of swimmers.

Our approach is based on two relationships verified through the data in Tab.
1. First, we plotted $A$ versus $L$ in Fig. 3a and verified that $A$ is
proportional to $L$, as originally shown by Bainbridge [5], but extended here
to 359 different specimens and four orders of magnitude in length. In fact the
best fit with a power law gives $A=(0.183\pm 0.002)L^{0.981\pm 0.005}$, which
shows that both quantities are proportional. The best proportionality relation
gives $A=(0.187\pm 0.001)L$. The geometric similarity for aquatic animals is
shown in Fig. 3b for 432 different individuals. The best fit of the data gives
$M=\left(12.14\pm 0.33\right)L^{3.04\pm 0.02}$, with $M$ in kilograms and $L$
in meters, which is consistent with the geometric similarity characterized by
an exponent of 3. This law extends Economos’s relation ($M=11.27L^{2.95}$
[86]) over four orders of magnitude in length or ten in mass. Forcing the
exponent of the relation to be exactly 3 gives $M=(12.09\pm 0.33)L^{3}$ or its
dimensional homogeneous form $M=(0.0121\pm 0.0003)\rho L^{3}$ with $\rho=1000$
kg.m-3.

a)  b)

Figure 3: a) Tail beat amplitude and b) animal mass as a function of animal
length. The solid lines represent the best power-law fits of the data.

###  4.2 Hill’s muscle model

Hill’s equations for tetanized muscle contraction (Eq. (1)) can be rewritten
to express the force in the muscle $F$ as a function of the swimming frequency
$f$:

$\frac{F}{F_{0}}=\frac{1-\frac{f}{f_{0}}}{1+\kappa\frac{f}{f_{0}}}.$ (8)

$F_{0}$ is the maximum isometric force generated in the muscle and $f_{0}$ is
the maximum tail beat frequency. In its dimensionless form, the force
$F/F_{0}$ is a decreasing and convex function of the frequency $f/f_{0}$,
whose degree of curvature is quantified by the parameter $\kappa$ (Fig. 4a).
From our analysis, very small animals ($L\ll L_{c}$) swim at maximum frequency
and negligible force ($f=f_{0}$ and $F\ll F_{0}$), while very large animals
($L\gg L_{c}$) swim at negligible frequency and maximum force ($f\ll f_{0}$
and $F=F_{0}$). Intermediate sized animals ($L\sim L_{c}$) swim at
intermediate frequency and force ($f\lessapprox f_{0}$ and $F\lessapprox
F_{0}$).

Figure 4: a) Muscle force and b) muscle power for locomotion as a function of
tail beat frequency, as predicted by Hill’s muscle model. The curves are drawn
with $\kappa=5$ and the quantities are represented in their dimensionless
form.

An estimate of muscle power for locomotion $P$ is obtained by multiplying the
muscle force by $Lf$, an estimate of the speed of muscle contraction.
$P\approx FLf$ can be expressed as a function of the tail beat frequency
following Hill’s equation:

$\frac{P}{F_{0}Lf_{0}}=\frac{1-\frac{f}{f_{0}}}{1+\kappa\frac{f}{f_{0}}}\frac{f}{f_{0}}.$
(9)

The power is drawn in Fig. 4b and $F_{0}Lf_{0}$ gives an estimate of the
maximum power that can be delivered by the muscle. This graph highlights the
fact that very small and large swimmers do not use the full capacity of muscle
power, unlike intermediate sized swimmers ($L\sim L_{c}$).

It is also possible to derive the expected scaling laws for muscle power as a
function of length. To do this, we take $F\sim\rho L^{4}f^{2}$ from the
interaction of the undulating swimmer with its environment and $f=f_{0}$ or
$\sqrt{\sigma_{0}/\rho}/L$ for very small or large swimmers, respectively.
This gives the following scaling laws in these two limits:

$\displaystyle P$ $\displaystyle\sim$ $\displaystyle\rho
f_{0}^{3}L^{5},\quad\mathrm{if}\quad L\ll L_{c}$ (10) $\displaystyle P$
$\displaystyle\sim$
$\displaystyle\rho^{-1/2}\sigma_{0}^{3/2}L^{2},\quad\mathrm{if}\quad L\gg
L_{c}.$ (11)

###  4.3 Characterization of the burst and sustained activity levels

All measurements are located within a band in the $L-f$ plane. The fast and
slow bounds corresponding to the burst and sustained activity levels,
respectively, are defined as follows. First, the length axis is divided into
$N$ intervals equally distributed in logarithmic scale. Second, in each
interval, the maximum and minimum frequencies are identified, as well as all
data points that are found between 90% and 100% of the maximum frequency value
and between 100% and 110% of the minimum frequency value. Finally, these two
sets of points are averaged independently to give one point for each bound
inside each interval.

For each bound, the three parameters of the model ($L_{c}$, $f_{0}$ and
$\kappa$) are obtained by fitting the $N$ averaged points (one per interval)
with Eq. (2) using the least absolute deviations (LAD) method. The robustness
and precision of the fit is probed by varying $N$ and taking the mean and
standard deviation to characterize each parameter. The number of intervals,
$N$, was varied between 10 and 50. This range ensured a statistically
significant number of points while at the same time avoiding empty intervals.
In practice, we selected nine values uniformly distributed between these two
extreme values.

###  4.4 Filtering of maximum speed data

We investigated the origin of the data gathered by Hirt et al. [49] to
establish objective criteria on data selection. Of all the measures reported
in the study, we found three classes of data that were not as reliable as the
others. First, we identified data coming from non-peer reviewed papers.
Second, some of the data are only estimates, not actual measurements. Third,
we found that all data obtained with rod-mounted devices are above the main
trend, which could be artificial due to the high fluctuations in these
measurements [53]. All these three classes of data points are summarized in
Tab. 2.

Animal | Mass (kg) | Length (m) |  Maximum speed (m.s-1) |  Speed reference |  Dismission reason
---|---|---|---|---|---
Aptenodytes patagonicus | 14 | 1.05 | 3.36 | [87] | Non-peer reviewed
Pygoscelis antarcticus | 4.5 | 0.72 | 8.89 | [87] | Non-peer reviewed
Pygoscelis papua | 8.2 | 0.88 | 10 | [87] | Non-peer reviewed
Acanthocybium solandri | 13.31 | 1.03 | 21.4 | [54] | Rod and reel
Acanthocybium solandri | 16.64 | 1.11 | 21.39 | [87] | Non-peer reviewed
Carcharodon carcharias | 800 | 4.04 | 11.11 | [87] | Non-peer reviewed
Galeocerdo cuvier | 550 | 3.57 | 8.89 | [87] | Non-peer reviewed
Istiompax indica | 150 | 2.32 | 36.11 | [50] | Estimation
Istiophorus albicans | 90 | 1.95 | 30.56 | [87] | Non-peer reviewed
Istiophorus albicans | 90 | 1.95 | 30 | [50] | Estimation
Isurus oxyrinchus | 105 | 2.06 | 18.8 | [88] | Non-peer reviewed
Isurus oxyrinchus | 300 | 2.92 | 13.89 | [87] | Non-peer reviewed
Makaira nigricans | 153.5 | 2.33 | 20.83 | [50] | Estimation
Sphyranea argentea | 4.5 | 0.72 | 12.22 | [87] | Non-peer reviewed
Sphyranea barracuda | 26.56 | 1.29 | 12.19 | [5] | Rod and reel
Tetrapturus audax | 163 | 2.38 | 22.5 | [87] | Non-peer reviewed
Thunnus albacares | 13.11 | 1.03 | 20.83 | [53] | Rod and reel
Thunnus orientalis | 250 | 2.74 | 19.44 | [87] | Non-peer reviewed
Thunnus thynnus | 27.22 | 1.31 | 19.67 | [5] | Rod and reel
Xiphias gladius | 98 | 2.01 | 26.94 | [87] | Non-peer reviewed
Balaenoptera musculus | 108400 | 20.75 | 10.30 | [31] | Estimation
Delphinus delphis | 95.32 | 1.99 | 10.30 | [31] | Estimation
Enhydra lutris | 30 | 1.35 | 2.5 | [87] | Non-peer reviewed
Megaptera novaeanglia | 30000 | 13.53 | 7.5 | [87] | Non-peer reviewed
Orcinus orca | 4100 | 6.97 | 13.33 | [87] | Non-peer reviewed
Orcinus orca | 4300 | 7.09 | 15.4 | [89] | Non-peer reviewed
Pusa hispida | 88.07 | 1.94 | 8.33 | [90] | Non-peer reviewed
Tursiops truncatus | 250 | 2.74 | 9.72 | [87] | Non-peer reviewed
Zalophus californianus | 158 | 2.36 | 11.11 | [91] | Non-peer reviewed
Dermochelys coriacea | 420 | 3.26 | 9.8 | [89] | Non-peer reviewed
Table 2: Data used in Hirt et al. [49] that are either provided by non-peer-
reviewed articles, given as estimates, or obtained using rod-mounted devices.
We list the type of animal together with the corresponding mass and speed, the
length using Methods 4.1 if not provided, and references.

#### Acknowledgements

We are grateful to François Gallaire, Guillaume Allibert and José Luis Trejo
for enlightning discussions. This work was supported by the French government,
through the UCAJEDI Investments in the Future project of the National Research
Agency (ANR-15-IDEX-01).

## References

  * [1] George V Lauder and Eric D Tytell.  Hydrodynamics of undulatory propulsion.  Fish physiology, 23:425–468, 2005.
  * [2] Michael Sfakiotakis, David M Lane, and J Bruce C Davies.  Review of fish swimming modes for aquatic locomotion.  IEEE Journal of oceanic engineering, 24(2):237–252, 1999.
  * [3] J. J. Videler.  Fish Swimming.  Springer Netherlands, 1993.
  * [4] Valentina Di Santo, Elsa Goerig, Dylan K. Wainwright, Otar Akanyeti, James C. Liao, Theodore Castro-Santos, and George V. Lauder.  Convergence of undulatory swimming kinematics across a diversity of fishes.  Proceedings of the National Academy of Sciences, 118(49):e2113206118, 2021.
  * [5] Richard Bainbridge.  The speed of swimming of fish as related to size and to the frequency and amplitude of the tail beat.  Journal of experimental biology, 35(1):109–133, 1958.
  * [6] Jim J Rohr and Frank E Fish.  Strouhal numbers and optimization of swimming by odontocete cetaceans.  Journal of Experimental Biology, 207(10):1633–1642, 2004.
  * [7] JR Hunter.  Swimming speed, tail beat frequency, tail beat amplitude and size in jack mackerel, trachurus symmetricus, and other fishes.  Fish. Bull., 69:253–266, 1971.
  * [8] M Saadat, Frank E Fish, AG Domel, V Di Santo, GV Lauder, and H Haj-Hariri.  On the rules for aquatic locomotion.  Physical Review Fluids, 2(8):083102, 2017.
  * [9] George S Triantafyllou, Michael S Triantafyllou, and Mark A Grosenbaugh.  Optimal thrust development in oscillating foils with application to fish propulsion.  Journal of Fluids and Structures, 7(2):205–224, 1993.
  * [10] Mattia Gazzola, Médéric Argentina, and Lakshminarayanan Mahadevan.  Scaling macroscopic aquatic locomotion.  Nature Physics, 10(10):758–761, 2014.
  * [11] J. J. Videler and C. S. Wardle.  Fish swimming stride by stride: speed limits and endurance.  Rev Fish Biol Fisheries, 1(1):23–40, September 1991.
  * [12] J. Videler and P. Kamermans.  Differences between upstroke and downstroke in swimming dolphins.  Journal of Experimental Biology, 119(1):265–274, November 1985\.
  * [13] Kristina Charlotte Curren.  Designs for swimming : morphometrics and swimming dynamics of several cetacean species.  Master’s thesis, Memorial University of Newfoundland, 1992.
  * [14] Morten B. S. Svendsen, Paolo Domenici, Stefano Marras, Jens Krause, Kevin M. Boswell, Ivan Rodriguez-Pinto, Alexander D. M. Wilson, Ralf H. J. M. Kurvers, Paul E. Viblanc, Jean S. Finger, and John F. Steffensen.  Maximum swimming speeds of sailfish and three other large marine predatory fish species based on muscle contraction time and stride length: a myth revisited.  Biology Open, 5(10):1415–1419, August 2016.
  * [15] Adrian Bejan and James H Marden.  Unifying constructal theory for scale effects in running, swimming and flying.  Journal of Experimental Biology, 209(2):238–248, 2006.
  * [16] Katsufumi Sato, Yutaka Watanuki, Akinori Takahashi, Patrick JO Miller, Hideji Tanaka, Ryo Kawabe, Paul J Ponganis, Yves Handrich, Tomonari Akamatsu, Yuuki Watanabe, et al.  Stroke frequency, but not swimming speed, is related to body size in free-ranging seabirds, pinnipeds and cetaceans.  Proceedings of the Royal Society B: Biological Sciences, 274(1609):471–477, 2007.
  * [17] Yuuki Y. Watanabe, Christian Lydersen, Aaron T. Fisk, and Kit M. Kovacs.  The slowest fish: Swim speed and tail-beat frequency of Greenland sharks.  Journal of Experimental Marine Biology and Ecology, 426-427:5–11, September 2012.
  * [18] Rahul Bale, Max Hao, Amneet Pal Singh Bhalla, and Neelesh A Patankar.  Energy efficiency and allometry of movement of swimming and flying animals.  Proceedings of the National Academy of Sciences, 111(21):7517–7521, 2014.
  * [19] William T Gough, Paolo S Segre, KC Bierlich, David E Cade, Jean Potvin, Frank E Fish, Julian Dale, Jacopo di Clemente, Ari S Friedlaender, David W Johnston, et al.  Scaling of swimming performance in baleen whales.  Journal of Experimental Biology, 222(20):jeb204172, 2019.
  * [20] Wen-Wen Zhao, Xu Pang, Jiang-Lan Peng, Zhen-Dong Cao, and Shi-Jian Fu.  The effects of hypoxia acclimation, exercise training and fasting on swimming performance in juvenile qingbo (spinibarbus sinensis).  Fish physiology and biochemistry, 38(5):1367–1377, 2012.
  * [21] Peter J Allen, Brian Hodge, Inge Werner, and Joseph J Cech, Jr.  Effects of ontogeny, season, and temperature on the swimming performance of juvenile green sturgeon (acipenser medirostris).  Canadian Journal of Fisheries and Aquatic Sciences, 63(6):1360–1369, 2006.
  * [22] Yuya Makiguchi, Hisaya Nii, Katsuya Nakao, and Hiroshi Ueda.  Sex differences in metabolic rate and swimming performance in pink salmon (o ncorhynchus gorbuscha): the effect of male secondary sexual traits.  Ecology of Freshwater Fish, 26(2):322–332, 2017.
  * [23] Theodore Y Wu.  Introduction to the scaling of aquatic animal locomotion.  Technical report, CALIFORNIA INST OF TECH PASADENA, 1977.
  * [24] J Roland Brett.  The respiratory metabolism and swimming performance of young sockeye salmon.  Journal of the Fisheries Board of Canada, 21(5):1183–1226, 1964\.
  * [25] J. R. Brett.  The metabolic demand for oxygen in fish, particularly salmonids, and a comparison with other vertebrates.  Respiration Physiology, 14(1):151–170, 1972.
  * [26] JJ Videler and F Hess.  Fast continuous swimming of two pelagic predators, saithe (Pollachius virens) and mackerel (Scomber scombrus): a kinematic analysis.  Journal of experimental biology, 109(1):209–228, 1984.
  * [27] Frank E Fish.  Comparative kinematics and hydrodynamics of odontocete cetaceans: morphological and ecological correlates with swimming performance.  The Journal of experimental biology, 201(20):2867–2877, 1998.
  * [28] Stefano Marras, Takuji Noda, John F Steffensen, Morten BS Svendsen, Jens Krause, Alexander DM Wilson, Ralf HJM Kurvers, James Herbert-Read, Kevin M Boswell, and Paolo Domenici.  Not so fast: swimming behavior of sailfish during predator–prey interactions using high-speed video and accelerometry.  Integrative and comparative biology, 55(4):719–727, 2015.
  * [29] Robert E Shadwick and George V Lauder.  Fish physiology: fish biomechanics.  Elsevier, 2006.
  * [30] Archibald Vivian Hill.  The heat of shortening and the dynamic constants of muscle.  Proceedings of the Royal Society of London. Series B-Biological Sciences, 126(843):136–195, 1938.
  * [31] A. V. HILL.  THE DIMENSIONS OF ANIMALS AND THEIR MUSCULAR DYNAMICS.  Science Progress (1933- ), 38(150):209–230, 1950.
  * [32] C. S. Wardle.  Limit of fish swimming speed.  Nature, 255(5511):725–727, June 1975.
  * [33] Yuiko Matsuura, Naoto Matsunaga, Satoshi Iizuka, Hiroshi Akuzawa, and Koji Kaneoka.  Muscle Synergy of the Underwater Undulatory Swimming in Elite Male Swimmers.  Frontiers in Sports and Active Living, 2:62, June 2020.
  * [34] Lev Davidovich Landau, Evgenij M Lifšic, Evegnii Mikhailovich Lifshitz, Arnold Markovich Kosevich, and Lev Petrovich Pitaevskii.  Theory of elasticity: volume 7, volume 7.  Elsevier, 1986.
  * [35] Mattia Gazzola, Médéric Argentina, and Lakshminarayanan Mahadevan.  Gait and speed selection in slender inertial swimmers.  Proceedings of the National Academy of Sciences, 112(13):3874–3879, 2015.
  * [36] Florine Paraz, Lionel Schouveiler, and Christophe Eloy.  Thrust generation by a heaving flexible foil: Resonance, nonlinearities, and optimality.  Physics of Fluids, 28(1):011903, 2016.
  * [37] Alexander P Hoover, Ricardo Cortez, Eric D Tytell, and Lisa J Fauci.  Swimming performance, resonance and shape evolution in heaving flexible panels.  Journal of Fluid Mechanics, 847:386–416, 2018.
  * [38] Robert H. Fitts, Kerry S. McDonald, and Jane M. Schluter.  The determinants of skeletal muscle force and power: Their adaptability with changes in activity pattern.  Journal of Biomechanics, 24(SUPPL. 1):111–122, 1991.  Cited by: 120.
  * [39] Lawrence C Rome, Roel P Funke, R McNeill Alexander, Gordon Lutz, Hugh Aldridge, Frank Scott, and Marvin Freadman.  Why animals have different muscle fibre types.  Nature, 335(6193):824–827, 1988.
  * [40] Ian A. Johnston and Richard Brill.  Thermal dependence of contractile properties of single skinned muscle fibres from Antarctic and various warm water marine fishes including Skipjack Tuna (Katsuwonus pelamis) and Kawakawa (Euthynnus affinis).  J Comp Physiol B, 155(1):63–70, January 1984.
  * [41] J. D. Altringham and I. A. Johnston.  The pCa-tension and force-velocity characteristics of skinned fibres isolated from fish fast and slow muscles.  The Journal of Physiology, 333(1):421–449, 1982.
  * [42] IA Johnston, Bruce D Sidell, and WR Driedzic.  Force-velocity characteristics and metabolism of carp muscle fibres following temperature acclimation.  Journal of Experimental Biology, 119(1):239–249, 1985.
  * [43] DR Wilkie.  The relation between force and velocity in human muscle.  The Journal of physiology, 110(3-4):249, 1949.
  * [44] IA Johnston.  Sustained force development: specializations and variation among the vertebrates.  Journal of Experimental Biology, 115(1):239–251, 1985.
  * [45] Scott L Hooper and Jeffrey B Thuma.  Invertebrate muscles: muscle specific genes and proteins.  Physiological reviews, 2005.
  * [46] C. S. Wardle.  Effects of Temperature on the Maximum Swimming Speed of Fishes.  In M. A. Ali, editor, Environmental Physiology of Fishes, NATO Advanced Study Institutes Series, pages 519–531. Springer US, Boston, MA, 1980.
  * [47] T Johnson and A Bennett.  The thermal acclimation of burst escape performance in fish: an integrated study of molecular and cellular physiology and organismal performance.  Journal of Experimental Biology, 198(10):2165–2175, October 1995\.
  * [48] Katsufumi Sato, Kozue Shiomi, Yuuki Watanabe, Yutaka Watanuki, Akinori Takahashi, and Paul J Ponganis.  Scaling of swim speed and stroke frequency in geometrically similar penguins: they swim optimally to minimize cost of transport.  Proceedings of the Royal Society B: Biological Sciences, 277(1682):707–714, 2010.
  * [49] Myriam R Hirt, Walter Jetz, Björn C Rall, and Ulrich Brose.  A general scaling law reveals why the largest animals are not the fastest.  Nature Ecology & Evolution, 1(8):1116–1122, 2017.
  * [50] Barbara A Block, David Booth, and Francis G Carey.  Direct measurement of swimming speeds and depth of blue marlin.  Journal of Experimental Biology, 166(1):267–284, 1992.
  * [51] FW Lane.  How fast do fish swim.  Country Life (London), 90:534–535, 1941.
  * [52] G Iosilevskii and D Weihs.  Speed limits on swimming of fishes and cetaceans.  Journal of The Royal Society Interface, 5(20):329–338, 2008.
  * [53] Vladimir Walters and Harry L Fierstine.  Measurements of swimming speeds of yellowfin tuna and wahoo.  Nature, 202(4928):208, 1964.
  * [54] John Tyler Bonner.  Size and cycle: an essay on the structure of biology, volume 2087\.  Princeton University Press, 2015.
  * [55] William T Gough, Hayden J Smith, Matthew S Savoca, Max F Czapanskiy, Frank E Fish, Jean Potvin, KC Bierlich, David E Cade, Jacopo Di Clemente, John Kennedy, et al.  Scaling of oscillatory kinematics and froude efficiency in baleen whales.  Journal of Experimental Biology, 224(13):jeb237586, 2021.
  * [56] Kagari Aoki, Masao Amano, Naoki Sugiyama, Hiroyuki Muramoto, Michihiko Suzuki, Motoi Yoshioka, Kyoichi Mori, Daisuke Tokuda, and Nobuyuki Miyazaki.  Measurement of swimming speed in sperm whales.  In 2007 Symposium on Underwater Technology and Workshop on Scientific Use of Submarine Cables and Related Technologies, pages 467–471. IEEE, 2007.
  * [57] Vincent Lebastard, Frédéric Boyer, and Sylvain Lanneau.  Reactive underwater object inspection based on artificial electric sense.  Bioinspir. Biomim., 11(4):045003, 2016.
  * [58] J. Zhu, C. White, D. K. Wainwright, V. Di Santo, G. V. Lauder, and H. Bart-Smith.  Tuna robotics: A high-frequency experimental platform exploring the performance space of swimming fishes.  Science Robotics, 4(34):eaax4615, 2019.
  * [59] J. Sánchez-Rodríguez, F. Celestini, C. Raufaste, and M. Argentina.  Proprioceptive Mechanism for Bioinspired Fish Swimming.  Phys. Rev. Lett., 126(23):234501, 2021.
  * [60] Robin Thandiackal, Kamilo Melo, Laura Paez, Johann Herault, Takeshi Kano, Kyoichi Akiyama, Frédéric Boyer, Dimitri Ryczko, Akio Ishiguro, and Auke J. Ijspeert.  Emergence of robust self-organized undulatory swimming based on local hydrodynamic force sensing.  Science Robotics, 6(57):eabf6354, August 2021.
  * [61] Keel Yong Lee, Sung-Jin Park, David G. Matthews, Sean L. Kim, Carlos Antonio Marquez, John F. Zimmerman, Herdeline Ann M. Ardoña, Andre G. Kleber, George V. Lauder, and Kevin Kit Parker.  An autonomously swimming biohybrid fish designed with human cardiac biophysics.  Science, 375(6581):639–647, 2022.
  * [62] Gen Li, Intesaaf Ashraf, Bill François, Dmitry Kolomenskiy, Frédéric Lechenault, Ramiro Godoy-Diana, and Benjamin Thiria.  Burst-and-coast swimmers optimize gait by adapting unique intrinsic cycle.  Communications biology, 4(1):1–7, 2021.
  * [63] Stefano Marras, Shaun S. Killen, Jan Lindström, David J. McKenzie, John F. Steffensen, and Paolo Domenici.  Fish swimming in schools save energy regardless of their spatial position.  Behav Ecol Sociobiol, 69(2):219–226, 2015.
  * [64] RICHARD J WASSERSUG and KARIN VON SECHENDORF HOFF.  The kinematics of swimming in anuran larvae.  Journal of experimental Biology, 119(1):1–30, 1985.
  * [65] Paul W Webb.  Kinematics of lake sturgeon, Acipenser fulvescens, at cruising speeds.  Canadian Journal of Zoology, 64(10):2137–2141, 1986.
  * [66] Paul W Webb.  Steady swimming kinematics of tiger musky, an esociform accelerator, and rainbow trout, a generalist cruiser.  Journal of Experimental Biology, 138(1):51–69, 1988.
  * [67] Lisa J Rosenberger and Mark W Westneat.  Functional morphology of undulatory pectoral fin locomotion in the stingray Taeniura lymma (Chondrichthyes: Dasyatidae).  Journal of Experimental Biology, 202(24):3523–3539, 1999.
  * [68] John J Magnuson and John H Prescott.  Courtship, locomotion, feeding, and miscellaneous behaviour of pacific bonito (sarda chiliensis).  Animal behaviour, 14(1):54–67, 1966.
  * [69] Anna-Maria Mueller, Deborah L Burwen, Kevin M Boswell, and Tim Mulligan.  Tail-beat patterns in dual-frequency identification sonar echograms and their potential use for species identification and bioenergetics studies.  Transactions of the American Fisheries Society, 139(3):900–910, 2010.
  * [70] JD Thiem, JW Dawson, AC Gleiss, EG Martins, A Haro, T Castro-Santos, AJ Danylchuk, RP Wilson, and SJ Cooke.  Accelerometer-derived activity correlates with volitional swimming speed in lake sturgeon (acipenser fulvescens).  Canadian Journal of Zoology, 93(8):645–654, 2015.
  * [71] PW Webb and Raymond S Keyes.  Swimming kinematics of sharks.  Fishery Bulletin, 80(4):803–812, 1982.
  * [72] Jeffrey B Graham, Heidi Dewar, NC Lai, William R Lowell, and Steve M Arce.  Aspects of shark swimming performance determined using a large water tunnel.  Journal of Experimental Biology, 151(1):175–192, 1990.
  * [73] G-J Yan, X-K He, Z-D Cao, and S-J Fu.  An interspecific comparison between morphology and swimming performance in cyprinids.  Journal of Evolutionary Biology, 26(8):1802–1815, 2013.
  * [74] Arif Malik, Kathryn A Dickson, Takashi Kitagawa, Ko Fujioka, Ethan E Estess, Charles Farwell, Kristy Forsgren, Jeannette Bush, and Kathryn A Schuller.  Ontogeny of regional endothermy in pacific bluefin tuna (thunnus orientalis).  Marine Biology, 167(9):1–20, 2020.
  * [75] DCT Dissanayake, EKV Samaraweera, and C Amarasiri.  Fishery and feeding habits of yellowfin tuna (thunnus albacares) targeted by coastal tuna longlining in the north western and north eastern coasts of sri lanka.  Sri Lanka J Aquat Sci, 13:1–21, 2008.
  * [76] Sahar F Mehanna and Alam Eldeen Farouk.  Length-weight relationship of 60 fish species from the eastern mediterranean sea, egypt (gfcm-gsa 26).  Frontiers in Marine Science, page 942, 2021.
  * [77] FS Motta, FP Caltabellotta, RC Namora, and OBF Gadig.  Technical contribution length-weight relationships of sharks caught by artisanal fisheries from southeastern brazil.  J. Appl. Ichthyol, 1:2, 2013.
  * [78] Frank E Fish.  Kinematics of undulatory swimming in the american alligator.  Copeia, pages 839–843, 1984.
  * [79] Brian D Clark and Willy Bemis.  Kinematics of swimming of penguins at the detroit zoo.  Journal of Zoology, 188(3):411–428, 1979.
  * [80] Tricia Kojeszewski and Frank E Fish.  Swimming kinematics of the florida manatee (trichechus manatus latirostris): hydrodynamic analysis of an undulatory mammalian swimmer.  Journal of Experimental Biology, 210(14):2411–2418, 2007.
  * [81] Jeremy A Goldbogen, John Calambokidis, Robert E Shadwick, Erin M Oleson, Mark A McDonald, and John A Hildebrand.  Kinematics of foraging dives and lunge-feeding in fin whales.  Journal of Experimental Biology, 209(7):1231–1244, 2006.
  * [82] Frank E Fish, S Innes, and K Ronald.  Kinematics and estimated thrust production of swimming harp and ringed seals.  Journal of Experimental Biology, 137(1):157–173, 1988.
  * [83] Masaharu Nishiwaki.  On the body weight of whales.  Scientific Reports of the Whales Research Institute, 4:184–209, 1950.
  * [84] Todd R Robeck, Steven L Monfort, Paul P Calle, J Lawrence Dunn, Eric Jensen, Jeffrey R Boehm, Skip Young, and Steven T Clark.  Reproduction, growth and development in captive beluga (delphinapterus leucas).  Zoo Biology: Published in affiliation with the American Zoo and Aquarium Association, 24(1):29–49, 2005.
  * [85] Fredrik Christiansen, Mariano Sironi, Michael J Moore, Matías Di Martino, Marcos Ricciardi, Hunter A Warick, Duncan J Irschick, Robert Gutierrez, and Marcela M Uhart.  Estimating body mass of free-living whales using aerial photogrammetry and 3d volumetrics.  Methods in Ecology and Evolution, 10(12):2034–2044, 2019.
  * [86] AC Economos.  Elastic and/or geometric similarity in mammalian design?  Journal of Theoretical Biology, 103(1):167–172, 1983.
  * [87] https://www.speedofanimals.com.
  * [88] Wikipedia contributors.  Shortfin mako shark — Wikipedia, the free encyclopedia.  https://en.wikipedia.org/w/index.php?title=Shortfin_mako_shark&oldid=1118544271, 2022\.
  * [89] www.livescience.com/32772-what-animal-is-the-fastest-swimmer.html.
  * [90] www.oceanwide-expeditions.com/to-do/wildlife/ringed-seal-1.
  * [91] www.nationalgeographic.com/animals/mammals/c/california-sea-lion.


11institutetext: Gentian Zavalani 22institutetext: Center for Advanced Systems
Understanding (CASUS), Görlitz, Germany 33institutetext: Technical University
Dresden, Faculty of Mathematics, Dresden, Germany 33email<EMAIL_ADDRESS>44institutetext: Michael Hecht 55institutetext: Center for Advanced Systems
Understanding (CASUS), Görlitz, Germany 66institutetext: University of
Wroclaw, Mathematical Institute, Wroclaw, Poland 66email<EMAIL_ADDRESS>
# High-order numerical integration on regular embedded surfaces

Gentian Zavalani and Michael Hecht

###### Abstract

We present a high-order surface quadrature (HOSQ) for accurately approximating
regular surface integrals on closed surfaces. The initial step of our approach
rests on exploiting _square-squeezing_ –a homeomorphic bilinear square-simplex
transformation, re-parametrizing any surface triangulation to a quadrilateral
mesh. For each resulting quadrilateral domain we interpolate the geometry by
tensor polynomials in Chebyshev–Lobatto grids. Posterior the tensor-product
Clenshaw-Curtis quadrature is applied to compute the resulting integral. We
demonstrate efficiency, fast runtime performance, high-order accuracy, and
robustness for complex geometries.

###### Keywords:

h

igh-order integration, spectral differentiation, numerical quadrature,
quadrilateral mesh

##  1 Introduction

Efficient numerical integration of surface integrals is an important
ingredient in applications ranging from _geometric processing_ Lachaud16 ,
_surface–interface and colloidal sciences_ zhou2005surface , optimization of
production processes Gregg_1967 ; riviere2009  to surface finite element
methods, solving partial differential equations on curved surfaces. Given a
conforming triangulation of a $2$-dimensional $C^{r+1}$-surface $S$, $r\geq
0$, a finite family of maps $\varrho_{i}$ and corresponding sets $V_{i}\subset
S$, $i=1,\dots,K$ such that

$\varrho_{i}:\triangle_{2}\to V_{i}\subseteq
S\,,\quad\bigcup_{i=1}^{K}\overline{V_{i}}=S\,,\quad\bigcap_{i\neq j}V_{i}\cap
V_{j}=\emptyset\,,$

and the restrictions of the $\varrho_{i}$ to the interior of the standard
simplex $\mathring{\triangle}_{2}$ are diffeomorphisms. For immersed manifolds
$S\subseteq\mathbb{R}^{3}$ we will write
$D\varrho_{i}(\mathrm{x}):\mathbb{R}^{2}\to\mathbb{R}^{3}$ for the Jacobian of
the parametrization $\varrho_{i}$ at $\mathrm{x}$. The surface integral of an
integrable function $f:S\to\mathbb{R}$ with $f\in C^{r}$, for each $V_{i}\in
S$ is

$\int_{S}f\,dS=\sum_{i=1}^{K}\int_{\triangle_{2}}f(\varrho_{i}(\mathrm{x}))g_{i}(\mathrm{x})\,d\mathrm{x}\,,$
(1)

here, the volume element is expressed as
$g_{i}(\mathrm{x})=\sqrt{\det((D\varrho_{i})^{T}(\mathrm{x})D\varrho_{i}(\mathrm{x}))}$.

Triangulations of an embedded manifold $S$ are frequently given as a set of
flat triangles in the embedding space $\mathbb{R}^{3}$, together with local
projections from these simplices onto $S$. More formally, let

$T_{i}\subseteq\mathbb{R}^{3}\,,\qquad i=1,\dots,K$ (2)

be a set of flat triangles. For each triangle $T_{i}$, we assume that there is
a well-defined $C^{r+1}$-embedding $\pi_{i}:T_{i}\to S$ and an invertible
affine transformation $\tau_{i}:\triangle_{2}\to T_{i}$, such that the maps
$\varrho_{i}=\pi_{i}\circ\tau_{i}:\triangle_{2}\to S$ form a conforming
triangulation. Commonly, _the closest-point projection_

$\pi^{*}:\mathcal{N}_{\delta}(S)\to
S\,,\quad\pi^{*}(\mathrm{x})=\operatorname*{arg\,min}_{\mathrm{y}\in
S}\mathrm{dist}(\mathrm{x},\mathrm{y})$

serves as a realization of the $\pi_{i}$. Recall from geometricpde  that given
an open neighborhood
$\mathcal{N}_{\delta}(S)=\\{\mathbf{x}\in\mathbb{R}^{m}:\mathrm{dist}(\mathbf{x},S)<\delta\\}$
of a $C^{r+1}$-surface $S$, $r\geq 2$, with $\delta$ bounded by the reciprocal
of the maximum of all principal curvatures on $S$, the closest-point
projection is well-defined on $\mathcal{N}_{\delta}(S)$ and of regularity
$\pi^{*}\in C^{r-1}(T,S)$.

Especially when it comes to handling complex geometries, in practice, simplex
meshes are typically much easier to obtain Persson . Here, we extend our
former approach zavalani2023highorder , involving the homeomorphic re-
parametrization

$\sigma_{*}:\square_{2}\to\triangle_{2}\,,\quad\triangle_{2}=\Big{\\{}\mathrm{x}\in\mathbb{R}^{2}:x_{1},x_{2}\geq
0\,,\sum_{i=1}^{2}|x_{i}|\leq 1\Big{\\}}\,,\square_{2}=[-1,1]^{2}$ (3)

referred to as _square-squeezing_ , for re-parametrizing initial triangular
meshes to quadrilateral ones. This enables pulling back interpolation and
integration tasks from the standard 2-simplex $\triangle_{2}$ to the square
$\square_{2}=[-1,1]^{2}$.

##  2 Square-squeezing re-parametrization

The main difficulty for providing a numerical approximation of (1) is
obtaining the unknown derivatives $D\varrho_{i}$ that appear in the volume
element. One classic approach, followed also by dziuk2013finite , is to
replace the Jacobians $D\varrho_{i}$ by the Jacobians of a polynomial
approximation, typically obtained by interpolation on a set of interpolation
points in $\triangle_{2}$. However, the question of how to distribute nodes in
simplices in order to enable stable high-order polynomial interpolation is
still not fully answered CHEN1995405 ; taylor2000generalized .

To bypass this issue, we instead propose to re-parametrize the curved
simplices $S_{i}$ over the square $\square_{2}=[-1,1]^{2}$. However, the
prominent _Duffy transformation_ Duffy82

$\sigma_{\text{Duffy}}:\square_{2}\to\triangle_{2}\,,\quad\sigma_{\text{Duffy}}(x,y)=\Big{(}\frac{1}{4}\left(1+x\right)\left(1-y\right),\frac{1+y}{2}\Big{)}\,,$
(4)

collapses the entire upper edge of the square to a single vertex, which
renders it to be no homeomorphism and, thus, to be infeasible for the re-
parametrization task.

To introduce the alternate _square-squeezing_ , we re-scale $\square_{2}$ to
$\square_{2}^{u}=[0,1]^{2}$, by setting $x_{1}=(x+1)/2$ and $x_{2}=(y+1)/2$,
and define

$\sigma:\square_{2}\to\triangle_{2}\,,\quad\sigma(x_{1},x_{2})=\Big{(}x_{1}-\frac{x_{1}x_{2}}{2},x_{2}-\frac{x_{1}x_{2}}{2}\Big{)}^{T}.$
(5)

The inverse map $\sigma^{-1}:\triangle_{2}\to\square_{2}$ is given by

$\sigma^{-1}(u,v)=\begin{pmatrix}1+\left(u-v\right)-\sqrt{\left(u-v\right)^{2}+4\left(1-u-v\right)}\\\
1-\left(u-v\right)-\sqrt{\left(u-v\right)^{2}+4\left(1-u-v\right)}\end{pmatrix}.$
(6)

Both $\sigma$ and $\sigma^{-1}$ are continuous, rendering square-squeezing to
be a homeomorphism on the closed set $\square_{2}$, see Fig. (1), and even a
smooth diffeomorphism in the interior, see zavalani2023highorder  for further
details.

######  Definition 1 (Quadrilateral re-parametrization)

Let $S$ be a $C^{r+1}$-surface, $r\geq 0$,
$\varrho_{i}=\pi_{i}\circ\tau_{i}:\triangle_{2}\to S$ a conforming
triangulation (2), and $\sigma:\square_{2}\to\triangle_{2}$ be a homeomorphism
whose restriction
$\sigma_{|\mathring{\square}_{2}}:\mathring{\square}_{2}\to\mathring{\triangle}_{2}$
to the interior is a $C^{r}$-diffeomorphism. We call

$\varphi_{i}:\square_{2}\to
S\,,\quad\varphi_{i}=\varrho_{i}\circ\sigma=\pi_{i}\circ\tau_{i}\circ\sigma\,,\quad
i=1,\dots,K\,,$ (7)

an $r$-regular quadrilateral re-parametrization, illustrated in Fig. 2.

(a) Standard square

(b) Duffy’s transformation

(c) Square-Squeezing

Figure 1: Bilinear square–simplex transformations: Deformations of equidistant
grids, under Duffy’s transformation (b) and square-squeezing (c)
$T_{i}\,$$\pi_{i}$$V_{i}\,$$\tau_{i}$$\triangle_{2}$$\sigma$$\square_{2}$
Figure 2: Construction of a surface parametrization over $\triangle_{2}$ by
closest-point projection from a piecewise affine approximate mesh, and re-
parametrization over the square $\Box_{2}$.

The quadrilateral re-parametrization enables interpolating the geometry
functions $\varphi_{i}=\varrho_{i}\circ\sigma:\square_{2}\to\mathbb{R}^{3}$ by
tensor-product polynomials, as we assert next.

##  3 Interpolation in the square

Our high-order format is based on interpolation in tensor-product
Chebyshev–Lobatto nodes in the square $\square_{2}$. We start by defining:

######  Definition 2 (Lagrange polynomials MIP ; minterpy )

Let $n\in\mathbb{N}$ and $G_{2,n}=\oplus_{i=1}^{2}\mathrm{Cheb}_{n}$ be the
tensorial Chebyshev–Lobatto grid, where
$\mathrm{Cheb}_{n}=\left\\{\cos\left(\frac{k\pi}{n}\right)\;:\;0\leq k\leq
n\right\\}$, indexed by a multi-index set
$A_{2,n}=\\{\alpha\in\mathbb{N}^{2}:\|\alpha\|_{\infty}\leq n\\}^{2}$. For
each $\alpha\in A_{2,n}$, the tensorial multivariate Lagrange polynomials are

$L_{\alpha}(x)=\prod_{i=1}^{2}l_{\alpha_{i},i}(x)\,,\quad
l_{j,i}(x)=\prod_{k=0,k\neq j}^{n}\frac{x_{i}-p_{k,i}}{p_{j,i}-p_{k,i}}\,.$
(8)

The Lagrange polynomials are a basis of the polynomial space
$\Pi_{2,n}=\mathrm{span}\\{\mathrm{x}^{\alpha}=x_{1}^{\alpha_{1}}x_{2}^{\alpha_{2}}:\alpha\in
A_{2,n}\\}$ induced by $A_{2,n}$. Since the $L_{\alpha}$ satisfy
$L_{\alpha}(p_{\beta})=\delta_{\alpha,\beta}$ for all $\alpha,\beta\in
A_{2,n}$, $p_{\alpha}\in G_{2,n}$, we deduce that given a function
$f:\square_{2}\to\mathbb{R}$, the interpolant $Q_{G_{2,n}}f\in\Pi_{2,n}$ of
$f$ in $G_{2,n}$ is

$Q_{G_{2,n}}f=\sum_{\alpha\in A_{2,n}}f(p_{\alpha})L_{\alpha}.$ (9)

######  Definition 3 ($k^{\text{th}}$-order quadrilateral re-parametrization)

Given an $r$-regular quadrilateral re-parametrization
$\varphi_{i}:\square_{2}\to
S\,,\quad\varphi_{i}=\varrho_{i}\circ\sigma=\pi_{i}\circ\tau_{i}\circ\sigma\,,\quad
i=1,\dots,K\,,$ mesh (7). We say that the mesh is of order $k$ if each element
has been provided as a set of nodes $\\{\varphi_{i}(p_{\alpha})\\}_{\alpha\in
A_{2,k}}$ sampled at $\\{p_{\alpha}\\}_{\alpha\in A_{2,k}}$.

This implies that on each element, we can approximate the $r$-regular
quadrilateral re-parametrization maps through interpolation using the nodes
$\\{\varphi_{i}(p_{\alpha})\\}_{\alpha\in A_{2,k}}$. This involves computing a
$k^{\text{th}}$-order vector-valued polynomial approximation:

$Q_{G_{2,k}}\varphi_{i}\left(\mathrm{x}\right)=\sum_{\alpha\in
A_{2,k}}\varphi_{i}(p_{\alpha})L_{\alpha},\quad i=1,\ldots,K$ (10)

The $r$-regular quadrilateral re-parametrization maps partial derivatives can
be computed using numerical spectral differentiation trefethen2000spectral .
These derivatives are then utilized to create an approximation of the metric
tensor, $g_{i}$, for each element.

To obtain partial derivatives, forming the Jacobian matrix
$DQ_{d,k}\varphi_{i}$ on the reference square $\square_{2}$, Kronecker
products are employed. Let $D\in\mathbb{C}^{(k+1)\times(k+1)}$ be the one-
dimensional spectral differentiation matrix trefethen2000spectral  associated
with Chebyshev–Lobatto nodes on the interval $[-1,1]$, and let
$I\in\mathbb{C}^{(k+1)\times(k+1)}$ be the identity matrix. The
differentiation matrices in the $x$ and $y$ directions on the reference square
$\square_{2}$ are defined as follows:

$D_{\mathrm{x}}Q_{G_{2,k}}\varphi_{i}=D\otimes I,\quad
D_{\mathrm{y}}Q_{G_{2,k}}\varphi_{i}=I\otimes D,$ (11)

where
$D_{\mathrm{x}}Q_{G_{2,k}}\varphi_{i},D_{\mathrm{y}}Q_{G_{2,k}}\varphi_{i}\in\mathbb{C}^{(k+1)^{2}\times(k+1)^{2}}$.

Representing the Jacobian as
$DQ_{d,k}\varphi_{i}=[D_{\mathrm{x}}Q_{d,k}\varphi_{i}\,|\,D_{\mathrm{y}}Q_{d,k}\varphi_{i}]$,
the ingredients above realise the HOSQ, computing the surface integral as:

$\displaystyle\int_{S}f\,dS$
$\displaystyle\approx\sum_{i=1}^{K}\int_{\square_{2}}f(\varphi(\mathrm{x}))\sqrt{\det((DQ_{d,k}\varphi_{i}(\mathrm{x}))^{T}DQ_{d,k}\varphi_{i}(\mathrm{x}))}\,d\mathrm{x}$
(12) $\displaystyle\approx\sum_{i=1}^{K}\sum_{\mathrm{p}\in
P}\omega_{\mathrm{p}}f(\varphi_{i}(\mathrm{p}))\sqrt{\det((DQ_{d,k}\varphi_{i}(p))^{T}DQ_{d,k}\varphi_{i}(p))}\,,$
(13)

where $\mathrm{p}\in P$, $\omega_{\mathrm{p}}\in\mathbb{R}^{+}$ can be the
points and weights of any quadrature of $\square_{2}$ e.g. the Gauss-Legendre
or Clenshaw-Curtis quadrature Xiang2012OnTC ; Trefethen2019 .

It is important to note that computing nodes and weights for a $n$-point
Gauss-Legendre rule requires $\mathcal{O}(n^{2})$ operations (error analysis
is given in zavalani2023highorder ), while the Clenshaw-Curtis method uses
$\mathcal{O}(n\log n)$ operations with the _Discrete Cosine Transform_ for
evaluation. We estimate the error of the latter:

###### Theorem 3.1

Consider a $C^{r}$ surface $S$, where $r\geq 1$, with $r$-regular
quadrilateral re-parametrization $\varphi_{i}:\square_{2}\to S$,

$\varphi_{i}=\varrho_{i}\circ\sigma=\pi_{i}\circ\tau_{i}\circ\sigma,\quad
i=1,\dots,K.$

Let $\mathrm{p}\in P$, and $\omega_{\mathrm{p}}$ be the points and weights of
the tensorial $(n+1)$-order Clenshaw-Curtis quadrature rule,
$f:S\to\mathbb{R}$ be a function with absolutely continuous derivatives up to
order $(r-1)$ and bounded variation $\|f^{(r)}\|_{T}=V<+\infty$, such that $f$
induces a negligible ”remainder of the remainder” elliott2011estimates . Then
the integration error can be estimated as

$\mathcal{E}_{f}=\sum_{i=1}^{K}\int_{\square_{2}}f(\varphi(\mathrm{x}))g_{i}\left(\mathrm{x}\right)d\mathrm{x}-\sum_{i=1}^{K}\sum_{\mathrm{p}\in
P}\omega_{\mathrm{p}}f(\varphi_{i}(\mathrm{p}))g_{i}\left(\mathrm{p}\right)\leq\frac{128\tilde{V}}{15\pi
r(n-r)^{r}},$ (14)

where
$g_{i}\left(\mathrm{x}\right)=\sqrt{\det((DQ_{d,k}\varphi_{i}(\mathrm{x}))^{T}DQ_{d,k}\varphi_{i}(\mathrm{x}))}$
and $\tilde{V}$ is defined as

$\tilde{V}=\max_{i=1,\ldots,K}\max_{\alpha,\beta\leq
r}\left\|\frac{\partial^{\alpha+\beta}}{\partial x^{\alpha}\partial
y^{\beta}}\Big{(}f(\varphi_{i}(\mathrm{x}))g_{i}(\mathrm{x})\Big{)}\right\|_{T}\;\text{with}\;\|\xi\|_{T}=\left\|\frac{\xi^{\prime}}{\sqrt{1-x^{2}}}\right\|_{1}.$
(15)

###### Proof

We denote $\psi(x,y):=f(\varphi(\mathrm{x}))g_{i}\left(\mathrm{x}\right)$ and

$R^{(2)}_{n}[\psi]=I^{(2)}[\psi]-Q^{(2)}_{n}[\psi]$ (16)

as the remainder of the exact integral and the $n+1$-order quadrature rule

$I^{(2)}[\psi]=\int_{\square_{2}}\psi(x,y)\,dx\,dy,\quad
Q^{(2)}_{n}[\psi]=\sum_{j=0}^{n}\sum_{i=0}^{n}\psi(x_{i},y_{j})\,w_{i}\,w_{j}\,.$
(17)

For a bivariate function $\psi(x,y)$,
$I_{x}[\psi]=\int_{-1}^{1}\psi(x,\cdot)\,dx$ denotes the integration with
respect to the $x$ variable only, yielding a function of $y$. The subscript
notation extends to $R_{n,x}[\psi]$ and $Q_{n,y}[\psi]$, and when replacing
the roles of $x$ and $y$. Fubini’s theorem Brezis2010FunctionalAS  implies:

$I^{(2)}[\psi]=\int_{\square_{2}}\psi(x,y)\,dx\,dy=\int_{-1}^{1}I_{x}[\psi]\,dy=I_{y}[I_{x}[\psi]]\,.$
(18)

Upon substitution into Eq. (18), we obtain

$\displaystyle I^{(2)}[\psi]$
$\displaystyle=I_{y}[R_{n,x}[\psi]+Q_{n,x}[\psi]]$
$\displaystyle=R_{n,y}[R_{n,x}[\psi]+Q_{n,x}[\psi]]+Q_{n,y}[R_{n,x}[\psi]+Q_{n,x}[\psi]]$
$\displaystyle=R_{n,y}[R_{n,x}[\psi]]+Q_{n,x}[R_{n,y}[\psi]]+Q_{n,y}[R_{n,x}[\psi]]+Q_{n,y}[Q_{n,x}[\psi]]$
(19)

Following elliott2011estimates  we assume the “remainder of a remainder” –
first term in Eq. (3) to contribute negligibly to the error, enabling us to
establish a sufficiently tight upper bound:

For large $n$, the quadrature rule approaches the value of the integral, i.e.,
$Q_{n,\beta}\approx I_{\beta}$ for $\beta=x,y$, we are left with:

$I^{(2)}[\psi]\approx
I_{x}[R_{n,y}[\psi]]+I_{y}[R_{n,x}[\psi]]+Q^{(2)}_{n}[\psi]\,.$ (20)

Hence,

$R^{(2)}_{n}[\psi]\leq I_{x}[R_{n,y}[\psi]]+I_{y}[R_{n,x}[\psi]],$ (21)

As noted in Trefethen2019 , considering a function $\theta(x)$ defined on the
interval $[-1,1]$, when computing $Q_{n}[\theta]$ using Clenshaw-Curtis
quadrature for $\theta\in C^{r}$ and $\|\theta^{(r)}\|_{T}<V$ for a real
finite value $V$, then for sufficiently large $n$, the subsequent inequality
holds:

$R_{n}[\theta]\leq\frac{32V}{15\pi r(n-r)^{r}}$ (22)

Consequently, applying Eq. (22) to Eq. (21) yields

$R^{(2)}_{n}[\psi]\leq\frac{32}{15\pi
r(n-r)^{r}}[I_{x}[V_{y}(x)]+I_{y}[V_{x}(y)]]\,,$ (23)

where $V_{y}(x)=\max_{y}\|\psi^{(r)}(x,y)\|_{T}$ and
$V_{x}(y)=\max_{x}\|\psi^{(r)}(x,y)\|_{T}$ for fixed $x$ and $y$,
respectively.

By utilizing the definition of the variation $\tilde{V}$ in Eq. (15), we are
left with Eq. (14). Applying Clenshaw-Curtis quadrature leads to an error
approximation for a 2D surface integral with an order of
$\mathcal{O}\big{(}n^{-r}\big{)}$, where $n$ denotes the quadrature’s order.

##  4 Numerical experiments

We demonstrate the effectiveness of HOSQ-CC (utilizing Clenshaw-Curtis
quadrature) for surface integration through two set of numerical examples. We
provide the initial flat surface triangulations by using Persson and Strang’s
algorithm Persson  and subsequently enhance them to curved triangulations
using Euclidean closest-point projections $\pi_{i}:T_{i}\to S$. Our
implementation of HOSQ is part of a Python package surfpy 111available at
https://github.com/casus/surfpy . The examples and results of this manuscript
using Dune-CurvedGrid are summarized and made available in the
repository.222https://github.com/casus/dune-surface_int

We compare HOSQ-CC with Dune-CurvedGrid integration algorithm (DCG) from the
surface-parametrization module dune-curvedgrid CurvedGrid , part of the Dune
finite element framework333www.dune-project.org. As outlined in
zavalani2024note , DCG performs total degree interpolation in uniform
(equidistant midpoint-rule) triangle-nodes for interpolating the closest-point
projection on each triangle, which makes it sensitive to Runge’s phenomenon
(overfitting).

(a) Unit sphere

(b) Torus with radii $R=2$ and $r=1$

Figure 3: Relative errors and runtimes between DCG and HOSQ-CC computing the
surface area of the unit sphere (3(a)) and torus (3(b)) respectively.

We start, by addressing integration tasks involving only constant integrands
$f:S\to\mathbb{R}$.

######  Experiment 1 (Area computation)

We integrate $f=1$ over the unit sphere $\mathrm{S}^{2}$ and the torus
$T^{2}_{r,R}$ with inner radius $r=1$ and outer radius $R=2$. The surface
areas are given due to $4\pi r^{2}$ and $4\pi^{2}Rr$, respectively, enabling
measuring the relative errors.

We choose initial triangulations of size $N_{\Delta}=124$ for the sphere and
of size $N_{\Delta}=256$ for the torus and apply Clenshaw-Curtis
quadratureTrefethen2019 , with a degree matching the geometry approximation.
We use symmetric Gauss quadrature on a simplex with a matching degree for the
geometry approximation in the case of DCG.

HOSQ-CC exhibits stable convergence to machine precision with exponential
rates, $0.05\times 10.8^{-n}$ and $0.025\times 7.5^{-n}$ fitted for the sphere
and torus, respectively, in accordance with the predictions from Theorem 3.1.
In contrast, DCG becomes unstable for orders larger than $\deg=10$. We
interpret the instability as the appearance of Runge’s phenomenon caused by
the choice of midpoint-triangle-refinements yielding equidistant interpolation
nodes for DCG.

In terms of execution time, the Python implementation of HOSQ-CC, empowered by
spectral differentiation, outperforms the $\text{C}++$ implementation of DCG.
Importantly, experimental findings reveal that DCG encounters constraints for
$\deg\geq 15$, attributable to the formation of a singular matrix during the
basis computation.

######  Experiment 2 (Gauss-Bonet validation)

We consider the Gauss curvature as a non-trivial integrand. Due to the
Gauss–Bonnet theorem spivak1999 , integrating the Gauss curvature over a
closed surface yields

$\int_{S}K_{\mathrm{Gauss}}\,dS=2\pi\chi\left(S\right)\,,$ (24)

where $\chi\left(S\right)$ denotes the Euler characteristic of the surface.
Here, we consider:

1) | Dziuk’s surface | $(x-z^{2})^{2}+y^{2}+z^{2}-1=0$
---|---|---
2) | Double torus |  $\big{(}x^{2}+y^{2})^{2}-x^{2}+y^{2}\big{)}^{2}+z^{2}-a^{2}=0$ , $a\in\mathbb{R}\setminus\\{0\\}$

The Gauss curvature is computed symbolically from the implicit surface
descriptions using Mathematica 11.3, enabling measuring errors of DCG and
HOSQ-CC when integrating the Gauss curvature. We maintain the experimental
design outlined in Experiment 1 and display error plots based on the
polynomial degree in Fig.4 and Fig.5.

(a) Dziuk’s surface with $8088$ triangles

(b) Integration errors and and runtimes

Figure 4: Gauss-Bonnet validation for Dziuk’s surface.

(a) Double torus, with $8360$ triangles and $a=0.2$

(b) Integration errors and and runtimes

Figure 5: Gauss-Bonnet validation for a double torus.

The HOSQ-CC consistently exhibits exponential convergence to the accurate
value $2\pi\chi(S)$ with exponential rates, $4.5\cdot 10^{-5}\times 3.7^{-n}$
and $5\cdot 10^{-4}\times 6.6^{-n}$ fitted for the for the Dziuk’s surface and
the double torus, respectively. The best fit of an algebraic rate, $3.95\times
n^{-12.16}$ for the Dziuk’s surface and $1.7\times n^{-12.3}$ for the double
torus, does not assert rapid convergence. As the order increases, the error in
HOSQ-CC tends to stabilize near the level of machine precision. In contrast,
DCG fails to achieve machine-precision approximations across all cases and
becomes unstable for degrees $k\geq 10$. In terms of execution time, once
more, HOSQ-CC, empowered by spectral differentiation, outperforms DCG.

##  5 Conclusion

The present HOSQ integration approach, utilizing the innovative _square-
squeezing_ transformation excels in speed, accuracy, and robustness for
integration task on complex geometries, suggesting its application potential
for triangular spectral element methods (TSEM) karniadakis2005spectral ;
heinrichs2001spectral  and fast spectral PDE solvers on surfaces
fortunato2022high .

## Acknowledgement

This research was partially funded by the Center of Advanced Systems
Understanding (CASUS), which is financed by Germany’s Federal Ministry of
Education and Research (BMBF) and by the Saxon Ministry for Science, Culture,
and Tourism (SMWK) with tax funds on the basis of the budget approved by the
Saxon State Parliament.

## References

  * [1] A. Bonito and R. H. Nochetto.  Geometric Partial Differential Equations — Part I.  Elsevier, 2020.
  * [2] H. Brezis.  Functional analysis, sobolev spaces and partial differential equations.  2010\.
  * [3] Q. Chen and I. Babuška.  Approximate optimal points for polynomial interpolation of real functions in an interval and in a triangle.  Computer Methods in Applied Mechanics and Engineering, 128(3):405–417, 1995.
  * [4] M. G. Duffy.  Quadrature over a pyramid or cube of integrands with a singularity at a vertex.  SIAM Journal on Numerical Analysis, 19:1260–1262, 1982.
  * [5] G. Dziuk and C. M. Elliott.  Finite element methods for surface PDEs.  Acta Numerica, 22:289–396, 2013.
  * [6] D. Elliott, P. R. Johnston, and B. M. Johnston.  Estimates of the error in gauss–legendre quadrature for double integrals.  Journal of Computational and Applied Mathematics, 236(6):1552–1561, 2011.
  * [7] D. Fortunato.  A high-order fast direct solver for surface PDEs.  arXiv preprint arXiv:2210.00022, 2022.
  * [8] S. J. Gregg, K. S. W. Sing, and H. W. Salzberg.  Adsorption surface area and porosity.  Journal of The Electrochemical Society, 114(11):279, nov 1967.
  * [9] M. Hecht, K. Gonciarz, J. Michelfeit, V. Sivkin, and I. F. Sbalzarini.  Multivariate interpolation in unisolvent nodes–lifting the curse of dimensionality.  arXiv preprint arXiv:2010.10824, 2020.
  * [10] W. Heinrichs and B. I. Loch.  Spectral schemes on triangular elements.  Journal of Computational Physics, 173(1):279–301, 2001.
  * [11] U. Hernandez Acosta, S. Krishnan Thekke Veettil, D. Wicaksono, and M. Hecht.  minterpy – Multivariate interpolation in Python.  https://github.com/casus/minterpy, 2021.
  * [12] G. E. Karniadakis, G. Karniadakis, and S. Sherwin.  Spectral/$hp$ element methods for computational fluid dynamics.  Oxford University Press, 2005.
  * [13] J.-O. Lachaud.  Convergent geometric estimators with digital volume and surface integrals.  In Discrete Geometry for Computer Imagery, pages 3–17. Springer, 2016.
  * [14] P.-O. Persson and G. Strang.  A simple mesh generator in MATLAB.  SIAM Review, 46(2):329–345, 2004.
  * [15] S. Praetorius and F. Stenger.  Dune-CurvedGrid – a Dune module for surface parametrization.  Archive of Numerical Software, page Vol. 1 No. 1 (2022), 2022.
  * [16] J. C. Riviere and S. Myhra.  Handbook of surface and interface analysis: methods for problem-solving.  CRC press, 2009.
  * [17] M. Spivak.  A Comprehensive Introduction to Differential Geometry, volume 1.  Publish or Perish Incorporated, 1999.
  * [18] M. A. Taylor and B. Wingate.  A generalized diagonal mass matrix spectral element method for non-quadrilateral elements.  Applied Numerical Mathematics, 33(1-4):259–265, 2000.
  * [19] L. N. Trefethen.  Spectral Methods in MATLAB.  Society for Industrial and Applied Mathematics, 2000.
  * [20] L. N. Trefethen.  Approximation theory and approximation practice, volume 164.  SIAM, 2019.
  * [21] S. Xiang and F. A. Bornemann.  On the convergence rates of gauss and clenshaw-curtis quadrature for functions of limited regularity.  SIAM J. Numer. Anal., 50:2581–2587, 2012.
  * [22] G. Zavalani, O. Sander, and M. Hecht.  High-order integration on regular triangulated manifolds reaches super-algebraic approximation rates through cubical re-parameterizations.  arXiv preprint arXiv:2311.13909, 2023.
  * [23] G. Zavalani, E. Shehu, and M. Hecht.  A note on the rate of convergence of integration schemes for closed surfaces.  Computational and Applied Mathematics, 43(2):1–17, 2024.
  * [24] Q. Zhou and P. Somasundaran.  Surface and Interfacial Tension: Measurement, Theory, and Applications. Surfactant Science Series, volume 119.  ACS Publications, 2005.