# On the trajectory of the nonlinear pendulum:
Exact analytic solutions via power series.

W. Cade Reinberger  School of Mathematical Sciences, Rochester Institute of
Technology, Rochester, NY 14623, USA  Morgan S. Holland  School of
Mathematical Sciences, Rochester Institute of Technology, Rochester, NY 14623,
USA  Nathaniel S. Barlow  School of Mathematical Sciences, Rochester Institute
of Technology, Rochester, NY 14623, USA  Steven J. Weinstein  Department of
Chemical Engineering, Rochester Institute of Technology, Rochester, NY 14623,
USA  School of Mathematical Sciences, Rochester Institute of Technology,
Rochester, NY 14623, USA

###### Abstract

We provide an exact infinite power series solution that describes the
trajectory of a nonlinear simple pendulum undergoing librating and rotating
motion for all time. Although the series coefficients were previously given in
[V. Fairén, V. López, and L. Conde, Am. J. Phys 56 (1), (1988), pp. 57?61],
the series itself – as well as the optimal location about which an expansion
should be chosen to assure series convergence and maximize the domain of
convergence – was not examined, and is provided here. By virtue of its
representation as an elliptic function, the pendulum function has
singularities that lie off of the real axis in the complex time plane. This,
in turn, imposes a radius of convergence on the physical problem in real time.
By choosing the expansion point at the top of the trajectory, the power series
converges all the way to the bottom of the trajectory without being affected
by these singularities. In constructing the series solution, we re-derive the
coefficients using an alternative approach that generalizes to other nonlinear
problems of mathematical physics. Additionally, we provide an exact
resummation of the pendulum series – Motivated by the asymptotic approximant
method given in [Barlow et al., Q. J. Mech. Appl. Math., 70 (1) (2017), pp.
21-48] – that accelerates the series’ convergence uniformly from the top to
the bottom of the trajectory. We also provide an accelerated exact resummation
of the infinite series representation for the elliptic integral used in
calculating the period of a pendulum’s trajectory. This allows one to preserve
analyticity in the use of the period to extend the pendulum series for all
time via symmetry.

##  I Introduction

The pendulum equation is often presented in undergraduate science and math
courses as a physically-relevant nonlinear ordinary differential equation
(ODE) whose general solution cannot be written in terms of elementary
functions. The ODE is used to motivate phase-plane and stability analyses
Simmons ; Strogatz  and, when linearized for small deflections, admits a
classical harmonic oscillator solution. It has been known for over a century
that the analytic solution to the pendulum equation, expressed as angle as a
function of time, may be written in terms of Jacobi elliptic functions
Helmholtz ; Appell ; Whittaker . The inverse solution (time as a function of
angle) is expressed in terms of elliptic integrals and has been well
documented for even longer; its application towards the evaluation of the
period of a pendulum’s trajectory was established alongside those of elliptic
integrals in the 18th century Stillwell ; Gray . The singularities of the
Jacobi elliptic functions (and thus the pendulum equation itself) were studied
during the early development of the field of complex analysis Gray  and are
now well-documented Byrd . The efficient representation of Jacobi elliptic
functions via power series is limited by the radius of convergence imparted by
such singularities. While the use of power series solutions for linear ODEs is
commonplace, the analogous treatment for their nonlinear counterparts is only
found in scattered pockets of the literature. In general, power series
solutions to nonlinear ODEs have been thought by many (in the words of Agnew
p. 263) to “ play minor or nonexistent roles in the lives of most people”. It
is perhaps not surprising that a power series solution – for the pendulum
angle as a function of time – is virtually absent from the literature,
especially since numerical algorithms for elliptic functions are readily
available in many scientific computing packages. An early attempt at finding
the power series solution to the pendulum equation is given in Agnew  (p.
264), where substitution of the initial conditions directly into the ODE and
then successive (and objectively arduous) applications the chain rule leads
one after four terms to (again in the author’s words) “feed the wastebasket
and go to work on the next chapter”. A few decades later, a methodology was
developed Fairen  to obtain the recursion for all coefficients of the pendulum
series through a clever decomposition of a system of recursions. That said,
the series of Fairen  was implied as being approximate (by the paper’s title),
and a study of the radius of convergence of the series was not provided. In
fact, the solution itself is incomplete for that reason, as the initial
conditions are left arbitrary, and convergence is not guaranteed. The aim of
this work is to show that the power series representation of the pendulum
function given by Fairen  is, in fact, exact and analytic over the entire
physical domain if one judiciously chooses the initial conditions (or
equivalently, expansion point) along a given trajectory to assure convergence.
Once a convergent series is constructed through this choice of initialization,
the series may be re-summed exactly to accelerate convergence, and thus can be
used as an alternative to elliptic functions. It is relevant to note here
that, although expansions for Jacobi elliptic functions are well-documented
Byrd , the number of arithmetic steps needed needed to convert such expansions
to a single power series for the pendulum solution is larger than those to
directly find the power series solution of the original nonlinear ODE; we
demonstrate this through the analysis in the current work. More generally,
insights provided in this paper regarding power series solutions to nonlinear
ODEs are relevant to problems that cannot be expressed by special functions
such as Jacobi elliptic functions.

In general, power series solutions to ODEs are limited by their radius of
convergence as mentioned above. In linear ODEs, the radius of convergence may
be anticipated based on the distance between the expansion point and
singularities in ODE coefficients. In nonlinear ODE’s, singularities may be
spontaneous and cannot be anticipated by inspection. As such, the radius of
convergence imparted to power series solutions of nonlinear ODEs can restrict
an analytic representation over the whole domain of interest - or, can lead to
slow convergence properties. Thus, power series solutions are often deficient,
and with the prevalence of numerical techniques nowadays, often ignored as a
viable solution technique. However, techniques have been developed to
accelerate convergence and even sum divergent series beyond the convergence
radius (i.e., analytic continuation). Approximant techniques, such as the
well-known Padé, have been utilized to this end Boyd1997 . A relatively new
approach, asymptotic approximants Barlow:2017 , has been successfully applied
to power series expansions to create highly accurate approximate solutions of
nonlinear ODEs arising in many areas of mathematical physics FS2020 ; SIR2020
; SEIR2020 ; Harkin ; Rame . In particular, if one knows the asymptotic
behavior in the vicinity of the boundaries of the physical domain, the method
of asymptotic approximants may be used to constrain analytic continuation by
enforcing these behaviors at both ends, thus providing efficient and accurate
analytic infinite series solutions which are in closed form when truncated. In
cases where a power series solution converges over the whole domain of
interest, but convergence is slow, the same approach may be used to re-sum the
series exactly and enhance convergence Harkin ; Rame . In this paper, we
establish that, with a judicious choice of initial conditions, a convergent
power series solution to the nonlinear pendulum ODE may be obtained. The
convergence of the series solution is then accelerated by exactly re-summing
it in a form motivated by the asymptotic approximant technique. We should note
that the word ‘approximant’ is not to be confused with ‘approximate’ and is
only used here to tie to previous literature. The infinite series solution we
present here, both in its original and re-summed form, is an exact solution.

While one goal of this work is to assemble the above-described solution
elements from across the centuries in order to give a full treatment of the
pendulum series for direct use, another aim here is to add this example to a
growing list of problems Barlow:2017 ; FS2020 ; SIR2020 ; SEIR2020 ; Harkin ;
Rame  in which power series solutions are demonstrated to be a viable solution
technique for nonlinear ODEs, even in the more typical cases where singularity
locations are unknown. The paper is organized as follows. In section II we
detail the key structural elements of the pendulum equation, as well as its
solution in terms of Jacobi elliptic functions. In section III, we demonstrate
an equivalent but alternative approach to Fairen  for deriving the power
series solution to the differential equation. Using the known placement of the
singularities in the complex time plane, we illustrate that the point of
expansion – or equivalently, initial conditions – can be chosen such that the
series solution encompasses the entire physical domain, using reflections
about a judiciously chosen effective period. In section IV, the method of
asymptotic approximants is used to motivate an exact resummation of the
pendulum series, and details are provided on its construction and use.
Appendix A provides an efficient way of implementing truncations of this
resummation. In section V, we develop an exact series resummation for the
complete elliptic integral of the first kind that is used in period
calculations. This period resummation, in turn, is utilized to periodically
extend the solution forms of sections III and IV. Concluding remarks are given
in section VI on the pendulum series solution provided here. This solution
adds to the growing body of literature in which nonlinear ODEs are solved
using series solutions, and so some closing remarks on this connection are
provided as well.

##  II Background

$m$$\ell$$\theta$ Figure 1: Setup of the simple pendulum problem. The mass,
$m$, at the end of the string (of length $\ell$) swings freely from some
initial angle and angular velocity. In this reference frame, the acceleration
of gravity, $g$, acts vertically downward.

A schematic of the problem to be examined is provided in figure 1, which shows
a pendulum bob of point mass $m$ anchored by a mass-less string of length
$\ell$. During motion, the bob is displaced from the vertical by an angle
$\theta$, acted upon by a downward gravitational force $mg$ ($g$ being the
acceleration of gravity), and moves with an angular acceleration of
$d^{2}\theta/dt^{2}$ where $t$ is time. Subject to the above assumptions, a
moment balance leads to the “simple pendulum equation” Agnew :

$\frac{d^{2}\theta}{dt^{2}}+\frac{g}{\ell}\sin(\theta)=0.$

Defining the dimensionless time $\bar{t}=t\sqrt{\frac{g}{\ell}}$ as done in
Strogatz , the initial value problem for the simple pendulum then becomes

$\frac{d^{2}\theta}{d\bar{t}^{2}}+\sin(\theta)=0$ (1a)
$\theta(0)=\theta_{0},~{}~{}\frac{d\theta}{d\bar{t}}(0)=\bar{\omega}_{0},$
(1b)

where $\theta_{0}$ and $\bar{\omega}_{0}$ are the initial angular displacement
and dimensionless angular velocity111The dimensionless angular velocity is
defined in terms of its dimensional equivalent, $\omega_{0}$, as
$\bar{\omega}_{0}=\omega_{0}\sqrt{\frac{\ell}{g}}$., respectively. An energy
constant for the simple pendulum is obtained by multiplying both sides of (1a)
by $d\theta/d\bar{t}$, integrating both sides with respect to $\bar{t}$, and
applying (1b) to obtain

$\frac{1}{2}\left(\frac{d\theta}{d\bar{t}}\right)^{2}+1-\cos(\theta)=\frac{1}{2}\bar{\omega}_{0}^{2}+1-\cos(\theta_{0})\equiv\bar{E}\geq
0.$ (2)

In (2), $\bar{E}$ is the dimensionless total energy222The dimensionless total
energy is defined in terms of its dimensional equivalent, $E$, as
$\bar{E}=\frac{E}{mg\ell}$., which is zero when the pendulum is resting at
zero angle.

Equation (1b) – and the types of solutions it admits – have been studied
extensively in prior literature Whittaker ; Lima ; Ochs ; chapter . All
possible moving (ignoring fixed points) solutions to (1b) are shown in figure
2. Initial conditions used for figure 2 are purposely chosen to be
$\theta_{0}$=0 while letting $\bar{\omega}_{0}$ vary to determine $\bar{E}$
according to (2) (i.e., pushing the pendulum from the bottom) since these
constraints are consistent with all possible pendulum motion. For other
combinations of initial conditions, one only needs to compute $\bar{E}$,
locate the initial conditions in figure 2 (shifting the $\bar{t}$ origin to
that point), and then follow the fixed $\bar{E}$ trajectory from that point
forward - this is consistent with the properties of autonomous ODEs. The
$\text{sgn}(\omega_{0})$ in figure 2 is used as a pre-factor on $\theta$ to
indicate that $\omega_{0}>0$ leads to the surface of solutions as shown while
$\omega_{0}<0$ leads to a reflection of the surface about the $t-\bar{E}$
plane.

Figure 2:  The solution to (1b) describing the trajectory of a pendulum over
time for several energy levels $\bar{E}$, encompassing the 3 possibilities:
libration ($\bar{E}<2$), the separatrix orbit ($\bar{E}$=2), and rotation
($\bar{E}>2$). The location of $\bar{t}=0$ is chosen to be at the bottom of
the pendulum’s trajectory, such that all energy states can be shown from the
same starting point. The $\text{sgn}(\omega_{0})$ is used as a pre-factor on
$\theta$ to indicate that $\omega_{0}>0$ leads to the surface of solutions as
shown while $\omega_{0}<0$ leads to a reflection of the surface about the
$t-\bar{E}$ plane. The solutions shown may be considered exact to within
machine precision and are generated using analytic formulae (23), (24d), and
(28) developed in sections III.3, IV, and V.

The types of pendulum solutions admitted by (1b) and shown in figure 2 are
classified by their value of $\bar{E}\in[0,\infty)$, approaching the standard
closed-form linear approximations as $\bar{E}\to 0$ and $\bar{E}\to\infty$. As
$\bar{E}\to 0$, simple harmonic motion arises in the form
$\theta\approx\theta_{0}\cos(\bar{t})+\bar{\omega}_{0}\sin(\bar{t})$, this
limit being accessed for small $\theta$ in (1b) such that
$\sin(\theta)\approx\theta$. As $\bar{E}\to\infty$, uniform circular motion
arises in the form $\theta\approx\bar{\omega}_{0}t$, this limit being accessed
by making $d\theta/d\bar{t}$ large in (2) such that the potential energy term
$1-\cos(\theta_{0})$ (always of $O$(1)) can be neglected, allowing direct
integration of (2).

In reality, the system (1b) admits 3 types of pendulum motion, as shown in
figure 2: the librating (i.e., “swinging” back and forth) trajectory for
$\bar{E}<2$, the separatrix orbit for $\bar{E}=2$ (dashed line in figure) in
which the pendulum rises towards the top, never reaching it yet never falling
back down, and the rotating or “overturning” case for $\bar{E}>2$.

For definiteness, we define the pendulum’s period $T$ as follows. For
libration ($\bar{E}<2$), we use the usual definition of period, such that
$\theta$ satisfies $\theta(\bar{t})=\theta(\bar{t}\pm T)$ while for rotation
($\bar{E}>2$) we use the unwound definition $\theta(\bar{t})=\theta(\bar{t}\pm
T)\mp\text{sgn}(\bar{\omega}_{0})2\pi$ where $\text{sgn}(\bar{\omega}_{0})$ is
1 for counter-clockwise rotation and $-1$ for clockwise rotation. For the
separatrix orbit ($\bar{E}=2$), $T$ is infinite and the solution is not
periodic. As mentioned in section I, the evaluation of $T$ itself dates back
to the 18th century Stillwell ; Gray ; it is obtained by applying separation
of variables to (2), solving for $\bar{t}$ in terms of an integral in
$\theta$, and then either: (for $\bar{E}<2$) letting $\bar{\omega}_{0}$=0,
integrating between the quarter-period limits ($\theta=0$ to
$\theta=\theta_{0}$), and multiplying by 4; or (for $\bar{E}>2$) letting
$\theta_{0}=0$, integrating between the half-period limits ($\theta=0$ to
$\theta=\pi$) and multiplying by 2. Although these integrals cannot be
evaluated in closed-form, they may be manipulated into elliptic integral form
by applying the half-angle identity, employing a variable substitution333For
$\bar{E}<2$, let $\sin\frac{\theta_{0}}{2}\sin u=\sin\frac{\theta}{2}$. For
$\bar{E}>2$, let $u=\frac{\theta}{2}$, and a rewriting the resulting
expression in terms of $\bar{E}$ to obtain

$T=\begin{cases}4~{}K\\!\left(\sqrt{\frac{\bar{E}}{2}}\right)&\text{if
$\bar{E}<2$ }\\\
2\sqrt{\frac{2}{\bar{E}}}~{}K\\!\left(\sqrt{\frac{2}{\bar{E}}}\right)&\text{if
$\bar{E}>2$}.\end{cases}$ (3)

where

$K(k)=\int_{0}^{\frac{\pi}{2}}\frac{d\phi}{\sqrt{1-k^{2}\sin^{2}(\phi)}}.$ (4)

In (4), $K(k)$ is the complete elliptic integral of the first kind Byrd .

While figure 2 shows trajectory solutions of (1b) over several periods, all
information above and below the separatrix for all $\bar{t}$ is, in fact,
obtainable from reflections and translations of the solution in a finite
interval $t\in[0,T^{*}]$, afforded by the symmetry of the problem. We define
the effective period $T^{*}$ to be the minimum portion of the period necessary
to map out the entire solution. In the case of libration, the pendulum sweeps
out the same (albeit negative) angles in time as it moves from $\theta=0$ to
$-\theta_{0}$ as it does going from $\theta=0$ to $\theta=\theta_{0}$. Thus
$T^{*}$ for a librating pendulum is the quarter-period. In the case of
rotation, instead of the back-and-forth motion described above, the pendulum
swings “over the top” and ends up back where it started from. Here, the
pendulum sweeps out the same angles (albeit with a $\pi$ offset) in time as it
moves from $\theta=\theta_{0}$ to
$\theta_{0}+~{}\text{sgn}(\bar{\omega}_{0})\pi$ as it does going from
$\theta=\theta_{0}+\text{sgn}(\bar{\omega}_{0})\pi$ to
$\theta_{0}+2~{}\text{sgn}(\bar{\omega}_{0})\pi$. Thus $T^{*}$ for a rotating
pendulum is the half-period. To summarize:

$T^{*}=\begin{cases}\frac{T}{4}&\text{if $\bar{E}<2$}\\\ \frac{T}{2}&\text{if
$\bar{E}>2$},\end{cases}$ (5)

where the separatrix orbit ($\bar{E}=2$) is omitted from the above definition
since it is not periodic.

Interestingly, initial conditions corresponding to $\bar{E}=2$ lead to
simplifications such that the elliptic integral for $\bar{t}(\theta)$
(obtained by separating (2)) can be evaluated in closed form and $\theta(t)$
can be solved for explicitly Lima 444Note that, in arriving at (6a) through
the route of elliptic integrals described here, $\theta_{0}$ is taken to be 0
in Lima , whereas we leave $\theta_{0}$ to be arbitrary. An alternative
derivation (starting from elliptic functions) is given in Whittaker  where
$\theta_{0}$ is also taken to be 0. as

$\theta=-\pi+4~{}\text{arctan}\left[e^{\bar{t}}\tan\left(\frac{\theta_{0}+\pi}{4}\right)\right],\text{
for }~{}~{}\bar{E}=2.$ (6a) For the libration and rotation cases, no known
solution exists for (1b) in terms of elementary functions; however, as
mentioned in section I, the solution can be written in terms of Jacobi
elliptic functions. These solutions are obtained by making the substitution
$y=\sin(\theta/2)$ into the energy equation (2), enforcing initial conditions
for either $\bar{E}<2$ or $\bar{E}>2$, and recognizing the resulting ODE as
that governing an elliptic function chapter . For libration and rotation, the
solutions are respectively,
$\theta=2\arcsin\left[\sqrt{\frac{\bar{E}}{2}}\operatorname{cd}\left(\bar{t},\sqrt{\frac{\bar{E}}{2}}\right)\right],\text{
for }~{}~{}\bar{E}<2,~{}\bar{\omega}_{0}=0,$ (6b) and
$\theta=2\arcsin\left[\operatorname{sn}\left(\sqrt{\frac{\bar{E}}{2}}\bar{t},\sqrt{\frac{2}{\bar{E}}}\right)\right],\text{
for }~{}~{}\bar{E}>2,~{}\theta_{0}=0$ (6c)

for a given $\bar{E}$. In (6a), (6b), and (6c), note that, once the energy
$\bar{E}$ is set, only one initial condition (either $\theta_{0}$ or
$\bar{\omega}_{0}$) is needed to specify the solution. The remaining initial
condition may then be extracted from (2). As stated previously, all initial
conditions can be accessed by shifting solutions of the same energy $\bar{E}$.
As there can be notational differences in the literature, for definiteness we
utilize the convention Byrd ; chapter  in (6b) and (6c) that the Jacobi
elliptic functions $\operatorname{cd}(u,k)$ and $\operatorname{sn}(u,k)$
satisfy the nonlinear initial value problems

$\left(\frac{d\operatorname{cd}}{du}\right)^{2}=\left(1-\operatorname{cd}^{2}\right)\left(1-k^{2}\operatorname{cd}^{2}\right),~{}~{}\operatorname{cd}(0,k)=1$

and

$\left(\frac{d\operatorname{sn}}{du}\right)^{2}=\left(1-\operatorname{sn}^{2}\right)\left(1-k^{2}\operatorname{sn}^{2}\right),~{}~{}\operatorname{sn}(0,k)=0.$

The solution of these equations is not known analytically in closed form;
however, efficient numerical algorithms exist to evaluate the Jacobi elliptic
functions in many scientific computing packages. As an alternative to the
numerical solution of (1b) or the solutions forms (6b) and (6c), in section
III of this paper we generate an exact solution to (1b) via a convergent power
series, whose convergence is subsequently accelerated in section IV. Although
the solution of the pendulum equation in this paper provide an alternative to
the numerical evaluation of the elliptic functions, the expressions (6b) and
(6c) are important when making exact statements about singularity locations in
section III. Note that there are several other equivalent representations of
the solution to (1b) involving other Jacobi elliptic functions. In fact, one
may also express the solution as a single elliptic function that governs all
pendulum movement Ochs .

##  III The Pendulum Series and Complex-Time Singularities

The power series solution of (1b) follows a standard approach Agnew , and is
developed in section III.1; this leads to a recurrence formula for its
coefficients, equivalent to that obtained via a different technique in Fairen
. However, the previous solution Fairen  leaves implementation to the reader
with arbitrary initial conditions and does not examine convergence properties
of the series, leaving the solution itself as incomplete. The radius of
convergence of the series based on closest singularities is determined in
section III.2; through the placement of judicious expansion points (or,
equivalently, initial conditions) this enables us to create a convergent
series in the time range $t\in[0,T^{*}]$ (see (5) and surrounding discussion).
This solution is subsequently utilized to construct the solution to (1b) for
all time in section III.3.

###  III.1 Power Series Solution to (1b) for $\theta(t)$

The power series solution may be generated by assuming

$\theta=\sum_{n=0}^{\infty}a_{n}\bar{t}^{n},$

substituting this into (1b), and equating like-powers of $\bar{t}$ on each
side. The result is a recurrence expressed as

$a_{n+2}=\frac{-s_{n}}{(n+1)(n+2)},~{}a_{0}=\theta_{0},~{}a_{1}=\bar{\omega}_{0}$
(7)

where $s_{n}$ corresponds to coefficients one obtains by expanding the sine of
an infinite series, i.e.,

$\sin\left(\sum_{n=0}^{\infty}a_{n}\bar{t}^{n}\right)=\sum_{n=0}^{\infty}s_{n}\bar{t}^{n}.$
(8)

The expansion (8) is obtained using the same method as that for raising a a
series to a power Henrici , taking the log of a series SIR2020 , and taking
the exponential of a series SEIR2020 . In fact, we make direct use of the
latter formula in what follows, and repeat its construction here. The
exponential of a series

$\text{exp}\left(\sum_{n=0}^{\infty}a_{n}\bar{t}^{n}\right)=\sum_{n=0}^{\infty}b_{n}\bar{t}^{n}$
(9)

is obtained by first noting that $b_{0}=e^{a_{0}}$ and then differentiating
both sides of (9) such that

$\left[\sum_{n=0}^{\infty}(n+1)~{}a_{n+1}~{}\bar{t}^{n}\right]\text{exp}\left(\sum_{n=0}^{\infty}a_{n}\bar{t}^{n}\right)=\sum_{n=0}^{\infty}(n+1)~{}b_{n+1}~{}\bar{t}^{n},$
(10)

which can be rewritten, using (9), as

$\left[\sum_{n=0}^{\infty}(n+1)~{}a_{n+1}~{}\bar{t}^{n}\right]\sum_{n=0}^{\infty}b_{n}\bar{t}^{n}=\sum_{n=0}^{\infty}(n+1)~{}b_{n+1}~{}\bar{t}^{n}.$
(11)

Cauchy’s product rule is then applied to the left-hand side of (11) to arrive
at

$\sum_{n=0}^{\infty}\left[\sum_{k=0}^{n}(k+1)~{}a_{k+1}~{}b_{n-k}\right]\bar{t}^{n}=\sum_{n=0}^{\infty}(n+1)~{}b_{n+1}~{}\bar{t}^{n},$
(12)

after which, equating like powers of $\bar{t}$ in (12) leads to a recursive
expression for the coefficients of the expansion of an exponential of an
infinite series, collected with (9) as

$\displaystyle\text{exp}\left(\sum_{n=0}^{\infty}a_{n}\bar{t}^{n}\right)=\sum_{n=0}^{\infty}b_{n}\bar{t}^{n},$
$\displaystyle
b_{n+1}=\frac{1}{(n+1)}\displaystyle\sum_{k=0}^{n}(k+1)~{}a_{k+1}~{}b_{n-k},~{}~{}b_{0}=e^{a_{0}}.$
(13)

To find the coefficients of (8), we now define the cosine of a series as

$\cos\left(\sum_{n=0}^{\infty}a_{n}\bar{t}^{n}\right)=\sum_{n=0}^{\infty}c_{n}\bar{t}^{n},$
(14)

and apply Euler’s identity

$\text{exp}\left(i\sum_{n=0}^{\infty}a_{n}\bar{t}^{n}\right)=\cos\left(\sum_{n=0}^{\infty}a_{n}\bar{t}^{n}\right)+i\sin\left(\sum_{n=0}^{\infty}a_{n}\bar{t}^{n}\right),$

directly to (13) (using the notation of (8) and (14)) to yield:

$c_{n+1}+i~{}s_{n+1}=\frac{1}{(n+1)}\displaystyle\sum_{k=0}^{n}(k+1)~{}i~{}a_{k+1}~{}\left(c_{n-k}+i~{}s_{n-k}\right).$
(15)

In (15), it is clear that $s_{0}=\sin(a_{0})$ and $c_{0}=\cos(a_{0})$ from (8)
and (14), respectively. Before continuing, it is useful to note that $a_{n}$,
$c_{n}$, and $s_{n}$ are real for all $n$. Now, equating the real part of each
side of (15) we obtain a recursion for the coefficients of (14). Similarly,
equating the imaginary part of each side of (15), we obtain a recursion for
the coefficients of (8). Collecting these (coupled) recursions with (7), we
may now write the series solution for the pendulum equation (1b) as

$\theta=\sum_{n=0}^{\infty}a_{n}\bar{t}^{n},~{}~{}|\bar{t}|<\bar{t}_{\text{ROC}}$
(16a) $\displaystyle a_{n+2}$
$\displaystyle=\frac{-s_{n}}{(n+1)(n+2)},~{}a_{0}=\theta_{0},~{}a_{1}=\bar{\omega}_{0},$
(16b) $\displaystyle s_{n+1}$
$\displaystyle=\frac{1}{n+1}\sum_{k=0}^{n}(k+1)a_{k+1}c_{n-k},~{}s_{0}=\sin(\theta_{0})$
(16c) $\displaystyle c_{n+1}$
$\displaystyle=\frac{-1}{n+1}\sum_{k=0}^{n}(k+1)a_{k+1}s_{n-k},~{}c_{0}=\cos(\theta_{0}).$
(16d)

where $\bar{t}_{\text{ROC}}$ is the radius of convergence of the series, to be
discussed in section III.2.

An alternative derivation of the recursion in series (16) is given in Fairen
by writing (II) as a coupled system of 4 ODEs in dependent variables $\theta$,
$d\theta/dt$, $\sin(\theta)$, and $\cos(\theta)$ and relating the coefficients
of their series expansions. As another alternative to the approach used here
to derive the series solution to (1b), one could substitute the well-known
series expansions Byrd ; Wrigge  of the Jacobi elliptic functions into (6b)
and (6c). However, an arcsine of the resulting series is then required to
extract the solution for $\theta(t)$, and for this reason, it is easier to
proceed as presented above.

A comparison of the numerical solution555All numerical solutions in this work
are generated using the 4th-order Runge-Kutta method with
$\Delta\bar{t}=10^{-5}$. Points on figures have been down-sampled for
readability. to (1b) and the $(N+1)$-term partial sums of the series solution
(16) are shown in figure 3 for both libration and rotation. Note that the
series diverges at a finite $\bar{t}_{\text{ROC}}$. In generating the
expansion (16) for figure 3, we chose initial conditions at the bottom of the
pendulum trajectory (i.e., $\theta_{0}=0$ in figure 3a). Figure 3 indicates
that for these choices of initial conditions, the pendulum series (16) cannot
fully describe the solution over the effective period $T^{*}$ due to
divergence. In the following section, we examine how $\bar{t}_{\text{ROC}}$ is
affected by initial conditions, and determine an optimal choice to ensure
convergence over the full range 0 to $T^{*}$. Note that such an exploration is
not restrictive, since the pendulum trajectory is invariant for a given choice
of $\bar{E}$ as discussed previously - the locus of all initial conditions
that achieve a given $\bar{E}$ yield trajectories that are simply shifted in
time but sweep through the same angles.

(a)  (b)
---
Figure 3: Series solution (16) (dashed lines) taken up to (and including) the
$\bar{t}^{N}$ term, shown for $N$=10, 102, and 103 and compared with the
numerical solution ($\circ$) to (1b) over 1 period. The solution is
initialized at the bottom of the pendulum’s trajectory for cases of (a)
libration ($\bar{E}<2$) with $\bar{E}=1.71$; and (b) rotation ($\bar{E}>2$)
with $\bar{E}=2.02$. Note that, for either case, the series does not converge
to the effective period $\bar{t}=T^{*}$ (indicated by a vertical line in the
plots).

###  III.2 Singularities of $\theta(t)$ and series convergence

Although any initial conditions $(\theta_{0},\omega_{0})$ may be chosen in
(16), the value of $\bar{t}_{\text{ROC}}$ is affected by this choice. Here, we
show that there is an optimal choice of $(\theta_{0},\omega_{0})$ such that
(16) converges beyond the effective period $(\bar{t}=T^{*})$ for all energy
levels $\bar{E}$ – and thus all possible motions – of the pendulum. This
enables the series solution (16) to be used as an exact solution that
describes the pendulum trajectory for all $\bar{t}$.

In order to determine the radius of convergence, $\bar{t}_{\text{ROC}}$, of
the pendulum series (16), the singularities of the pendulum function in the
complex-time plane are examined. The value of $\bar{t}_{\text{ROC}}$ is equal
to the distance between the expansion point, chosen to be $\bar{t}$=0 in (16),
and its closest singularity in the complex plane. For the case of the
separatrix orbit ($\bar{E}$=2), the arctan in (6a) admits logarithmic branch
point singularities, the smallest of which are at $\bar{t}=\pm i\pi/2$ (for
$\theta_{0}=0$ in (6a)), such that $\bar{t}_{\text{ROC}}=\pi/2$. However,
since an exact solution is known for the separatrix case in terms of
elementary functions, this restriction is illustrative but not relevant.
Unlike the separatrix case for which a closed-form solution exists, or for the
inverse solution $t(\theta)$ where the locations of singularities are known
from the integrands in the elliptic integrals involved (such as (4)), the time
singularities of $\theta(t)$ for $\bar{E}\neq 2$ are not obvious. As mentioned
in section I, it is typical that nonlinear ODEs such as (1b) contain
“spontaneous” singularities that are not discernible by inspection of the ODE
itself bender ; Barlow:2017 ; FS2020 ; SIR2020 ; SEIR2020 . In such problems,
an estimate for $\bar{t}_{\text{ROC}}$ may be obtained from numerical versions
of the root or ratio tests, the latter of which is often represented
graphically as a Domb-Sykes plot DombSykes .

Since the Jacobi elliptic functions have been well studied, their pole
locations are known exactly, and thus the singularities in the pendulum
function are also known as logarithmic branch points666The change of
singularity type from a pole to a logarithmic branch point is due to the
$\arcsin$ in (6b) and (6c). at the same locations through (6b) and (6c). This
enables an exact determination of $\bar{t}_{\text{ROC}}$ for the pendulum
series (16). For libration, the Jacobi elliptic function
$\operatorname{cd}(\bar{t},\sqrt{\bar{E}/2})$ appearing in (6b) has an
infinite number of poles Byrd  at the complex time values

$\bar{t}_{\mathrm{poles}}=n_{o}\underbrace{K\left(\sqrt{\frac{\bar{E}}{2}}\right)}_{T^{*}=T/4}+n_{o}^{\prime}~{}iK^{\prime}\left(\sqrt{\frac{\bar{E}}{2}}\right),\text{
for }~{}~{}\bar{E}<2,~{}\bar{\omega}_{0}=0$ (17)

where $n_{o}$ and $n_{o}^{\prime}$ are odd integers, $K(k)$ is given by (4),
and $K^{\prime}(k)$ is (in general) defined as

$K^{\prime}(k)\equiv K(\sqrt{1-k^{2}}).$ (18)

In (18), $K^{\prime}(k)$ is the magnitude of the “imaginary quarter-period”
chapter ; Appell . The singularities prescribed by (17) are shown in the
complex $\bar{t}$-plane in figure 4 for a constant value of $\bar{E}<2$; note
that the singularity structure relative to $\bar{t}=0$ is qualitatively the
same for all values of $\bar{E}<2$ so one figure suffices to draw general
conclusions. The series expansion about $t=0$ given by (16) encounters the
closest singularity when $n_{o}=n_{o}^{\prime}=1$ in (17) and the distance to
it (i.e., the radius of convergence), $\bar{t}_{\text{ROC}}$, is given as
follows:

$\bar{t}_{\text{ROC}}=\sqrt{{T^{*}}^{2}+\left[K^{\prime}\left(\sqrt{\frac{\bar{E}}{2}}\right)\right]^{2}}>T^{*}\text{;
valid for }~{}~{}\bar{E}<2,~{}\bar{\omega}_{0}=0.$ (19)

As noted in (19), the radius of convergence of (16) is always larger than the
effective period $T^{*}$ (defined in (5)) for the given initial conditions,
since the integral $K^{\prime}(k)$ is never 0. This circle of convergence is
shown as a solid curve in figure 4 centered around the expansion point
$\bar{t}$=0 with initial conditions at the top of the pendulum’s trajectory
($\theta_{0}\neq 0,\bar{\omega}_{0}=0$); this is an optimal expansion point
with respect to convergence for a given $\bar{E}<2$, as indicated in the
figure. If one were to expand about $\bar{t}=\pm T^{*}$ or, equivalently,
redefine $\bar{t}=0$ on the real-axis and choose initial conditions at the
bottom of the trajectory ($\theta_{0}=0,\bar{\omega}_{0}\neq 0$) for the same
$\bar{E}<2$, the radius of convergence of (16) is less than the effective
period $T^{*}$, as shown by the smaller dashed circles in figure 4; this is
observed in figure 3a, where the series diverges prior to reaching $t=T^{*}$.
Note that the radius of convergence decreases monotonically for a given $E<2$
as the initial conditions are moved from the top of the trajectory
($\theta_{0}\neq 0,\bar{\omega}_{0}=0$) to the bottom of the trajectory
($\theta_{0}=0,\bar{\omega}_{0}\neq 0$) for the case of a librating pendulum.
This statement is made in the context of Equation (2), where $\theta_{0}$ and
$\bar{\omega}_{0}$ are not independent quantities for a given $\bar{E}$ –
hence the omission of $\theta_{0}$ in (19).

Figure 4: Schematic showing periodic structure of the singularities (shown as
$\bullet$’s) of (6b) for the librating pendulum at a constant value of
$\bar{E}<2$, as specified by (17). The integrals $K(k)$ and $K^{\prime}(k)$
are defined by (4) and (18), respectively. The radius (19) and placement of
the solid circle corresponds to the original initial conditions given by (6b),
initializing the pendulum at the top of its trajectory with zero speed; this
is the optimal radius of convergence $\bar{t}_{\text{ROC}}$ for the series
solution (16) of a librating pendulum. By contrast, dashed circles show the
radius of convergence when a librating pendulum is initialized at the bottom
of the trajectory with nonzero angular speed.

For the case of rotation ($E>2$), a similar expression for pole locations is
available Byrd  for the Jacobi elliptic function
$\operatorname{sn}(\sqrt{\bar{E}/2}~{}\bar{t},\sqrt{2/\bar{E}})$ that appears
in (6c); here, the infinite number of singularities in the pendulum function
are placed in the complex $\bar{t}$-plane at

$\bar{t}_{\text{poles}}=n_{e}\underbrace{\tilde{K}\left(\sqrt{\frac{2}{\bar{E}}}\right)}_{T^{*}=T/2}+n_{o}^{\prime}~{}i\tilde{K}^{\prime}\left(\sqrt{\frac{2}{\bar{E}}}\right),\text{
for }~{}~{}\bar{E}>2,~{}\theta_{0}=0,$ (20)

where $n_{e}$ and $n_{o}^{\prime}$ respectively denote even and odd integers
and we define $\tilde{K}(k)\equiv(\sqrt{2/\bar{E}})~{}K(k)$ and
$\tilde{K^{\prime}}(k)\equiv(\sqrt{2/\bar{E}})~{}K^{\prime}(k)$ with
$K^{\prime}(k)$ given by (18). The singularities prescribed by (20) are shown
in the complex $\bar{t}$-plane in figure 5 for a constant value of
$\bar{E}>2$. The singularity structure for rotation is qualitatively similar
for all values of $\bar{E}>2$ so one figure suffices to draw general
conclusions. The circle of convergence for the series expansion (16) taken
about $t=0$, with initial conditions specified by (20) at the bottom of the
pendulum trajectory, is shown as the leftmost dashed circle in figure 5. A
circle of convergence of the same size is shown centered around $t=T^{*}$,
where the pendulum is also at the bottom of its trajectory. In both cases, the
radius of convergence777The precise value of the radius of the smaller circles
in figure 5 may be obtained by letting $n_{e}=0$ and $n_{0}^{\prime}=1$ in
(20) and taking its modulus; this radius is only larger than $T^{*}$ if
$\bar{E}>4$. Thus, the pendulum series for rotation cannot converge beyond
$\bar{t}=T^{*}$ for $\bar{E}\in[2,4]$ if initialized at the bottom of the
trajectory. is shown to be less than $T^{*}$, such that the series diverges
prior to reaching the effective period. This is observed in figure 3b, where
the series (16) is initialized with $\theta_{0}=2\pi$ and diverges prior to
reaching $t=T^{*}$. All other points initialized at the bottom of a rotating
pendulum’s trajectory are equally detrimental to convergence for a given
$\bar{E}>2$. As shown in figure 5, an optimal expansion point (or initial
condition) for a given energy $\bar{E}>2$ with respect to convergence is one
that is again taken at the top of the pendulum’s trajectory (e.g.,
$\theta_{0}=\pi,\bar{\omega}_{0}\neq 0$). Using the solid circle centered at
$t=T^{*}$ and the surrounding lattice in figure 5, the optimal radius of
convergence and initial conditions for rotation ($\bar{E}>2$) is given by

$\bar{t}_{\text{ROC}}=\sqrt{{T^{*}}^{2}+\left[\tilde{K}^{\prime}\left(\sqrt{\frac{2}{\bar{E}}}\right)\right]^{2}}>T^{*}\text{;
valid for }~{}~{}\bar{E}>2,~{}\theta_{0}=\pi,$ (21)

where again, the radius of convergence is always larger than the effective
period $T^{*}$, since $\tilde{K}^{\prime}(k)$ is nonzero for finite $\bar{E}$.
As a reminder when interpreting these equations, note that, by inspection of
(2), $\theta_{0}$ and $\bar{\omega}_{0}$ are not independent quantities for a
given $\bar{E}$, allowing for the omission of $\bar{\omega}_{0}$ in (21).

Figure 5: Schematic showing the periodic structure of the singularities (shown
as $\bullet$’s) of (6b) of (6c) for the rotating pendulum at a constant value
of $\bar{E}>2$, as specified by (20).
$\tilde{K}(k)\equiv\sqrt{2/\bar{E}}~{}K(k)$ and
$\tilde{K^{\prime}}(k)\equiv\sqrt{2/\bar{E}}~{}K^{\prime}(k)$, where the
integrals $K(k)$ and $K^{\prime}(k)$ are defined by (4) and (18),
respectively. The leftmost dashed circle corresponds to the original initial
conditions given by (6c), initializing the pendulum at the bottom of its
trajectory. The radius (21) and placement of the solid circle correspond to a
shifted initialization at the top of the pendulum’s trajectory, which is the
optimal radius of convergence $\bar{t}_{\text{ROC}}$ for the series solution
(16) of a rotating pendulum.

In both the case of libration and rotation, initializing the series solution
to the pendulum at the highest point along the trajectory leads to the largest
radius of convergence. This allows us to exploit the symmetry of the problem
such that the series solution (16), when prescribed initial conditions as
follows,

$\theta(\bar{t})=\sum_{n=0}^{\infty}a_{n}t^{n},~{}(\theta_{0},~{}\bar{\omega}_{0})=\begin{cases}\begin{rcases}\left(\arccos(1-\bar{E}),~{}0\right),&\text{$\bar{E}<2$}\\\
\left(\pi,~{}\pm\sqrt{2\bar{E}-4}\right),&\text{$\bar{E}>2$},\end{rcases},\end{cases},|\bar{t}|\leq
T^{\star}$ (22)

is an exact representation of the pendulum trajectory for all $t$ between 0
and $T^{*}$ and for any $\bar{E}$ describing librating and rotating motions of
the pendulum. In (22), the $a_{n}$ coefficients are given in (16), the $\pm$
symbol specifies either counterclockwise or clockwise rotation, and (2) has
been used to represent the initial conditions (taken at the top of the
trajectory) in terms of $\bar{E}$.

###  III.3 Periodic extension of the series solution

The convergent series solution (22), justified in section III.2 for
$\bar{t}\leq T^{*}$, allows for a convergent solution to (1b) for all
$\bar{t}$ as follows. For any given initial condition
$(\theta_{0},\bar{\omega}_{0}$), one first computes the energy $\bar{E}$ given
by (2). Using this value of $\bar{E}$, the convergent series for the solution
up to the effective period $T^{*}$ is constructed from (22) (or the
resummation (24d) provided in section IV) as

$\tilde{\theta}(\bar{t})\equiv\theta(\bar{t}\leq T^{*}).$ (23a) The solution
may then be extended for all $\bar{t}$ using the symmetry of the problem
described in section II, such that the periodic solution is given by
$\displaystyle\theta(\bar{t})=$
$\displaystyle\begin{cases}\begin{rcases}\tilde{\theta}(\hat{t}),&0\leq\hat{t}\leq
T^{*}\\\ -\tilde{\theta}(2T^{*}-\hat{t}),&T^{*}\leq\hat{t}\leq 2T^{*}\\\
-\tilde{\theta}(\hat{t}-2T^{*}),&2T^{*}\leq\hat{t}\leq 3T^{*}\\\
\tilde{\theta}(4T^{*}-\hat{t}),&3T^{*}\leq\hat{t}\leq
4T^{*}\end{rcases},&\bar{E}<2\\\ ~{}\\\
\begin{rcases}\pm\frac{2\pi}{T}(\bar{t}-\hat{t})+\tilde{\theta}(\hat{t}),&0\leq\hat{t}\leq
T^{*}\\\
\pm\frac{2\pi}{T}(\bar{t}-\hat{t})-\tilde{\theta}(2T^{*}-\hat{t})+2\pi\pm
2\pi,&T^{*}\leq\hat{t}\leq 2T^{*}\\\ \end{rcases},&\bar{E}>2\end{cases}$
$\displaystyle\text{ where }\hat{t}=t\text{ modulo }T,~{}~{}~{}~{}\bar{t}\geq
0.$ (23b)

As in (22), the $\pm$ symbol in (23b) specifies either counterclockwise or
clockwise rotation. The modulo operation888The modulo operation is available
in most programming languages. For example, the MATLAB syntax is
mod$(\bar{t},T)$. in (23b) is used such that $\hat{t}$ may be interpreted as
the additional time beyond the smallest integer multiple of the period $T$. In
section V, we provide an exact and rapidly converging series for the
computation of $T$ (and by extension, $T^{*}$), thus preserving the
analyticity of the solution when used in (23). The re-initializing and period
extension given by (23) is shown in figures 6 through 8. Note that the
formerly-divergent series solutions of figures 3a and 3b have been re-
initialized using (22) such that they are now convergent, and whose prior
domains are now contained within the periodically extended domains of figures
6 and 7. A final step in this solution process is to shift the origin such
that the originally intended initial conditions are at $\bar{t}=0$.

The norm of the error for $\bar{t}\in[0,T^{*}]$ incurred by series (22) is
plotted versus $\bar{E}$ (shown as dashed lines) in figure 9 for various
values of series truncation. Although the error increases as $\bar{E}\to 2$,
for any given $\bar{E}$, the error can be reduced to machine precision by
increasing the number of terms in the series. This is expected, since there
are no singularities of the pendulum function in the circle $|\bar{t}|<T^{*}$
for the initial conditions used in (22).

Figure 6: Series solution (23) (using (22) for $\tilde{\theta}(t)$) taken up
to (and including) the $\bar{t}^{N}$ term, shown for $N$=2, 4, and 6 and
compared with the numerical solution ($\bullet$) to (1b) for the case of
libration ($\bar{E}<2$). The energy $\bar{E}=1.71$ is same as that in figure
3a but here the series is initialized at the top of the trajectory such that
it converges. The solution is shown over 2 periods above, encompassing the
solution in figure 3a.  Figure 7: Series solution (23) (using (22) for
$\tilde{\theta}(t)$) taken up to (and including) the $\bar{t}^{N}$ term, shown
for $N$=3, 6, and 36 and compared with the numerical solution ($\bullet$) to
(1b) for the case of rotation ($\bar{E}>2$). The energy $\bar{E}=2.02$ is same
as that in figure 3b but here the series is initialized at the top of the
trajectory such that it converges. The solution is shown over 2 periods above,
encompassing the solution in figure 3b.

Figure 8: Series solution (23) (using (22) for $\tilde{\theta}(t)$) taken up
to (and including) the $\bar{t}^{N}$ term, shown for $N$=5, 10, and 20 and
compared with the numerical solution ($\bullet$) to (1b) for the librating
case of $\bar{E}=1.9998$.

##  IV Exact resummation of the series solution

Although (22) is an exact analytic solution to the pendulum problem, its
convergence may be accelerated by a similarly exact resummation technique. The
resummation is motivated by the method of asymptotic approximants Barlow:2017
that assures that it asymptotically approaches the exact solution at each end
of the domain, taken here to be $\bar{t}\in[0,T^{*}]$. We match the low-order
behavior at $\bar{t}=T^{*}$ (up to the 2nd derivative) and up any any desired
order at $t=0$. Using the choice of initial conditions prescribed by (22) to
preserve analyticity, the resummation of (22) is

$\displaystyle\theta(\bar{t})=\bar{\omega}^{*}\left(\bar{t}-T^{*}\right)+\left(\bar{t}-T^{*}\right)^{2}\sum_{n=0}^{\infty}\hat{a}_{n}\bar{t}^{n},~{}~{}|\bar{t}|\leq
T^{\star},$
$\displaystyle(\theta_{0},~{}\bar{\omega}_{0})=\begin{cases}\begin{rcases}\left(\arccos(1-\bar{E}),~{}0\right),&\text{$\bar{E}<2$}\\\
\left(\pi,~{}\pm\sqrt{2\bar{E}-4}\right),&\text{$\bar{E}>2$},\end{rcases},\end{cases}$
(24a) where $\bar{\omega}^{*}\equiv\begin{cases}\sqrt{2\bar{E}}&\bar{E}<2\\\
\pm\sqrt{2\bar{E}}&\bar{E}>2\end{cases}$ (24b) enforces the correct angular
velocity at $\bar{t}=T^{*}$ and the $\pm$ symbol specifies either
counterclockwise or clockwise rotation. In (24a), the coefficients
$\hat{a}_{n}$ are related to the known coefficients $a_{n}$ (given by (16b))
as follows. We set (24a) equal to (16a) to obtain
$\sum_{n=0}^{\infty}a_{n}\bar{t}^{n}=\bar{\omega}^{*}\left(t-T^{*}\right)+\left(t-T^{*}\right)^{2}\sum_{n=0}^{\infty}\hat{a}_{n}t^{n},$
and then isolate the $\hat{a}_{n}$ series onto one side to obtain
$\overbrace{\left[\sum_{n=0}^{\infty}\left((n+1)\left(\frac{1}{T^{*}}\right)^{n+2}\right)\bar{t}^{n}\right]}^{(\bar{t}-T^{*})^{-2}}\underbrace{\sum_{n=0}^{\infty}b_{n}\bar{t}^{n}}_{\left[-\bar{\omega}^{*}\left(t-T^{*}\right)+\displaystyle\sum_{n=0}^{\infty}a_{n}\bar{t}^{n}\right]}=\sum_{n=0}^{\infty}\hat{a}_{n}\bar{t}^{n}.$
As indicated in the above, $(\bar{t}-T^{*})^{-2}$ is expanded explicitly as an
arithmetico-geometric series, and the coefficients
$b_{0}=a_{0}+T^{*}\bar{\omega}^{*},~{}b_{1}=a_{1}-\bar{\omega}^{*},~{}b_{n>1}=a_{n}$
(24c) are introduced to enable the left-hand side of (24) to be combined into
a single Cauchy product series. The result is an expression for the
$\hat{a}_{n}$ coefficients as
$\hat{a}_{n}=\sum_{k=0}^{n}b_{n-k}\left((k+1)\left(\frac{1}{T^{*}}\right)^{k+2}\right).$
(24d)

Together, (24a), (24c), and (24d) form an exact resummation of (22) that
explicitly enforces the asymptotic behaviors as $\bar{t}\to 0$ (from the power
series (22)) and $\bar{t}\to T^{*}$, whereas (22) does not have this latter
feature. We remind the reader that the restriction on initial conditions in
(24d) is not a restriction on energy; one can use (24d) to provide convergent
solutions for all librating and rotating energy states of the pendulum.

Figure 9: The error here is taken to be the difference between the numerical
solution to (1b) and either the series solution (22) or exact resummation
(24d) (as indicated in the figure) taken up to (and including) the
$\bar{t}^{N}$ term and shown for various $N$. The infinity norm of this error
is taken over the interval $\bar{t}\in[0,T^{*}]$ for each value of $\bar{E}$
shown, except at $\bar{E}=2$ where an exact solution exists. The difference
between the exact solution (6a) and the numerical solution for $\bar{E}=2$ is
shown as a horizontal line, indicating the effective machine precision, as a
baseline of comparison with all the other error calculations shown (for which
the same numerical scheme and time-step is used).

Taking the resummation (24d) as $\tilde{\theta}(t)$ in the periodically
extended solution (23), the ($N+1$)-term partial sums are shown in figure 10.
Inspection of figure 8 and figure 10, generated for the same physical
parameters, indicates that the resummation (24d) accelerates convergence, in
that less terms are needed to achieve the same accuracy. This convergence
acceleration extends to other energy values, as shown in figure 9, where the
norm of the error in (24d) is shown (by a solid line) for various ($N+1$)-term
partial sums of the resummation series (24d).

Figure 10: Exact resummation (24d) taken up to (and including) the
$\bar{t}^{N}$ term, shown for $N$=5, 10, and 20, periodically extended using
(23), and compared with the numerical solution ($\bullet$) to (1b) for same
conditions of figure 8. A direct comparison with figure 8 shows the
convergence advantage imparted by the resummation series (24d) over the
original series (22).

Note that further manipulations, shown in Appendix A, may be used to increase
computational speed for any finite truncation of the resummation series (24d).

##  V Exact resummation of the pendulum period

In both the series solution (22) and its resummation (24d), the series
converge for $\bar{t}\in[0,T^{*}]$, and $T^{*}$ explicitly appears in both
(22) and the periodic extension formulae given in (23); thus, the value of
$T^{*}$ (and by extension, the value of $T$) must be determined accurately. As
can be seen in (3) and (5), in every case of period motion, $T^{*}$ is written
in terms of the complete elliptic integral of the first kind, $K(k)$, where
$k$ depends on the energy of the pendulum. Thus we are faced with computing
the period either by numerical integration, using an approximation (of which
there are many Hinrichsen ), or using a truncation of the long-known infinite
series representation of $K(k)$ Byrd ; Abramowitz , namely:

$K(k)=\frac{\pi}{2}\sum_{n=0}^{\infty}\left(\frac{(2n)!}{n!^{2}}\right)^{2}16^{-n}k^{2n},~{}~{}k^{2}<1.$
(25)

In (25), note that the $k^{2}<1$ constraint is always satisfied for the
pendulum problem (see (3), for example). The error in calculating $T^{*}$
using (25) in (5) is shown by the dashed lines in figure 11, plotted versus
$\bar{E}$ for partial sums up to and including the $k^{2N}$ term (i.e, $N+1$
terms in (25)). Note that the error increases as $\bar{E}\to 2$ where the
period becomes infinite, as it corresponds to the separatrix in figure 2.

Figure 11: The error is shown for evaluating a pendulum’s effective period
$T^{*}$ (and by extension, actual period) given by (5). In evaluating
$T^{*}_{\mathrm{series}}$, the elliptic integral $K(k)$ given by (4), is taken
to be either the original series (25) (dashed line) or its resummation (28)
(solid line) taken up to (and including) the $k^{2N}$ term. The exact
solution, denoted here as $T^{*}_{\mathrm{exact}}$, uses a numerical
evaluation of (4) that has converged to within machine precision.

Figure 12 shows results close to $\bar{E}=2$ (above and below) that show the
effect on the periodic extension of solution (24d) via (23) when (25) is used
to evaluate $T^{*}$ (see dashed curves). In the figure, a truncation of
$N$=100 is used in (25) in order in illustrate that one would need to take
more than 100 terms in the series (25) before a solution using this value of
$T^{*}$ can be considered accurate.

Figure 12:  Here, we demonstrate the effect of period accuracy on the
periodically extending (24d) via (23) (taken up to 20th order), shown for
librating and rotating cases near $\bar{E}=2$ as indicated in the figure. The
numerical solution is given by $\bullet$s. The long-dashed curves show a
solution that uses a period computed using (3) with (25) and letting $N=100$.
The short-dashed and solid curves show a solution that uses a period computed
using (3) with (28) and letting $N=1$ and $N=10$, respectively.

In order to make practical use of (25), we accelerate its convergence using a
resummation approach. To do so, we begin with an equivalent definition to (4)
for the complete elliptic integral of the first kind Byrd :

$K(k)=\int_{0}^{1}\left[(1-v^{2})(1-k^{2}v^{2})\right]^{-1/2}\ \text{d}v.$
(26)

We follow a similar approach to that used in Xue2012  to strip off the
singular behavior of (26) by first rewriting the integral as

$K(k)=\underbrace{\int_{0}^{1}\frac{\sqrt{1-k^{2}v^{2}}-\sqrt{1-v^{2}}}{(1-k^{2}v^{2})\sqrt{1-v^{2}}}\
\text{d}v}_{I_{1}}+\underbrace{\int_{0}^{1}\frac{1}{1-k^{2}v^{2}}\
\text{d}v}_{I_{2}}$

so that $I_{2}$ naturally handles the apparent integrand singularity at $k=1$.
Using partial fractions, this integral may be evaluated exactly as

$I_{2}=\frac{1}{2k}\ln\left(\frac{1+k}{1-k}\right),$

and Taylor expanded as

$I_{2}=\frac{1}{2k}\ln\left(\frac{1+k}{1-k}\right)=\frac{\operatorname{arctanh}(k)}{k}=\sum_{n=0}^{\infty}\frac{1}{2n+1}k^{2n}.$
(27)

Thus, we may use (27) to add 0 to (25) in such a way that the singular
behavior (responsible for slow convergence) is exactly summed, leading to an
improved series for $K(k)$, given as

$K(k)=\sum_{n=0}^{\infty}\left.\vphantom{\frac{\pi}{2}\left(\frac{(2n)!}{n!^{2}}\right)^{2}16^{-n}}\right[\frac{\pi}{2}\left(\frac{(2n)!}{n!^{2}}\right)^{2}16^{-n}-\frac{1}{2n+1}\left.\vphantom{\frac{\pi}{2}\left(\frac{(2n)!}{n!^{2}}\right)^{2}16^{-n}}\right]k^{2n}+\frac{1}{2k}\ln\left(\frac{1+k}{1-k}\right),~{}~{}k^{2}<1.$
(28)

When used to compute $T^{*}$ (or $T$) in (5) (or (3)) the convergence of the
resummation (28) is significantly improved, compared with that of (25), as
seen in figure 11 for a range of $\bar{E}$ and in figure 12 for $\bar{E}$
close to 2.

##  VI Conclusions

The exact series solution to the pendulum problem, whether librating or
rotating, attains optimal convergence when the initial conditions are chosen
to be at the top of the pendulum’s trajectory, allowing the series to converge
beyond the bottom of the trajectory. Through symmetry, this convergent
solution is then used to construct the solution for all time. Taking
inspiration from the method of asymptotic approximants, the exact power series
solution is re-summed to accelerate convergence, based on the asymptotic
behavior at the bottom of the trajectory. The resulting exact power series
representation is shown to be more accurate than the original series for any
given series truncation and pendulum energy.

We close this paper with some general comments about the methodology used here
as well as the solution itself. When finding the series solution to a
nonlinear ODE, there are (at least) 3 scenarios: (i) the radius of convergence
is larger than the physical domain, (ii) the radius of convergence is
precisely the physical domain length (e.g., series for Rayleigh’s equation
Harkin ), and (iii) the radius of convergence lies within the physical domain
(e.g, Blasius series Blasius ). The exact pendulum solution provided in this
paper falls into case (i), once the initial conditions are chosen with care.
Here our analysis benefited from the properties of the Jacobi elliptic
functions whose singularity locations are known. In case (ii), a singularity
is known at the boundary of the domain and this is the closest singularity in
the complex plane to the expansion point. Case (iii) arises from spontaneous
singularities and is commonplace in the solution of nonlinear ODEs. Note that
even when case (iii) occurs, resummation methods such as asymptotic
approximants may be used to analytically continue a series over the whole
domain to obtain a useful convergent solution (see Barlow:2017 ; FS2020 ;
SIR2020 ; SEIR2020 ).

In contrast to the problem examined here, there is often no pre-knowledge of
singularity locations in the solution of nonlinear ODEs. For such cases, one
can find the radius of convergence for a series around a given expansion point
numerically using a Domb-Sykes plot DombSykes  or root test. For autonomous
equations such as this one, then, an optimal expansion point would need to be
deduced numerically; whether or not this expansion point leads an exact
solution depends on the final radius of convergence achieved. Note, however,
that even in cases where the radius of convergence is not larger than the
domain of interest, approximant techniques may be utilized to analytically
continue such series and make them useful. Coupled with approximants (when
needed), the power series solution method is a prevailing technique that can
lead to efficient, useful, and many times exact solutions to nonlinear ODEs.

##  Appendix A Efficient computation of (24d)

Recognizing that the exact resummation of the pendulum series, given by (24d),
will be implemented in any code as a partial sum, we provide an efficient
means of calculating the deflection angle. One may observe that the truncation
of (24a), denoted here as $\theta_{N}(t)$, has the following structure

$\theta_{N}(t)=\underbrace{\bar{\omega}^{*}\left(t-T^{*}\right)}_{\text{degree
1 polynomial}}+\underbrace{\left(t-T^{*}\right)^{2}}_{\text{degree 2
polynomial}}\underbrace{\hphantom{\int}\sum_{n=0}^{N}\hat{a}_{n}t^{n}}_{\text{degee
$N$ polynomial}}.$ (29a) Since the approximant given by (29a) is just a degree
$N+2$ polynomial, we may rewrite it as
$\theta_{N}(t)=\sum_{n=0}^{N}a_{n}t^{n}+\alpha t^{N+1}+\beta t^{N+2}.$ (29b)
Upon equating the above to the original series (16a), one finds that the
$a_{n}$ coefficients are precisely the original series coefficients given in
(16). Hence, in order to obtain the approximant in the form of (29b), we first
determine the $a_{n}$ from the formulae give in (16), and then use those to
compute $\alpha$ and $\beta$ so that the first and second derivatives of (29a)
and (29b) match at $\bar{t}=T^{*}$. This leads to a system of two linear
equations in $\alpha$ and $\beta$, whose solution is

$\displaystyle\alpha$
$\displaystyle=-(N+2)\sigma_{N}({T^{*}}){T^{*}}^{-N-1}-\left(\bar{\omega}^{*}-\sigma_{N}^{\prime}({T^{*}})\right){T^{*}}^{-N}$
(29c) $\displaystyle\beta$
$\displaystyle=(N+1)\sigma_{N}({T^{*}}){T^{*}}^{-N-2}+\left(\bar{\omega}^{*}-\sigma_{N}^{\prime}({T^{*}})\right){T^{*}}^{-N-1}$
(29d)

where

$\sigma_{N}(t)\equiv\sum_{n=0}^{N}a_{n}t^{N}.$ (29e)

For comparison, the entire process to compute the resummation coefficients
(for a given truncation $N$) via the formulae given in (24d) (including the
computation of series coefficients $a_{n}$ from (16)) requires two quadratic
time steps. The first step is the recursion in (16) used to get the series
solution, and the second is Cauchy’s product rule, used in (24d). The method
given in this appendix - recasting (24d) as (29e) - removes the latter of
those quadratic time steps, which leads to a substantial reduction in
computational time - such as when the series is re-used for different values
of $\bar{E}$ to construct the solution surface in figure 2.

## References

  * [1] G. F. Simmons.  Differential Equations with Applications and Historical Notes.  McGraw-Hill, 2 edition, 1994.
  * [2] S. H. Strogatz.  Nonlinear Dynamics and Chaos.  Westview, 1994.
  * [3] H. von Helmholtz.  Vorlesungen über theoretische Physik.  Leipzig, 1898.
  * [4] P. Appell.  Sur une interprétation des valeurs imaginaires du temps en mécanique.  Comptes. Rendus., 87(1):1074–1077, 1878.
  * [5] E.T. Whittaker.  A Treatise on the Analytical Dynamics of Particles and Rigid Bodies; With an Introduction to the Problem of Three Bodies.  Cambridge University Press, 2 edition, 1917.
  * [6] J. Stillwell.  Mathematics and its History.  Springer, 2 edition, 2002.
  * [7] J. Gray.  The Real and the Complex: A History of Analysis in the 19th Century.  Springer, 2015.
  * [8] P. F. Byrd and M. D. Friedman.  Handbook of Elliptic Integrals for Engineers and Scientists.  Springer-Verlag, 2 edition, 1971.
  * [9] R. P. Agnew.  Differential Equations.  McGraw-Hill, 1960.
  * [10] V. Fairén, V. López, and L. Conde.  Power series approximation to solutions of nonlinear systems of differential equation.  Am. J. Phys, 56(1):57–61, 1988.
  * [11] J. P. Boyd.  Padé approximant algorithm for solving nonlinear ordinary differential equation boundary value problems on an unbounded domain.  Comp. Phys., 11:299–303, 1997.
  * [12] N. S. Barlow, C. R. Stanton, N Hill, S. J. Weinstein, and A. G. Cio.  On the summation of divergent, truncated, and underspecified power series via asymptotic approximants.  Q. J. Mech. Appl. Math., 70(1):21–48, 2017.
  * [13] E. R. Belden, Z. A. Dickman, S. J. Weinstein, A. D. Archibee, E. Burroughs, and N. S. Barlow.  Asymptotic approximant for the falkner-skan boundary-layer equation.  Q. J. Mech. Appl. Math, 73(1):36–50, 2020.
  * [14] N. S. Barlow and S. J. Weinstein.  Accurate closed-form solution of the SIR epidemic model.  Physica D, 408:132540:1–4, 2020.
  * [15] S. J. Weinstein, M. S. Holland, K. E. Rogers, and N. S. Barlow.  Analytic solution of the SEIR epidemic model via asymptotic approximant.  Physica D, 411:132633:1–6, 2020.
  * [16] A. Harkin, A. Giammarese, N. S. Barlow, and S. J. Weinstein.  The rayleigh collapse of two spherical bubbles.  arXiv:2102.05222.
  * [17] E. Ramé, S. J. Weinstein, and N. S. Barlow.  Free surface shapes in rigid body rotation: Exact solutions, asymptotics and approximants.  IMA J. Appl. Math., doi:10.1093/imamat/hxab038, 2021.
  * [18] F. M. S. Lima.  Analytical study of the critical behavior of the nonlinear pendulum.  Am. J. Phys., 78(11):1146–1151, 2010.
  * [19] K. Ochs.  A comprehensive analytical solution of the nonlinear pendulum.  Eur. J. Phys., 32:479–490, 2011.
  * [20] W. A. Schwalm.  Lectures on Selected Topics in Mathematical Physics: Elliptic Functions and Elliptic Integrals.  Morgan & Claypool, 2015.
  * [21] P. Henrici.  Automatic computations with power series.  JACM, 3:10–15, 1956.
  * [22] S. Wrigge.  Calculation of the Taylor series expansion coefficients of the Jacobian elliptic function sn($x$, $k$).  Math. Comp., 36(154):555–564, 1981.
  * [23] C. M. Bender and S. A. Orszag.  Advanced Mathematical Methods for Scientists and Engineers I: Asymptotic Methods and Perturbation Theory.  McGraw-Hill, 1978.
  * [24] C. Domb and M. F. Sykes.  On the susceptibility of a ferromagnetic above the curie point.  Proc. Roy. Soc. London A, 240(1221):214–228, 1957.
  * [25] P. F. Hinrichsen.  Review of approximate equations for the pendulum period.  Eur. J. Phys., 42:015005:1–19, 2021.
  * [26] M. Abramowitz and I Stegun.  Handbook of Mathematical Functions, page 298.  Dover, 1972.
  * [27] De-Sheng Xue, Zhao Zhou, and Mei-Zhen Gao.  Accurate period approximation for any simple pendulum amplitude.  Chinese Physics Letters, 29(4):044601, April 2012.
  * [28] H. Blasius.  Grenzschichten in flussigkeiten mit kleiner reibung.  Zeitschrift fur Mathematik und Physik, 56:1–37, 1908.


11institutetext: Department of Computer Science and Engineering,
Indian Institute of Technology Jammu, Jammu & Kashmir 181221, India.
11email<EMAIL_ADDRESS><EMAIL_ADDRESS>

# Envy-free Trip Planning in Group Trip
Planning Query Problem ††thanks: Both the authors have contributed equally in
this work.

Mayank Singhal 11 Suman Banerjee 11

###### Abstract

In recent times, Group Trip Planning Query (henceforth referred to as GTP
Query) is one of the well-studied problems in Spatial Databases. The inputs to
the problem are a road network where the vertices represent the Point-of-
Interests (mentioned as POIs henceforth) and they are grouped into different
categories, edges represent the road segments, and edge weight represents the
distance and a group of users along with their source and destination
location. This problem asks to return one POI from every category such that
the aggregated distance traveled by the group is minimized. As the objective
is to minimize the aggregated distance, the existing solution methodologies do
not consider the individual distances traveled by the group members. To
address this issue, we introduce and study the Envy Free Group Trip Planning
Query Problem. Along with the inputs of the GTP Query Problem, in this
variant, we also have a threshold distance $D$ such that aggregated distance
traveled by the group is minimized and for any member pairs the difference
between their individual distance traveled is less than equal to $D$. However,
it may so happen that a given $D$ value no such set POIs are found. To tackle
this issue, we introduce the surrogate problem Envy Free Group Trip Planning
Query with Minimum Additional Distance Problem which asks what is the minimum
distance to be added with $D$ to obtain at least one solution. For these
problems, we design efficient solution approaches and experiment with real-
world datasets. From the experiments, we observe that the proposed solution
approaches lead to less aggregated distance compared to baseline methods with
reasonable computational overhead.

###### Keywords:

Spatial Database Trip Planning Query Road Network.

##  1 Introduction

Due to the recent advances of mobile devices and wireless networks, storing
the location information of moving objects become easy. We call them as
trajectory datasets and they are freely available in different data
repositories [8]. A spatial database contains the location information of
different POIs of a large geographical area (e.g. city) and provides services
to different kinds of user queries [3]. One of the well-studied problems in
the domain of spatial databases is the Group Trip Planning Query [2, 7, 4].
The inputs to the problem are a road network where the vertices represent the
Point-of-Interests (mentioned as POIs henceforth) and they are grouped into
different categories, edges represent the road segments, and edge weight
represents the distance and a group of users along with their source and
destination location. This problem asks to return one POI from every category
such that the aggregated distance traveled by the group is minimized [6]. As
the objective is to minimize the aggregated distance, the existing solution
methodologies do not consider the individual distances traveled by the group
members. It may so happen that the individual distance traveled by a group
member is much more than the other one. In such case, the previous group
member will be envious to the second one and such solutions even may not be
accepted in a realistic situation. To tackle this problem, in this paper we
introduce the notion of Envy-free ness in the GTP Query Problem and call this
problem as the Envy-free Group Trip Planning Query Problem.

In this new variant, along with the inputs of the GTP Query Problem, we also
have a threshold distance $D$ and the goal is to find one POI from each
category such that the aggregated distance is minimized as well as for any
pair of group members, the difference between the individual distances
traveled by them is less than equal to $D$. In practical situations, it is
desirable to keep the value of $D$ as much low as possible. However, if we
choose the value of $D$ too low then there may be no feasible solution at all.
To address this issue, we formulate the Envy-free GTP Query with Minimum
Additional Distance Problem. Here the question is how much minimum distance to
be added with $D$ such that the feasible solution set for the Envy-free GTP
Query Problem will be non-empty. In this paper, we study both these problems
and, in particular, we make the following contributions in this paper:

  * •

We incorporate the concept of Envy-freeness in the GTP Query Problem and
introduce the Envy-Free Group Trip Planning Query Problem.

  * •

To take care of the situation when the feasible solution set of the Envy-Free
GTP Query Problem is empty, we study the Envy-Free GTP Query with Minimum
Additional Distance Problem.

  * •

We design a solution methodology for solving these problems along with a
detailed analysis. Furthermore, we also explore the situation when the spatial
data is indexed with $\mathcal{R}$-Tree

  * •

We conduct an extensive set of experiments with three publicly available
trajectory datasets to show the effectiveness and efficiency of the proposed
solution approaches.

Rest of the paper of is organized as follows. Section 2 describes the
background and defines the problem formally. Section 3 describes the proposed
solution methodologies. Section 4 contains the experimental evaluation of the
proposed solution approaches. Section 5 concludes our study and gives future
research directions.

##  2 Background and Problem Definition

In this section, we give the background and defines the problem formally.
Initially, we start by defining the road network.

######  Definition 1 (Road Network)

A road network is denoted by weighted graph
$\mathcal{G}(\mathcal{V},\mathcal{E},\mathcal{W})$ where the vertex set
$\mathcal{V}(\mathcal{G})=\\{v_{1},v_{2},\ldots,v_{n}\\}$ are the set of $n$
POIs and the edge set
$\mathcal{E}(\mathcal{G})=\\{e_{1},e_{2},\ldots,e_{m}\\}$ are the $m$ road
segments connecting the POIs. The edge weight function $\mathcal{W}$ assigns
each edge to a real number which is basically the corresponding distance
between the two POIs of the edge; i.e.;
$\mathcal{W}:\mathcal{E}(\mathcal{G})\longrightarrow\mathbb{R}^{+}$.

We denote the number of vertices and edges of the road network by $n$ and $m$,
respectively. For any edge $(v_{i}v_{j})\in\mathcal{E}(\mathcal{G})$, For any
two POIs $v_{i}$ and $v_{j}$, we denote its weight by
$\mathcal{W}(v_{i},v_{j})$. For any two POIs $v_{i}$ and $v_{j}$,
$dist(v_{i},v_{j})$ denotes the shortest path distance between $v_{i}$ and
$v_{j}$. A path between the POIs $v_{i}$ and $v_{j}$, is a sequence of
vertices $<v_{i},v_{i+1},v_{i+2},\ldots,v_{j}>$ such that for all
$p\in\\{i,i+1,\ldots,j-1\\}$, $(v_{p}v_{p+1})\in\mathcal{E}(G)$. Also, we
assume that the road network is connected that means between any two POIs
there exist a path in $\mathcal{G}$. For any two POIs $v_{i}$ and $v_{j}$ the
shortest path distance is denoted by $dist(v_{i},v_{j})$. Most of the symbols
and notations related to graph theory has been taken from [1].

Now, a group of $b$ friends $\mathcal{U}=\\{u_{1},u_{2},\ldots,u_{b}\\}$ wants
to plan a trip in the city with the road network $\mathcal{G}$. Their starting
and destination locations are
$\mathcal{U}_{\mathcal{S}}=\\{v^{s}_{1},v^{s}_{2},\ldots,v^{s}_{b}\\}$ and
$\mathcal{U}_{\mathcal{D}}=\\{v^{d}_{1},v^{d}_{2},\ldots,v^{d}_{b}\\}$,
respectively. They want to travel from their respective starting locations
through different POIs. As an example, they may decide to first meet at a
coffee shop, then to watch movie in a movie theater, subsequently in a hotel
for lunch, then in a park for gossiping, and finally to a icecream parlour
before heading towards their respective destinations. It is quite natural that
in a large city there are multiple locations of each category. Now, which one
POI from each category to choose such that total aggregated distance traveled
is minimized? This problem is formally referred to as the Group Trip Planning
Query Problem.

######  Definition 2 (Group Trip Planning Query Problem)

Given a road network $\mathcal{G}(\mathcal{V},\mathcal{E},\mathcal{W})$ where
the vertices are marked into $k$ different groups $\mathcal{V}_{1}$,
$\mathcal{V}_{2}$, $\ldots$, $\mathcal{V}_{k}$, a group of $b$ friends along
with their source and destination locations, this problem asks to return one
POI for every category such that the total aggregated distance traveled by all
the group members is minimized.

Assume that the POIs $v_{1}$, $v_{2}$, $\ldots$, $v_{k}$ constitute the
solution for the GTP Query Problem. Now, it can be observed that the total
aggregated distance traveled by the group $\mathcal{U}$ is given by the
Equation 1.

$\mathcal{D}_{\mathcal{U}}=\underset{i\in[b]}{\sum}\ dist(v^{s}_{i},v_{1})\ +\
b\cdot\underset{i\in[k-1]}{\sum}\
dist(v_{i},v_{i+1})+\underset{i\in[b]}{\sum}\ dist(v_{k},v^{d}_{i})$ (1)

For any $i\in\mathbb{Z}^{+}$, $[i]$ denotes the set $\\{1,2,\ldots,i\\}$. The
goal here is the $v_{i}$s such that $\mathcal{D}_{\mathcal{U}}$ is minimized.
Hence, mathematically this problem can be expressed as follows.

$(v^{*}_{1},v^{*}_{2},\ldots,v^{*}_{k})=\underset{v_{1}\subseteq\mathcal{V}_{1},v_{2}\subseteq\mathcal{V}_{2},\ldots,v_{k}\subseteq\mathcal{V}_{k}}{argmin}\
\mathcal{D}_{\mathcal{U}}$

Now, it can be observed that this problem is only concerned with the
minimization of $\mathcal{D}_{\mathcal{U}}$ not the individual distance
traveled by a group member. So, it may so happen that as per the recommended
POI the distances traveled by only $u_{i}$ (denoted by $\mathcal{D}_{u_{i}}$)
is much more than that of $u_{j}$. In this situation, $u_{i}$ may become
envious towards $v_{j}$ which is not at all desirable. This motivates us to
introduce and study the Envy-free Group Trip Planning Query Problem which is
stated in Definition 3.

######  Definition 3 (Envy-Free Group Trip Planning Query Problem)

Given a road network $\mathcal{G}(\mathcal{V},\mathcal{E},\mathcal{W})$ where
the vertices are marked into $k$ different categories $\mathcal{V}_{1}$,
$\mathcal{V}_{2}$, $\ldots$, $\mathcal{V}_{k}$, a group of $b$ friends along
with their source and destination locations, and a threshold distance $D$,
this problem asks to select one POI from each category such that the
aggregated distance is minimized as well as for any pair of members the
individual distance traveled is bounded by the given threshold distance.

From the computational perspective, this problem can be written as follows:

Envy-free group trip planning query problem Input: Road Network
$\mathcal{G}(\mathcal{V},\mathcal{E},\mathcal{W})$; The Group
$\mathcal{U}=\\{u_{1},u_{2},\ldots,u_{b}\\}$ with their source and destination
locations respectively, and a threshold distance $D$. Problem: Return the POIs
$v^{*}_{1}\in\mathcal{V}_{1}$, $v^{*}_{2}\in\mathcal{V}_{2}$, $\ldots$,
$v^{*}_{k}\in\mathcal{V}_{k}$ such that the the quantity
$\mathcal{D}_{\mathcal{U}}$ is minimized and also for all
$u_{i},u_{j}\in\mathcal{U}$,
$|\mathcal{D}_{u_{i}}-\mathcal{D}_{u_{j}}|\leq\mathcal{D}$.

Now, its an important question how to choose the value of $D$. It is easy to
observe that if we choose a low value of $D$, then it may so happen that we do
not get any solution. So, the immediately next problem is the how much minimum
distance to be added with $D$ to obtain a solution for the Envy-free GTP Query
Problem. Formally, we call this problem as the Envy-free GTP Query with
Minimum Additional Distance Problem. It is important to observe that the
second problem comes into existence only when the feasible solution set of the
first problem becomes empty. Hence, we call the second problem as the
surrogate of the first problem.

##  3 Proposed Solution Approaches

In this section, we describe the proposed solution methodologies for the Envy-
free GTP Query Problem. Our first approach is based on the exhaustive search
method which has been described in the following subsection.

###  3.1 Exhaustive Search Method

In this method, we exhaustively search among all category of POIs to check
whether their exist at least one feasible solution or not. If there exists
more than one feasible solution then the one with the least aggregated
distance will be returned. If the feasible solution set is empty then the
Envy-free GTP Query with Minimum Additional Distance problem is solved. It is
important to observe that the number of possible combinations of the POIs can
be of $\mathcal{O}(\underset{i\in[k]}{\prod}\ n_{i})$ where $n_{i}$ denotes
the number of POIs of $i$-th category. The set of all possible POI
Combinations is denoted by $\mathcal{Q}$. Now, the worst case of Algorithm 1
will arise when the value of this quantity is maximum.Lemma 1 tells the
criteria of maximizing this quantity. Due to space limitation we are unable to
give the proves.

###### Lemma 1

Let $f:\mathbb{R}^{k}\longrightarrow\mathbb{R}$ be a function over a
$k$-dimensional real space and
$f(x_{1},x_{2},\ldots,x_{k})=\underset{i\in[k]}{\prod}\ x_{i}$. This function
attains its maximum value subject to the constraint
$\underset{i\in[k]}{\sum}x_{i}\leq n$ when
$x_{1}=x_{2}=\ldots=x_{k}=\frac{n}{k}$.

In particular, Lemma 1 talks about the distribution of the POIs among the
categories and it shows that the worst case of Algorithm 1 will occur when
POIs are equally distributed, and in that case the number of combinations will
be of $\mathcal{O}(n^{k})$. In realistic situations, no. of categories of POIs
is a constant, and hence, in this study also, we consider that $k$ is a
constant.

Now, we describe the working principle of Algorithm 1 maintains one two
dimensional array $Route\\_Matrix$ which is a two dimensional array of size
$(\underset{i\in[k]}{\prod}n_{i})\times(k+3)$ and its content is as follows:

  * •

Corresponding to every POI combination (and certainly there can be at most
$\underset{i\in[k]}{\prod}n_{i}$ many) there will be a row in the array
$Route\\_Matrix$.

  * •

In each row, among $(k+3)$ many cells, the first $k$ many cells contains the
POI combination corresponding to that row.

  * •

$(k+1)$-th cell contains the aggregated distance of the group considering the
POI Combinations as solution.

  * •

$(k+2)$-th entry stores the maximum difference between the individual
distances traveled among the all possible group member pairs.

$d^{\rho}_{max}=\underset{(u_{i}u_{j})\in\binom{\mathcal{U}}{2}}{max}\
|\mathcal{D}^{\rho}_{u_{i}}-\mathcal{D}^{\rho}_{u_{j}}|$ (2)

Here, $\rho$ denotes arbitrary POI Combination and $\binom{\mathcal{U}}{2}$
denotes the $2$-element subsets of $\mathcal{U}$.

  * •

$(k+3)$-th entry stores only a bit. If this entry is $1$, then it signifies
that the corresponding POI combination ( which is stored in first $k$ cells in
the row) is a feasible solution, else it is not a feasible solution.

Input: $\text{The Road Network
}\mathcal{G}(\mathcal{V},\mathcal{E},\mathcal{W})$, Threshold Distance $D$.

Output: Optimal

1 $Route\\_Matrix=Create\\_Array(\underset{i\in[k]}{\prod}\ n_{i},k+3,0)$

2

3for  _$\text{All }v_{1}\in\mathcal{V}_{1}$_ do

4 for  _$\text{All }v_{2}\in\mathcal{V}_{2}$_ do

5 $\ldots$

6 for  _$\text{All }v_{k}\in\mathcal{V}_{k}$_ do

7 $\text{The current POIs are }\rho=<v_{1},v_{2},\ldots,v_{k}>$

8 $DIS=Create\\_Array(b,0)$

9 for  _$\text{All }u\in\mathcal{U}$_ do

10
$DIS(u)=dist(u^{s}_{1},v_{1})+\underset{i\in[k-1]}{\sum}dist(v_{i},v_{i+1})+dist(v_{k},u^{d}_{1})$

11

12 $Aggr\\_Dis\longleftarrow\underset{u\in\mathcal{U}}{\sum}\ DIS(u)$

13 $Route\\_Matrix\longleftarrow(<v_{1},v_{2},\ldots,v_{k}>;Aggr\\_Dis;--;--)$

14 $Count\leftarrow 0$; $max\longleftarrow-1$

15 for  _$\text{All }(u_{i}u_{j})\in\binom{\mathcal{U}}{2}$_ do

16 if  _$|DIS(u_{i})-DIS(u_{j})|\leq D$_ then

17 $Count\leftarrow Count+1$

18

19 if  _$|DIS(u_{i})-DIS(u_{j})| >max$_ then

20 $max\longleftarrow|DIS(u_{i})-DIS(u_{j})|$

21

22

23 $Route\\_Matrix\longleftarrow(--;--;max;--)$

24 if  _$Count==\binom{b}{2}$_ then

25 $Route\\_Matrix\longleftarrow(--;--;--;1)$

26

27

28 $\ldots$

29

30

31$flag\longleftarrow FALSE$

32 for  _$i=1\text{ to }|Route\\_Matrix|$_ do

33 if  _$Route\\_Matrix[i,k+3]==1$_ then

34 $flag\longleftarrow TRUE$

35 $break$

36

37

38if  _$flag==TRUE$_ then

39 $<v^{*}_{1},v^{*}_{2},\ldots,v^{*}_{k}>=argmin\ (Route\\_Matrix[;*;;])$

40 $\text{return }<v^{*}_{1},v^{*}_{2},\ldots,v^{*}_{k}>$

41

42else

43 $d\longleftarrow min(Route\\_Matrix[;;*;])$

44 $\epsilon\longleftarrow d-D$

45 $\text{return ``No feasible solution. However, }\epsilon\text{ is added to
}D\text{ to have a solution}$.

46

Algorithm 1 Construction of the optimal solution of the Envy Free GTP Query
Problem from the GTP Query Problem

Now, for every combination POIs, we perform the following operations. We
compute the individual distances traveled by the group members which requires
$\mathcal{O}(b+k)\approx\mathcal{O}(b)$ time. Next, we compute the aggregated
distance for the current POI Combinations and stores it along with the current
POI combinations in the array $Route\\_Matrix$. This step requires
$\mathcal{O}(b)$ time. Initializing two variables count and max requires
$\mathcal{O}(1)$ time. As there are $b$ many members hence, the number of
member pairs will be of $\mathcal{O}(b^{2})$. Now, for each member pair the
time requirement for checking the envy-free ness constraint, updating the
count, condition checking of the if statement in Line No. $16$, and
subsequently updating the max will take $\mathcal{O}(1)$ time. Hence, time
requirement from Line No. $13$ to $17$ will take $\mathcal{O}(b^{2})$ time.
Inserting max into its the $Route\\_Matrix$ array will take $\mathcal{O}(1)$
time. It is important to observe that for one combination of POIs if for all
member pairs, the envy-free ness constraint is satisfied then only we can say
that the current POI combination is a feasible solution. We have incorporated
this in the following way. We maintain a count variable and initialized it
with zero. Now, we check the envy-free ness constraint for every user pairs
one by one and we keep on incrementing count by one. A the end, we check
whether the value of count is equals to $\binom{b}{2}$ or not. If it is so,
then it clearly signifies that for all the member pairs, the envy-free ness
constraint is satisfied. This implies that the current POI combination is a
feasible solution for the Envy-free GTP Query Problem, and hence, $(k+3)$-th
entry of the current row of $Route\\_Matrix$ array is set to $1$. Now, we can
easily observe that for one POI combination, time requirement is of
$\mathcal{O}(b^{2})$. Hence, total time requirement from Line No. $2$ to $21$
is of $\mathcal{O}(n^{k}b^{2})$.

In the $Route\\_Matrix$ array, if we find at least one row such that its
$(k+3)$-th entry is $1$ then the feasible solution set of the input instance
of the Envy-free GTP Query Problem is non-empty. This has been stated from
Line No. $22$ to $25$. First, we set a flag to FALSE. Next, we scan through
the rows of the $Route\\_Matrix$ array, and if we find one row whose
$(k+3)$-th entry is $1$ then we update the flag to TRUE and breaking out from
the loop. It is easy to observe that the time requirement for this step is of
$\mathcal{O}(n^{k})$. After this step if the flag is TRUE that means the
feasible solution set is non-empty. In that case we need to only scan the
$(k+1)$-th column of those rows of the $Route\\_Matrix$ array whose $(k+3)$-th
entry is $1$ and returns the POI Combination corresponding to the minimum
value of the $(k+1)$-th column. It is easy to observe that the time
requirement for this step is of $\mathcal{O}(n^{k})$ time. If the flag remains
FALSE after Line No. $26$ then there does not exist any feasible solution for
the given instance of the Envy-free GTP Query Problem. In that case, as
mentioned previously, Envy-free GTP Query Problem with Minimum Additional
Distance needs to be solved. For this, we find out minimum value among
$(k+2)$-th entries of all the rows of the $Route\\_Matrix$ array. Assume that
this minimum value is $d$. Mathematically, $d$ can be expressed as follows:

$d=\underset{<v_{1},v_{2},\ldots,v_{k}>\in\mathcal{Q}}{min}\ \
\underset{(u_{i}u_{j})\in\binom{\mathcal{U}}{2}}{max}\
|\mathcal{D}^{<v_{1},v_{2},\ldots,v_{k}>}_{u_{i}}-\mathcal{D}^{<v_{1},v_{2},\ldots,v_{k}>}_{u_{j}}|$
(3)

Without loss of generality, assume that corresponding to one of the minimum
values (if there are multiple minimum values) the POI Combination is
$<v^{{}^{\prime}}_{1},v^{{}^{\prime}}_{2},\ldots,v^{{}^{\prime}}_{k}>$. Then,
if we add $(d-D)$ (represented by $\epsilon$ in Algorithm 1) to $D$ then we
are going to have at least one feasible solution (and it will be optimal
solution as well). Also, as we are considering the minimum value among all the
entries of the $(k+2)$-th entries, so $\epsilon$ is the minimum distance to be
added to have at least one solution. It is easy to see that this step also
requires $\mathcal{O}(n^{k})$ time. Hence, the total time requirement of
Algorithm 1 is of $\mathcal{O}(n^{k}\cdot
b^{2}+n^{k}+n^{k})=\mathcal{O}(n^{k}\cdot b^{2})$. The extra space consumed by
Algorithm 1 is to store the $2$d-Array $Route\\_Matrix$ which will take
$\mathcal{O}((k+3)\cdot n^{k})=\mathcal{O}(n^{k})$ space and to store the
array $DIS$ which will take $\mathcal{O}(b)$ space. Hence, the total space
requirement is of $\mathcal{O}(n^{k}+b)$.

###### Theorem 3.1

Total time and space requirement of Algorithm 1 is of $\mathcal{O}(n^{k}\cdot
b^{2})$ and $\mathcal{O}(n^{k}+b)$, respectively.

Now, we state the correctness statement of Algorithm 1 in Theorem 3.2. Due to
the space constraint, we are unable to provide its proof.

###### Theorem 3.2

Given an Envy-free GTP Query Problem instance if its feasible solution set is
non-empty, then Algorithm 1 correctly returns an optimal solution.

Though this approach is easy to understand and simple to implement, the key
problem is excessive number of POI Combinations and hence it takes huge
computational time.

###  3.2 Solving the Problem Approximately

To solve our problem approximately, we use the Group Nearest Neighbor (GNN)
Search Approach. Given the source location of the group members
$\\{u^{s}_{1},u^{s}_{2},\ldots,u^{s}_{b}\\}$ among the POIs in the set
$\mathcal{V}_{1}$, one which is the GNN of the POIs
$\\{u^{s}_{1},u^{s}_{2},\ldots,u^{s}_{b}\\}$ is returned as the POI (say, in
this case $v^{{}^{\prime\prime}}_{1}$) of the first category. Next, among the
POIs in $\mathcal{V}_{2}$, one with the nearest neighbor of
$v^{{}^{\prime\prime}}_{1}$ is returned as the POI for the second category
(say, in this case $v^{{}^{\prime\prime}}_{2}$). Subsequently, among the POIs
in $\mathcal{V}_{3}$, one with the nearest neighbor of
$v^{{}^{\prime\prime}}_{2}$ is returned as the POI for the third category and
so on. This process is continued to search the POIs upto the $(k-1)$-th
category. Again, for the POI of $k$-th category again we use GNN Query
approach. So the POI for the $k$-th category $v^{{}^{\prime\prime}}_{k}$ is
basically among the POIs of $\mathcal{V}_{k}$ which is GNN of
$\\{u^{d}_{1},u^{d}_{2},\ldots,u^{d}_{b}\\}$. It is easy to observe that the
number of GNN Queries and NN Queries required are $2$ and $(k-2)$,
respectively. The entire process can be expressed as the following set of
equations.

$v^{{}^{\prime\prime}}_{1}=GNN(\\{u^{s}_{1},u^{s}_{2},\ldots,u^{s}_{b}\\},\mathcal{V}_{1})$

$v^{{}^{\prime\prime}}_{2}=NN(v^{{}^{\prime\prime}}_{1},\mathcal{V}_{2})$

$v^{{}^{\prime\prime}}_{3}=NN(v^{{}^{\prime\prime}}_{2},\mathcal{V}_{3})$

$\ldots$

$\ldots$

$v^{{}^{\prime\prime}}_{k-1}=NN(v^{{}^{\prime\prime}}_{k-2},\mathcal{V}_{k-1})$

$v^{{}^{\prime\prime}}_{k}=GNN(\\{u^{d}_{1},u^{d}_{2},\ldots,u^{d}_{b}\\},\mathcal{V}_{k})$

As we mentioned that between every POIs, the shortest path distance has been
computed and stored in a matrix. Now, to analyze this methodology we need to
find out how many table loop ups are required. It is easy to observe that for
finding $1$-st and $k$-th POI using GNN Queries the number of look ups will be
of $\mathcal{O}(b(n_{1}+n_{k}))$, and for finding the POIs for $2$-nd onwards
upto the $(k-1)$-th category is of
$\mathcal{O}(\underset{j\in\\{2,3,\ldots,k-1\\}}{\sum}n_{j})$. As the number
of POIs in the road network is $n$, hence the quantity
$\mathcal{O}(\underset{j\in\\{2,3,\ldots,k-1\\}}{\sum}n_{j})$ is upper bounded
by $\mathcal{O}(n)$. Hence, the total time requirement to solve the Envy-free
Trip Planning Problem approximately is of
$\mathcal{O}(b\cdot(n_{1}+n_{k})+n)$. It is easy to follow that this
methodology can be implemented in $\mathcal{O}(1)$ space. Hence, Theorem 3.3
holds.

###### Theorem 3.3

To solve the Envy-free Trip Planning Query Problem approximately using GNN
Query-based Approach, time and space requirement is of
$\mathcal{O}(b\cdot(n_{1}+n_{k})+n)$ and $\mathcal{O}(1)$, respectively.

##  4 Experimental Evaluation

In this section we describe the experimental evaluations of the proposed
methodologies. Initially, we start by describing the datasets.

#### Datasets

In this study, we use the following two datasets: Europe Road Network and
Minnesota Road Network. Both these datasets are publicly available and have
been downloaded from https://networkrepository.com/road.php [5]. The basic
statistics of the datasets have been mentioned in Table 1.

Table 1: Basic Statistics of the Datasets Dataset Name | $n$ | $m$ | Density | Avg. Degree
---|---|---|---|---
Europe Road Network | 1.2K | 1.4K | 0.00205794 | 2
Minnesota Road Network | 2.6K | 3.3K | 0.000946754 | 2

#### Goals of Our Experiments

In this study, we are interested to answer the following three research
questions:

  * •

With the increase of $D$, how the number of feasible solutions is changing?

  * •

In case of empty feasible solution set of the Envy-free GTP Query problem is
empty, how the $\epsilon$ value is changing with the change of $D$ value and
number of categories.

  * •

How much fast the approximate solution approach is compared to the exhaustive
search approach and what is the additional distance incurred for that.

#### Experiment Setup

In our study, we fix the group size to be $12$ and we consider the number of
categories as $2$, $3$, and $4$. For $D$, we start with the value $4400$ and
go upto $6000$ with an increment of $400$ every time.

#### Experimental Results with Discussions

In this section, we describe our experimental results. Figure 1 shows three
different plots for each dataset. The first one is $D$ Vs. No. of feasible
solutions. The second one is the $D$ Vs. $\epsilon$, and the third one is No.
of categories Vs. Computational Time. From the figure, we observe that for all
the datasets when the number of categories are fixed and the threshold
distance $D$ is increased then the number of feasible solutions are also
increased. As an example, for the Oldenburg Road Network dataset when the
number of categories are $2$ and threshold distance is increased from $4400$
to $6000$, the number of feasible solutions are increased from $3$ to $95$.
This is due to the following reason. As the threshold distance is increased,
then it is natural that more number of POI combinations will satisfy the Envy-
free ness Constraint. Also, for the same threshold distance, when the number
of categories are increased then also the number of feasible solutions are
increasing. As an example, for the Oldenburg Road Network when the threshold
distance is $4400$ and the number of categories are increased from $2$ to $4$,
then the number of feasible solutions are increased from $3$ to $6084$. When
the threshold distance is kept fixed and the number of categories are
increased then the number of POI Combinations are also increased. So, it is
quite natural that more number of POI Combinations will qualify for the Envy-
free ness constraint.

|
---|---
(a) Europe Road Network | (b) Minnesota Road Network
|
(c) Europe Road Network | (d) Minnesota Road Network
|
(e) Europe Road Network | (f) Minnesota Road Network
Figure 1: (Group Size, No. of POIs) Vs. Aggregated Distance Traveled for all
the methods on three datasets

Next, we describe our observations related to the second research question. It
is easy to observe that for a given problem instance of the Envy-free GTP
Query Problem is empty, if its feasible solution set is empty then there exist
one value for $d$ (as defined in Equation 3). Now, if the value of $D$ is
increased then as $\epsilon$ is defined as $(d-D)$, hence $\epsilon$ value
will be decreased. In our experiments also, we observe similar phenomenon. As
an example, for the Oldenburg Road Network Dataset when the number of
categories are fixed and the value of $D$ is increased from $2400$ to $2800$
then the $\epsilon$ value is decreased from $2159.48$ to $1759.67$. In other
words, the sum of $D$ and $\epsilon$ is constant. Similar observations are
made in case of the San Joqine Road Network dataset also. When the number of
categories are $3$ and the value of $D$ is increased from $3600$ to $4000$ the
value for $\epsilon$ has been decreased from $687.78$ to $282.56$. For the
Oldenburg Road Network dataset, when the $R$-Tree method is used and the
number of categories are increased from $2$ to $4$ for the fixed distance
threshold $2400$ then the $\epsilon$ value is decreased from $2369.49$ to
$2019.48$.

Finally, we describe our observation regarding the third research question. It
can be observed that the exhaustive search approach scans all possible POI
Combinations whether the envy-free ness constraint is satisfied or not.
However, in case of approximately solving this problem using Nearest Neighbor
Search-based approach we are exploiting the property of the $\mathcal{R}$-Tree
and we are not checking all the POI Combinations. Hence, it is expected that
the exhaustive search approach will take more time compared to the
$\mathcal{R}$-Tree based approach. In our experiments, we observe similar
phenomenon in particular for larger problem instance. Hence, from our
experiments we can conclude that if the spatial database is indexed with
$\mathcal{R}$-Trees then the problems can be solved in much faster way,
particularly when the problem instance size is larger.

##  5 Concluding Remarks and Future Research Directions

In this paper, we have introduced and studied two problems related to the trip
planning query namely Envey Free Group Trip Planning Query and Envy Free Group
Trip Planning Query with Minimum Additional Distance. With one example we have
shown that these problems are more relevant compared to their classical
version in practice. First, we relate the GTP Query Problem and Envey Free
Group Trip Planning Query Problem and show given an optimal solution for the
first problem how we can obtain the same for the second one. Subsequently, we
propose two efficient algorithms to solve both these problems with detailed
analysis and illustration with an example. Experimentation with real-world
datasets show that the proposed methodologies can recommend better quality
solution compared the baseline methods. In particular, the indexed based
approach leads to the best quality solution with reasonable computational
time. Now, our study can be extended in the following directions.

## References

  * [1] Diestel, R.: Graph theory 3rd ed. Graduate texts in mathematics 173 (2005)
  * [2] Hashem, T., Hashem, T., Ali, M.E., Kulik, L.: Group trip planning queries in spatial databases. In: International Symposium on Spatial and Temporal Databases. pp. 259–276. Springer (2013)
  * [3] Huang, H., Gartner, G., Krisp, J.M., Raubal, M., Van de Weghe, N.: Location based services: ongoing evolution and research agenda. Journal of Location Based Services 12(2), 63–93 (2018)
  * [4] Li, W., Zhu, H., Liu, W., Yin, J., Xu, J.: Optimal sequenced route query with poi preferences. In: International Conference on Database Systems for Advanced Applications. pp. 457–473. Springer (2021)
  * [5] Rossi, R., Ahmed, N.: The network data repository with interactive graph analytics and visualization. In: Twenty-ninth AAAI conference on artificial intelligence (2015)
  * [6] Singhal, M., Banerjee, S.: Group trip planning queries on road networks using geo-tagged textual information. In: International Conference on Advanced Data Mining and Applications. pp. 243–257. Springer (2022)
  * [7] Tabassum, A., Barua, S., Hashem, T., Chowdhury, T.: Dynamic group trip planning queries in spatial databases. In: Proceedings of the 29th international conference on scientific and statistical database management. pp. 1–6 (2017)
  * [8] Zheng, Y.: Trajectory data mining: an overview. ACM Transactions on Intelligent Systems and Technology (TIST) 6(3), 1–41 (2015)