# Inequalities for the overpartition function arising from determinants

Gargi Mukherjee  Institute for Algebra, Johannes Kepler University,
Altenberger Straße 69, A-4040 Linz, Austria.  gargi.mukherjee@dk-
compmath.jku.at

###### Abstract.

Let $\overline{p}(n)$ denote the overpartition funtion. This paper presents
the $2$-$\log$-concavity property of $\overline{p}(n)$ by considering a more
general inequality of the following form

$\begin{vmatrix}\overline{p}(n)&\overline{p}(n+1)&\overline{p}(n+2)\\\
\overline{p}(n-1)&\overline{p}(n)&\overline{p}(n+1)\\\
\overline{p}(n-2)&\overline{p}(n-1)&\overline{p}(n)\end{vmatrix}>0,$

which holds for all $n\geq 42$.

Mathematics Subject Classifications. Primary 05A20; 11C20.
Keywords. Overpartition; determinant; 2-log-concavity.

##  1\. Introduction

A sequence $a_{0},a_{1},\dots,a_{n}$ of real numbers is said to be unimodal if
for some $0\leq k\leq n$ we have $a_{0}\leq a_{1}\leq\dots\leq a_{k-1}\leq
a_{k}\geq a_{k+1}\geq\dots\geq a_{n}$ and $\log$-concave if $a^{2}_{j}\geq
a_{j-1}a_{j+1}$ for all $1\leq j\leq n-1$. Note that a $\log$-concave sequence
of positive numbers is unimodal. We say that the sequence
$a_{0},a_{1},\dots,a_{n}$ has no internal zeros if there do not exist integers
$0\leq i<j<k\leq n$ satisfying $a_{i}\neq 0$, $a_{j}=0$ and $a_{k}\neq 0$.
Then a nonnegative $\log$-concave sequence with no internal zeroes is
unimodal. The study of $\log$-concavity problems of sequences shares an
intimate connection with zeros of polynomials, due to Newton through the
following result; see for example [11, p. 52].

######  Theorem 1.1.

Let

$P(x)=\sum_{j=0}^{n}\binom{n}{j}a_{j}x^{j}$

be a (real) polynomial with real zeros. Then $a^{2}_{j}\geq a_{j-1}a_{j+1}$
for all $j$.

We call a real polynomial $P(x)=\sum_{j=0}^{n}a_{j}x^{j}$ $\log$-concave if
its coefficient sequence is $\log$-concave. More generally, the following
theorem provides a necessary condition for a real polynomial to have only real
zeros.

######  Theorem 1.2.

$($[2, Theorem 1.2.1]$)$ Let $P(x)=\sum_{i=0}^{d}a_{j}x^{j}$ be a polynomial
with nonnegative coefficients with only real zeros. Then the sequence
$\\{a_{k}\\}_{0\leq k\leq d}$ is $\log$-concave with no internal zeros; in
particular, it is unimodal.

Before we state a more general version of Theorem 1.2, let us introduce the
theory of total positivity in brevity. A matrix $A$ with entries in real
number is called totally positive if the determinant of each of its minors is
nonnegative. For a sequence $(a_{k})_{k\geq 0}$, define its Toeplitz matrix
$T$ by

$T=\begin{pmatrix}a_{0}&&&&\\\ a_{1}&a_{0}&&&\\\ a_{2}&a_{1}&a_{0}&&\\\
a_{3}&a_{2}&a_{1}&a_{0}\\\ \vdots&\vdots&\vdots&\vdots&\ddots\end{pmatrix}.$
(1.1)

We say that $(a_{k})_{k\geq 0}$ is a totally positive sequence or Pólya
frequency sequence if its associated Toeplitz matrix is totally positive. The
following theorem shares a deep entanglement between combinatorics of zeros of
polynomials and total positivity of its coefficient sequence.

######  Theorem 1.3 (Aissen et al. [1]).

Let $P(x)=\sum_{i=0}^{d}a_{j}x^{j}$ be a real polynomial with nonnegative
coefficients. Then $P(x)$ has only real zeros if and only if its coefficient
sequence is a Pólya frequency sequence.

Associated to a real polynomial $P(x)=\sum_{i=0}^{d}a_{j}x^{j}$ with only real
zeros, a sequence $(\gamma_{k})_{k\geq 0}$ is called multiplier sequence if
the corresponding polynomial
$P_{\gamma}(x)=\sum_{i=0}^{d}\gamma_{j}a_{j}x^{j}$ also has only real zeros.
This theory laid its foundation in the seminal work of Pólya and Schur [18].
In this context, Craven and Csordas obtained the following theorem that
further reduces to $2$-$\log$-concavity $(\text{or double Tur\'{a}n
inequality})$ of a multiplier sequence.

######  Theorem 1.4 (Theorem 2.13, [7]).

If $(\gamma_{k})_{k\geq 0}$ with $\gamma_{k}>0$, is a multiplier sequence,
then

$\begin{vmatrix}\gamma_{k}&\gamma_{k+1}&\gamma_{k+2}\\\
\gamma_{k-1}&\gamma_{k}&\gamma_{k+1}\\\
\gamma_{k-2}&\gamma_{k-1}&\gamma_{k}\\\ \end{vmatrix}\geq 0,\ \ \text{for}\
k\in\mathbb{Z}_{\geq 2}.$ (1.2)

Since $\gamma_{k}>0$, (1.2) is equivalent to the following inequality

$(\gamma^{2}_{k}-\gamma_{k-1}\gamma_{k+1})^{2}-(\gamma^{2}_{k-1}-\gamma_{k-2}\gamma_{k})(\gamma^{2}_{k+1}-\gamma_{k}\gamma_{k+2})\geq
0,$

which amounts to say that $(\gamma_{k})_{k\geq 0}$ is $2$-$\log$-concave.
Now we turn to discuss briefly about partitions through the lens of
combinatorial analysis, as discussed in the previous paragraph. A partition of
a positive integer $n$ is a nonincreasing sequence
$(\lambda_{1},\lambda_{2},\dots,\lambda_{r})$ of positive integers with
$\lambda_{1}+\lambda_{2}+\dots+\lambda_{r}=n$ and $p(n)$ denotes the number of
partitions of $n$. The systematic study of partitions dates back to Euler.
Rigorous analytic approach comes into play in the theory of partitions since
the foundational work of Hardy and Ramanujan [12]. Hardy and Ramanujan
employed the celebrated circle method in order to explicitly describe the
asymptotics of $p(n)$, specifically, given by

$p(n)\sim\frac{1}{4\sqrt{3}n}e^{\pi\sqrt{\frac{2n}{3}}}\ \
(n\rightarrow\infty).$ (1.3)

Later Rademacher [19] refined the formulation of Hardy and Ramanujan to set a
convergent series expression for $p(n)$ and an error bound was given due to
Lehmer [15]. Log-concavity of $p(n)$ has been studied independently by Nicolas
[17] and by DeSalvo and Pak [8] by confirming a conjecture of Chen [3]. Since
then the study on inequalities of the partition function from combinatorial
analysis perspective has been documented in the works of Chen et al. [4], [5].
Recently Griffin, Ono, Rolen and Zagier [10] consider a more general paradigm
to trace the zeros of a certain polynomial, called Jensen polynomials
associated with $p(n)$, defined by

$J^{d,n}_{p}(x)=\sum_{j=0}^{d}\binom{d}{j}p(n+j)x^{j}.$

For a more detail study on hyperbolicity of $J^{d,n}_{p}(x)$, see [14].
Corteel and Lovejoy [6] initiated a broad generalization of partitions, called
overpartition that offers a panorama of combinatorial perspective of basic
hypergeometric series. An overpartition of $n$ is a nonincreasing sequence of
natural numbers whose sum is $n$ in which the first occurrence (equivalently,
the final occurrence) of a number may be overlined and $\overline{p}(n)$
denotes the number of overpartitions of $n$. For convenience, define
$\overline{p}(0)=1$. For example, there are $8$ overpartitions of $3$
enumerated by
$3,\overline{3},2+1,\overline{2}+1,2+\overline{1},\overline{2}+\overline{1},1+1+1,\overline{1}+1+1$.
Similar to the Hardy-Ramanujan-Rademacher type formula for $\overline{p}(n)$,
Zuckerman [21] showed that

$\overline{p}(n)=\frac{1}{2\pi}\underset{2\nmid
k}{\sum_{k=1}^{\infty}}\sqrt{k}\underset{(h,k)=1}{\sum_{h=0}^{k-1}}\dfrac{\omega(h,k)^{2}}{\omega(2h,k)}e^{-\frac{2\pi
inh}{k}}\dfrac{d}{dn}\biggl{(}\dfrac{\sinh\frac{\pi\sqrt{n}}{k}}{\sqrt{n}}\biggr{)},$
(1.4)

where

$\omega(h,k)=\text{exp}\Biggl{(}\pi
i\sum_{r=1}^{k-1}\dfrac{r}{k}\biggl{(}\dfrac{hr}{k}-\biggl{\lfloor}\dfrac{hr}{k}\biggr{\rfloor}-\dfrac{1}{2}\biggr{)}\Biggr{)}$

for positive integers $h$ and $k$. In order to prove $\log$-concavity of
$\overline{p}(n)$, Engel [9] provided an error term for $\overline{p}(n)$

$\overline{p}(n)=\frac{1}{2\pi}\underset{2\nmid
k}{\sum_{k=1}^{N}}\sqrt{k}\underset{(h,k)=1}{\sum_{h=0}^{k-1}}\dfrac{\omega(h,k)^{2}}{\omega(2h,k)}e^{-\frac{2\pi
inh}{k}}\dfrac{d}{dn}\biggl{(}\dfrac{\sinh\frac{\pi\sqrt{n}}{k}}{\sqrt{n}}\biggr{)}+R_{2}(n,N),$
(1.5)

where

$\bigl{|}R_{2}(n,N)\bigr{|}<\dfrac{N^{5/2}}{\pi
n^{3/2}}\sinh\biggl{(}\dfrac{\pi\sqrt{n}}{N}\biggr{)}.$ (1.6)

Along the lines of works of Chen et al. in context of the partition function,
somewhat similar research works on inequalities for $\overline{p}(n)$ has
already been recorded in [20] and [16].
Jia and Wang [13] examined determinantal inequalities for $p(n)$ arising from
the theory of total positivity, by set up the following theorem. As a
corollary of Theorem 1.5, they proved $2$ $\log$-concavity of $p(n)$.

######  Theorem 1.5 (Theorem 1.5, [13]).

Let $p(n)$ denote the partition function and

$M_{3}(p(n))=\begin{pmatrix}p(n)&p(n+1)&p(n+2)\\\ p(n-1)&p(n)&p(n+1)\\\
p(n-2)&p(n-1)&p(n)\end{pmatrix}$ (1.7)

Then for all $n\geq 222$, we have

$\displaystyle\det M_{3}(p(n))$ $\displaystyle>$ $\displaystyle 0.$ (1.8)

In this paper, our primary goal is to prove $2$-$\log$-concavity of
$\overline{p}(n)$. In order to prove this, we set up a similar device as that
of Theorem 1.5 but in context of overpartitions. In particular, we shall prove
the following result.

######  Theorem 1.6.

Let $\overline{p}(n)$ denote the overpartition funtion. Then for all $n\geq
42$, we have

$\begin{vmatrix}\overline{p}(n)&\overline{p}(n+1)&\overline{p}(n+2)\\\
\overline{p}(n-1)&\overline{p}(n)&\overline{p}(n+1)\\\
\overline{p}(n-2)&\overline{p}(n-1)&\overline{p}(n)\end{vmatrix}>0.$ (1.9)

Theorem 1.6 straight away implies the $2$-$\log$-concavity of
$\overline{p}(n)$, precisely

######  Theorem 1.7.

For all $n\geq 42$,

$(\overline{p}(n)^{2}-\overline{p}(n-1)\overline{p}(n+1))^{2}-(\overline{p}(n-1)^{2}-\overline{p}(n-2)\overline{p}(n))\
(\overline{p}(n+1)^{2}-\overline{p}(n)\overline{p}(n+2))>0.$ (1.10)

We organize this paper in the following format. In Section 2, we set up the
premise first by introducing an inequality for $\overline{p}(n)$ and
reformulate Theorem 1.6, given in Theorem 2.5, followed by the documentation
of Theorem 2.6 and 2.7. This foundation enables us to provide the proof of
Theorem 2.5 in Section 3. In the end, we discuss how one can guess an infinite
number of inequalities for the overpartition function by considering totally
positive matrix of order $k\times k$ with $k\in\mathbb{Z}_{\geq 2}$, described
in Problem 4.1.

##  2\. Inequality for $\overline{p}(n)$ and its consequences

The principal aim of this section is to construct the machinery in order to
prove Theorem 1.6, the primary objective of this paper. To some extent, we
follow a similar line of argument as described in the work of Jia and Wang
[13]. We will see that the Theorem 2.6 and 2.7 are the key tools to prove
Theorem 2.5, a reprise version of Theorem 1.6. First we need to estimate the
quotient
$u_{n}=\frac{\overline{p}(n-1)\overline{p}(n+1)}{\overline{p}(n)^{2}}$ by
showing its upper and lower bound $g(n)$ and $f(n)$ respectively (cf. Theorem
2.3), derived from the inequality $B_{1}(n)<\overline{p}(n)<B_{2}(n)$ with
$B_{1}(n),B_{2}(n)$ given in (2) and as an immediate consequence, we get the
inequality (2.15) for $s(n)=u_{n-1}+u_{n+1}-u_{n-1}u_{n+1}$ as follows; for
all $n\geq 3$,

$s_{1}(n)<s(n)<s_{2}(n),$

where $s_{1}(n)$ and $s_{2}(n)$ are combinations of $f(n+1),f(n-1),g(n+1)$,
and $g(n-1)$. Therefore, to prove Theorem 2.6 and 2.7, the principal idea
behind it is to approximate $s_{2}(n)$ (cf. (2.14)) and
$s_{1}(n)g(n)^{2}-2g(n)+1$ (cf. (2.45)) by rational functions in $y\
(=y(n)=\pi\sqrt{n})$ (cf. (2.40) and (2.58)). In order to arrive at such
estimation to ease the computation, it is necessary to bound the error term
$\widetilde{T}(n)$ by $\mu(n)^{-2m}$ for some $m\in\mathbb{Z}_{\geq 1}$
because our estimation turns out to get a suitable polynomial approximation of
$r,\ x,\ z$ and $w$ (cf. (2.18)) in terms of $y$. In our case, it is
sufficient to consider $m=3$, as stated in Lemma 2.1.
We denote $\mu(n)=\pi\sqrt{n}$ and define

$\displaystyle B_{1}(n)$ $\displaystyle=$
$\displaystyle\frac{e^{\mu(n)}}{8n}\Bigl{(}1-\frac{1}{\mu(n)}-\frac{1}{\mu(n)^{6}}\Bigr{)}$
$\displaystyle\text{and}\ \ B_{2}(n)$ $\displaystyle=$
$\displaystyle\frac{e^{\mu(n)}}{8n}\Bigl{(}1-\frac{1}{\mu(n)}+\frac{1}{\mu(n)^{6}}\Bigr{)}.$
(2.1)

######  Lemma 2.1.

For all $n\geq 94$, we have

$B_{1}(n)<\overline{p}(n)<B_{2}(n).$ (2.2)

  * Proof.

From [16, eqn. (3.5)], it follows that

$\overline{p}(n)=\frac{e^{\mu(n)}}{8n}\Bigl{(}1-\frac{1}{\mu(n)}+\widetilde{T}(n)\Bigr{)}$
(2.3)

where

$\widetilde{T}(n)=\Bigl{(}1+\dfrac{1}{\mu(n)}\Bigr{)}e^{-2\mu(n)}+\dfrac{8n}{e^{\mu(n)}}R_{2}(n,2),$

and $R_{2}(n,2)$ is the error term of (1.4), given in (1.6). By [16, eqn.
(3.6)], we have

$|\widetilde{T}(n)|<10\ e^{-\frac{1}{2}\mu(n)}.$ (2.4)

Now,

$10\ e^{-\frac{1}{2}\mu(n)}<\frac{1}{\mu(n)^{6}}\ \ \ \text{for all}\ n\geq
275.$ (2.5)

(2.3)-(2.5) altogether imply (2.2) for all $n\geq 275$. We finish the proof by
confirming (2.2) by checking numerically in Mathematica for all $94\leq n\leq
274$. ∎

######  Remark 2.2.

We note that the upper bound of the absolute value of error term
$\widetilde{T}(n)$ can be improved by considering a more generalized version
of (2.5) of the following form: there exists $N(m)\in\mathbb{Z}_{\geq 1}$,
such that for all $n\geq N(m)$,

$10\ e^{-\frac{1}{2}\mu(n)}<\frac{1}{\mu(n)^{m}}.$ (2.6)

Let

$u_{n}=\frac{\overline{p}(n-1)\overline{p}(n+1)}{\overline{p}(n)^{2}}$ (2.7)

and consequently, denote

$s(n)=u_{n-1}+u_{n+1}-u_{n-1}u_{n+1}.$ (2.8)

Following the notations as given in [13], we set

$r=\mu(n-2),\ \ x=\mu(n-1),\ \ y=\mu(n),\ \ z=\mu(n+1),\ \ w=\mu(n+2),$ (2.9)

and

$\displaystyle f(n)$ $\displaystyle=$ $\displaystyle
e^{x-2y+z}\frac{(x^{6}-x^{5}-1)y^{16}(z^{6}-z^{5}-1)}{x^{8}(y^{6}-y^{5}+1)^{2}z^{8}},$
(2.10) $\displaystyle g(n)$ $\displaystyle=$ $\displaystyle
e^{x-2y+z}\frac{(x^{6}-x^{5}+1)y^{16}(z^{6}-z^{5}+1)}{x^{8}(y^{6}-y^{5}-1)^{2}z^{8}}.$
(2.11)

As an immediate consequence of Lemma 2.1, we have the following theorem.

######  Theorem 2.3.

For all $n\geq 94$,

$f(n)<u_{n}<g(n).$ (2.12)

We begin with the following setup. Define

$s_{1}(n)=f(n-1)+f(n+1)-g(n-1)g(n+1),$ (2.13)

and

$s_{2}(n)=g(n-1)+g(n+1)-f(n-1)f(n+1).$ (2.14)

As a corollary of Theorem 2.3, we arrive at the following inequality for
$s(n)$.

######  Corollary 2.4.

For $n\geq 91$, we have

$s_{1}(n)<s(n)<s_{2}(n).$ (2.15)

Now we interpret the Theorem 1.6 in terms of a polynomial expression in $s(n)$
and $u_{n}$ (cf. (2.7) and (2.8)), given as follows

######  Theorem 2.5.

For all $n\geq 42$, we have

$s(n)u^{2}_{n}-2u_{n}+1>0.$ (2.16)

To prove (2.16), first it is required to estimate upper bound of $s(n)$, given
by studying $s_{2}(n)$, as follows;

######  Theorem 2.6.

For all $n\geq 3$, we have

$s_{2}(n)<1.$ (2.17)

  * Proof.

For $n\geq 3$, rewriting (2.9), we have

$r=\sqrt{y^{2}-2\pi^{2}},\ \ x=\sqrt{y^{2}-\pi^{2}},\ \
z=\sqrt{y^{2}+\pi^{2}},\ \ w=\sqrt{y^{2}+2\pi^{2}}.$ (2.18)

Expanding $r,\ x,\ z\ \text{and}\ w$ in terms of $y$, it follows that

$\displaystyle r$ $\displaystyle=$ $\displaystyle
y-\frac{\pi^{2}}{y}-\frac{\pi^{4}}{y^{3}}-\frac{\pi^{6}}{2y^{5}}-\frac{5\pi^{8}}{8y^{7}}-\frac{7\pi^{10}}{8y^{9}}-\frac{21\pi^{12}}{16y^{11}}-\frac{33\pi^{14}}{16y^{13}}+O\Bigl{(}\frac{1}{y^{15}}\Bigr{)},$
$\displaystyle x$ $\displaystyle=$ $\displaystyle
y-\frac{\pi^{2}}{2y}-\frac{\pi^{4}}{8y^{3}}-\frac{\pi^{6}}{16y^{5}}-\frac{5\pi^{8}}{128y^{7}}-\frac{7\pi^{10}}{256y^{9}}-\frac{21\pi^{12}}{1024y^{11}}-\frac{33\pi^{14}}{2048y^{13}}+O\Bigl{(}\frac{1}{y^{15}}\Bigr{)},$
$\displaystyle z$ $\displaystyle=$ $\displaystyle
y+\frac{\pi^{2}}{2y}-\frac{\pi^{4}}{8y^{3}}+\frac{\pi^{6}}{16y^{5}}-\frac{5\pi^{8}}{128y^{7}}+\frac{7\pi^{10}}{256y^{9}}-\frac{21\pi^{12}}{1024y^{11}}+\frac{33\pi^{14}}{2048y^{13}}+O\Bigl{(}\frac{1}{y^{15}}\Bigr{)},$
$\displaystyle w$ $\displaystyle=$ $\displaystyle
y+\frac{\pi^{2}}{y}-\frac{\pi^{4}}{y^{3}}+\frac{\pi^{6}}{2y^{5}}-\frac{5\pi^{8}}{8y^{7}}+\frac{7\pi^{10}}{8y^{9}}-\frac{21\pi^{12}}{16y^{11}}+\frac{33\pi^{14}}{16y^{13}}+O\Bigl{(}\frac{1}{y^{15}}\Bigr{)}.$

It can be easily verified that for all $n\geq 59$,

$r_{1}<r<r_{2},$ (2.19) $x_{1}<x<x_{2},$ (2.20) $z_{1}<z<z_{2},$ (2.21)

and

$w_{1}<w<w_{2},$ (2.22)

with

$\displaystyle r_{1}$ $\displaystyle=$ $\displaystyle
y-\frac{\pi^{2}}{y}-\frac{\pi^{4}}{y^{3}}-\frac{\pi^{6}}{2y^{5}}-\frac{5\pi^{8}}{8y^{7}}-\frac{7\pi^{10}}{8y^{9}}-\frac{21\pi^{12}}{16y^{11}}-\frac{34\pi^{14}}{16y^{13}},$
$\displaystyle r_{2}$ $\displaystyle=$ $\displaystyle
y-\frac{\pi^{2}}{y}-\frac{\pi^{4}}{y^{3}}-\frac{\pi^{6}}{2y^{5}}-\frac{5\pi^{8}}{8y^{7}}-\frac{7\pi^{10}}{8y^{9}}-\frac{21\pi^{12}}{16y^{11}}-\frac{33\pi^{14}}{16y^{13}},$
$\displaystyle x_{1}$ $\displaystyle=$ $\displaystyle
y-\frac{\pi^{2}}{2y}-\frac{\pi^{4}}{8y^{3}}-\frac{\pi^{6}}{16y^{5}}-\frac{5\pi^{8}}{128y^{7}}-\frac{7\pi^{10}}{256y^{9}}-\frac{21\pi^{12}}{1024y^{11}}-\frac{34\pi^{14}}{2048y^{13}},$
$\displaystyle x_{2}$ $\displaystyle=$ $\displaystyle
y-\frac{\pi^{2}}{2y}-\frac{\pi^{4}}{8y^{3}}-\frac{\pi^{6}}{16y^{5}}-\frac{5\pi^{8}}{128y^{7}}-\frac{7\pi^{10}}{256y^{9}}-\frac{21\pi^{12}}{1024y^{11}}-\frac{33\pi^{14}}{2048y^{13}},$
$\displaystyle z_{1}$ $\displaystyle=$ $\displaystyle
y+\frac{\pi^{2}}{2y}-\frac{\pi^{4}}{8y^{3}}+\frac{\pi^{6}}{16y^{5}}-\frac{5\pi^{8}}{128y^{7}}+\frac{7\pi^{10}}{256y^{9}}-\frac{21\pi^{12}}{1024y^{11}},$
$\displaystyle z_{2}$ $\displaystyle=$ $\displaystyle
y+\frac{\pi^{2}}{2y}-\frac{\pi^{4}}{8y^{3}}+\frac{\pi^{6}}{16y^{5}}-\frac{5\pi^{8}}{128y^{7}}+\frac{7\pi^{10}}{256y^{9}}-\frac{21\pi^{12}}{1024y^{11}}+\frac{33\pi^{14}}{2048y^{13}},$
$\displaystyle w_{1}$ $\displaystyle=$ $\displaystyle
y+\frac{\pi^{2}}{y}-\frac{\pi^{4}}{y^{3}}+\frac{\pi^{6}}{2y^{5}}-\frac{5\pi^{8}}{8y^{7}}+\frac{7\pi^{10}}{8y^{9}}-\frac{21\pi^{12}}{16y^{11}},$
$\displaystyle w_{2}$ $\displaystyle=$ $\displaystyle
y+\frac{\pi^{2}}{y}-\frac{\pi^{4}}{y^{3}}+\frac{\pi^{6}}{2y^{5}}-\frac{5\pi^{8}}{8y^{7}}+\frac{7\pi^{10}}{8y^{9}}-\frac{21\pi^{12}}{16y^{11}}+\frac{33\pi^{14}}{16y^{13}}.$

Following the definition of $s_{2}(n)$ given in (2.14), we see that it
suffices to estimate $f(n-1)$, $f(n+1)$, $g(n-1)$ and $g(n+1)$. Now, we
observe that each of these four functions consists of two factors, the
exponential factor and the rational function in variables $x,y$ and $z$ (cf.
(2.10) and (2.11)). This suggests that it is enough to estimate $e^{r-2x+y}$,
$e^{y-2z+w}$, $h(n-1)$, $h(n+1)$, $q(n-1)$ and $q(n+1)$ individually, where

$f(n)=e^{x-2y+z}h(n),\ \ g(n)=e^{x-2y+z}q(n),$ (2.23)

with

$h(n)=\frac{(x^{6}-x^{5}-1)y^{16}(z^{6}-z^{5}-1)}{x^{8}(y^{6}-y^{5}+1)^{2}z^{8}},$
(2.24)

and

$q(n)=\frac{(x^{6}-x^{5}+1)y^{16}(z^{6}-z^{5}+1)}{x^{8}(y^{6}-y^{5}-1)^{2}z^{8}}.$
(2.25)

First, let us consider the exponential factors $e^{r-2x+y}$ and $e^{y-2z+w}$.
By (2.19)-(2.22), for all $n\geq 59$, it follows that

$e^{r_{1}-2x_{2}+y}<e^{r-2x+y}<e^{r_{2}-2x_{1}+y},$ (2.26)
$e^{y-2z_{2}+w_{1}}<e^{y-2z+w}<e^{y-2z_{1}+w_{2}}.$ (2.27)

Next, we estimate (2.26) and (2.27) by Taylor expansion of the exponential
function in order to get bounds in terms of rational function in $y$. For
convenience, set

$\Phi(t)=\sum_{i=0}^{6}\frac{t^{i}}{i!},$ (2.28)

and

$\phi(t)=\sum_{i=0}^{7}\frac{t^{i}}{i!}.$ (2.29)

For $t\in\mathbb{R}_{<0}$,

$\phi(t)<e^{t}<\Phi(t).$ (2.30)

We note that

$r_{2}-2x_{1}+y=-\frac{\pi^{4}\Bigl{(}1039\pi^{10}+651\pi^{8}y^{2}+420\pi^{6}y^{4}+280\pi^{4}y^{6}+192\pi^{2}y^{8}+128y^{10}\Bigr{)}}{512y^{13}}<0$

and for all $n\geq 2$,

$y-2z_{1}+w_{2}=\frac{\pi^{4}\Bigl{(}1056\pi^{10}-651\pi^{8}y^{2}+420\pi^{6}y^{4}-280\pi^{4}y^{6}+192\pi^{2}y^{8}-128y^{10}\Bigr{)}}{512y^{13}}<0.$

Putting (2.30) into (2.26) and (2.27), we get for $n\geq 59$,

$\phi(r_{1}-2x_{2}+y)<e^{r-2x+y}<\Phi(r_{2}-2x_{1}+y),$ (2.31)

and

$\phi(y-2z_{2}+w_{1})<e^{y-2z+w}<\Phi(y-2z_{1}+w_{2}).$ (2.32)

Finally, it remains to estimate $h(n-1)$, $h(n+1)$, $q(n-1)$ and $q(n+1)$. We
rewrite these four functions as

$\begin{split}h(n-1)&=\frac{x^{16}\beta(r)\beta(y)}{r^{8}y^{8}\alpha(x)^{2}},\
\ \ h(n+1)=\frac{z^{16}\beta(y)\beta(w)}{w^{8}y^{8}\alpha(z)^{2}},\\\
q(n-1)&=\frac{x^{16}\alpha(r)\alpha(y)}{r^{8}y^{8}\beta(x)^{2}},\ \ \
q(n+1)=\frac{z^{16}\alpha(y)\alpha(w)}{w^{8}y^{8}\beta(z)^{2}},\end{split}$

where

$\alpha(t)=t^{6}-t^{5}+1\ \ \text{and}\ \ \beta(t)=t^{6}-t^{5}-1.$ (2.33)

Using (2.19)-(2.22), for $n\geq 59$, we put down a list of inequalities as
follows

$\displaystyle r^{6}-r_{2}r^{4}+1$ $\displaystyle<\alpha(r)$
$\displaystyle<r^{6}-r_{1}r^{4}+1,$ $\displaystyle x^{6}-x_{2}x^{4}+1$
$\displaystyle<\alpha(x)$ $\displaystyle<x^{6}-x_{1}x^{4}+1,$ $\displaystyle
z^{6}-z_{2}z^{4}+1$ $\displaystyle<\alpha(z)$
$\displaystyle<z^{6}-z_{1}z^{4}+1,$ $\displaystyle w^{6}-w_{2}w^{4}+1$
$\displaystyle<\alpha(w)$ $\displaystyle<w^{6}-w_{1}w^{4}+1,$ $\displaystyle
r^{6}-r_{2}r^{4}-1$ $\displaystyle<\beta(r)$
$\displaystyle<r^{6}-r_{1}r^{4}-1,$ $\displaystyle x^{6}-x_{2}x^{4}-1$
$\displaystyle<\beta(x)$ $\displaystyle<x^{6}-x_{1}x^{4}-1,$ $\displaystyle
z^{6}-z_{2}z^{4}-1$ $\displaystyle<\beta(z)$
$\displaystyle<z^{6}-z_{1}z^{4}-1,$ $\displaystyle w^{6}-w_{2}w^{4}-1$
$\displaystyle<\beta(w)$ $\displaystyle<w^{6}-w_{1}w^{4}-1,$ $\displaystyle
x^{12}-2x_{2}x^{10}+x^{10}+2x^{6}-2x_{2}x^{4}+1$ $\displaystyle<\alpha(x)^{2}$
$\displaystyle<x^{12}-2x_{1}x^{10}+x^{10}+2x^{6}-2x_{1}x^{4}+1,$
$\displaystyle z^{12}-2z_{2}z^{10}+z^{10}+2z^{6}-2z_{2}z^{4}+1$
$\displaystyle<\alpha(z)^{2}$
$\displaystyle<z^{12}-2z_{1}z^{10}+z^{10}+2z^{6}-2z_{1}z^{4}+1,$
$\displaystyle x^{12}-2x_{2}x^{10}+x^{10}-2x^{6}+2x_{1}x^{4}+1$
$\displaystyle<\beta(x)^{2}$
$\displaystyle<x^{12}-2x_{1}x^{10}+x^{10}-2x^{6}+2x_{2}x^{4}+1,$
$\displaystyle z^{12}-2z_{2}z^{10}+z^{10}-2z^{6}+2z_{1}z^{4}+1$
$\displaystyle<\beta(z)^{2}$
$\displaystyle<z^{12}-2z_{1}z^{10}+z^{10}-2z^{6}+2z_{2}z^{4}+1.$

By application of the above inequalities, it follows that

$h(n-1)>\frac{(r^{6}-r_{2}r^{4}-1)x^{16}(y^{6}-y^{5}-1)}{r^{8}y^{8}(x^{12}-2x_{1}x^{10}+x^{10}+2x^{6}-2x_{1}x^{4}+1)},$
(2.35)
$h(n+1)>\frac{(w^{6}-w_{2}w^{4}-1)z^{16}(y^{6}-y^{5}-1)}{w^{8}y^{8}(z^{12}-2z_{1}z^{10}+z^{10}+2z^{6}-2z_{1}z^{4}+1)},$
(2.36)
$q(n-1)<\frac{(r^{6}-r_{1}r^{4}+1)x^{16}(y^{6}-y^{5}+1)}{r^{8}y^{8}(x^{12}-2x_{2}x^{10}+x^{10}-2x^{6}+2x_{1}x^{4}+1)},$
(2.37)
$q(n+1)<\frac{(w^{6}-w_{1}r^{4}+1)z^{16}(y^{6}-y^{5}+1)}{r^{8}y^{8}(z^{12}-2z_{2}z^{10}+z^{10}-2z^{6}+2z_{1}z^{4}+1)}.$
(2.38)

Invoking (2.31)-(2.32) and (2.35)-(2.38) into (2.23), for $n\geq 59$, we have

$\displaystyle g(n-1)$ $\displaystyle<$ $\displaystyle
R_{1}(y)=\frac{(r^{6}-r_{1}r^{4}+1)x^{16}(y^{6}-y^{5}+1)\Phi(r_{2}-2x_{1}+y)}{r^{8}y^{8}(x^{12}-2x_{2}x^{10}+x^{10}-2x^{6}+2x_{1}x^{4}+1)},$
$\displaystyle g(n+1)$ $\displaystyle<$ $\displaystyle
R_{2}(y)=\frac{(w^{6}-w_{1}r^{4}+1)z^{16}(y^{6}-y^{5}+1)\Phi(y-2z_{1}+w_{2})}{w^{8}y^{8}(z^{12}-2z_{2}z^{10}+z^{10}-2z^{6}+2z_{1}z^{4}+1)},$
$\displaystyle f(n-1)$ $\displaystyle>$ $\displaystyle
R_{3}(y)=\frac{(r^{6}-r_{2}r^{4}-1)x^{16}(y^{6}-y^{5}-1)\phi(r_{1}-2x_{2}+y)}{r^{8}y^{8}(x^{12}-2x_{1}x^{10}+x^{10}+2x^{6}-2x_{1}x^{4}+1)},$
$\displaystyle f(n+1)$ $\displaystyle>$ $\displaystyle
R_{4}(y)=\frac{(w^{6}-w_{2}w^{4}-1)z^{16}(y^{6}-y^{5}-1)\phi(y-2z_{2}+w_{1})}{w^{8}y^{8}(z^{12}-2z_{1}z^{10}+z^{10}+2z^{6}-2z_{1}z^{4}+1)}.$

By definition of $s_{2}(n)$ (cf. (2.14)), it suffices to prove that

$R_{1}(y)+R_{2}(y)-R_{3}(y)R_{4}(y)-1<0$ (2.39)

We can reduce $R_{1}(y)+R_{2}(y)-R_{3}(y)R_{4}(y)-1$ into a rational function
in $y$; i.e,

$R_{1}(y)+R_{2}(y)-R_{3}(y)R_{4}(y)-1=\frac{N_{1}(y)}{D_{1}(y)},$ (2.40)

where $N_{1}(y)$ and $D_{1}(y)$ are polynomials in $y$ with respective degree
$324$ and $330$. In order to prove (2.39), it is equivalent to prove
$N_{1}(y)D_{1}(y)<0$. We write

$N_{1}(y)D_{1}(y)=\sum_{i=0}^{654}a_{i}y^{i},$ (2.41)

where

$a_{654}=2^{348}\cdot 3^{8}\cdot 5^{4}\cdot 7^{4}\cdot(2^{8}-\pi^{8})<0.$

We observe that if a polynomial, say
$P(x)=\sum_{i=0}^{m}a_{i}x^{i}\in\mathbb{R}[x]$ of degree $m$ with its leading
coefficient $a_{m}\in\mathbb{R}_{<0}$, then $P(x)$ is a decreasing function in
$x$ and consequently, $P(x)<0$ for all $x\geq x_{0}$ where
$x_{0}\in\mathbb{R}$. So the only undetermined factor left over is the
explicit value of $y_{0}$ such that $N_{1}(y)D_{1}(y)<0$ for all $y\geq
y_{0}$, checked by Mathematica that $y_{0}=6$. We conclude the proof by
numerical verification that $s_{2}(n)<1$ holds for $3\leq n\leq 59$. ∎

Next, using the bound of $s_{2}(n)$ given in Theorem 2.6, we propose an upper
bound for $g(n)$ in terms of a function of $s(n)$ that enables us to get into
the proof of Theorem 2.5.

######  Theorem 2.7.

For $t\in(0,1)$, define

$\varphi(t)=\frac{1-\sqrt{1-t}}{t}.$ (2.42)

Then for all $n\geq 30$, we have

$g(n)<\varphi(s(n)).$ (2.43)

  * Proof.

We observe that for $t\in(0,1)$, $\varphi(t)$ is an increasing function in
$t$. From Corollary 2.4 and Theorem 2.6, it suggests that we need to prove for
$n\geq 91$,

$g(n)<\varphi(s_{1}(n)),$ (2.44)

or equivalently,

$s_{1}(n)g(n)^{2}-2g(n)+1>0.$ (2.45)

Recalling the definition of $\alpha(t)$ and $\beta(t)$ (cf. (2.33)),
$s_{1}(n)g(n)^{2}-2g(n)+1$ can be written in the following form

$s_{1}(n)g(n)^{2}-2g(n)+1=\frac{-g_{1}e^{r+w-2y}+g_{2}e^{w+2x-3y}+g_{3}-2g_{4}e^{x-2y+z}+g_{5}e^{r-3y+2z}}{r^{8}w^{8}x^{16}z^{16}(x^{6}-x^{5}-1)^{2}(y^{6}-y^{5}-1)^{4}(z^{6}-z^{5}-1)^{2}},$
(2.46)

with

$\displaystyle g_{1}$ $\displaystyle=$ $\displaystyle
x^{16}y^{16}z^{16}\alpha(r)\alpha(w)\alpha(x)^{2}\alpha(y)^{2}\alpha(z)^{2},$
(2.47) $\displaystyle g_{2}$ $\displaystyle=$ $\displaystyle
r^{8}y^{24}z^{16}\beta(w)\beta(x)^{2}\alpha(x)^{2}\beta(y)\beta(z)^{2},$
(2.48) $\displaystyle g_{3}$ $\displaystyle=$ $\displaystyle
r^{8}w^{8}x^{16}z^{16}\beta(x)^{2}\beta(y)^{4}\beta(z)^{2},$ (2.49)
$\displaystyle g_{4}$ $\displaystyle=$ $\displaystyle
r^{8}w^{8}x^{8}y^{16}z^{8}\alpha(x)\beta(x)^{2}\beta(y)^{2}\alpha(z)\beta(z)^{2},$
(2.50) $\displaystyle g_{5}$ $\displaystyle=$ $\displaystyle
w^{8}x^{16}y^{24}\beta(r)\beta(x)^{2}\beta(y)\beta(z)^{2}\alpha(z)^{2}.$
(2.51)

Since the denominator of (2.46) is a perfect square and hence positive,
therefore it is required to prove that

$G(y):=-g_{1}e^{r+w-2y}+g_{2}e^{w+2x-3y}+g_{3}-2g_{4}e^{x-2y+z}+g_{5}e^{r-3y+2z}>0.$
(2.52)

Following a similar method as used in Theorem 2.6, we first estimate the
exponential terms in (2.52). It is straightforward to observe that

$\displaystyle r_{2}+w_{2}-2y$ $\displaystyle=$
$\displaystyle-\frac{\pi^{4}(21\pi^{8}+10\pi^{4}y^{4}+8y^{8})}{8y^{11}}<0,$
$\displaystyle w_{1}+2x_{1}-3y$ $\displaystyle=$
$\displaystyle-\frac{\pi^{4}(17\pi^{10}+693\pi^{8}y^{2}-420\pi^{6}y^{4}+360\pi^{4}y^{6}-192\pi^{2}y^{8}+384y^{10})}{512y^{13}}<0\
\ \text{for}\ \ n\geq 1,$ $\displaystyle x_{2}-2y+z_{2}$ $\displaystyle=$
$\displaystyle-\frac{\pi^{4}(21\pi^{8}+40\pi^{4}y^{4}+128y^{8})}{512y^{11}}<0,$
$\displaystyle r_{1}+2z_{1}-3y$ $\displaystyle=$
$\displaystyle-\frac{\pi^{4}(1088\pi^{10}+693\pi^{8}y^{2}+420\pi^{6}y^{4}+360\pi^{4}y^{6}+192\pi^{2}y^{8}+384y^{10})}{512y^{13}}<0.$

As a consequence, by (2.19)-(2.22) and the monotonicity property of the
exponential function, for $n\geq 59$ we have

$e^{r+w-2y}<e^{r_{2}+w_{2}-2y}<\Phi(r_{2}+w_{2}-2y),$ (2.53)
$e^{x+z-2y}<e^{x_{2}+z_{2}-2y}<\Phi(x_{2}+z_{2}-2y),$ (2.54)
$e^{w+2x-3y}>e^{w_{1}+2x_{1}-3y}>\phi(w_{1}+2x_{1}-3y),$ (2.55)
$e^{r+2z-3y}>e^{r_{1}+2z_{1}-3y}>\phi(r_{1}+2z_{1}-3y).$ (2.56)

Substituting (2.53)-(2.56) into (2.52) implies that for $n\geq 59$,

$G(y)>-g_{1}\Phi(r_{2}+w_{2}-2y)+g_{2}\phi(w_{1}+2x_{1}-3y)+g_{3}-2g_{4}\Phi(x_{2}+z_{2}-2y)+g_{5}\phi(r_{1}+2z_{1}-3y).$
(2.57)

The right hand side of the above equation can be simplified further by
obtaining its lower bound with the aid of employing (Proof.) and (2.18) into
the definition of $\\{g_{\ell}\\}_{1\leq\ell\leq 5}$. More precisely, we have
that for $n\geq 59$,

$G(y)>\frac{N_{2}(y)}{D_{2}(y)}:=\frac{\sum_{i=0}^{223}b_{i}y^{i}}{2^{101}\cdot
3^{2}\cdot 5\cdot 7\ y^{119}},$ (2.58)

where $b_{223}=2^{101}\cdot 3^{2}\cdot 5\cdot 7$.
Due to similar remark as before; i.e, if a polynomial, say
$P(x)=\sum_{i=0}^{m}a_{i}x^{i}\in\mathbb{R}[x]$ of degree $m$ with its leading
coefficient $a_{m}\in\mathbb{R}_{>0}$, then $P(x)$ is an increasing function
in $x$ and consequently, $P(x)>0$ for all $x\geq x_{0}$ where
$x_{0}\in\mathbb{R}$. As an immediate consequence, we note that $G(y)>0$ by
verifying that $N_{2}(y)>0$ for all $y\geq 5$ or equivalently for $n\geq 3$.
It remains to prove (2.43) for $30\leq n\leq 91$ which is done by numerical
checking in Mathematica. ∎

##  3\. Proof of Theorem 2.5

_Proof of Theorem 2.5_. By corollary 2.4 and Theorem 2.6, we have $s(n)<1$ for
$n\geq 91$. Define

$Q(t)=s(n)t^{2}-2t+1.$ (3.1)

To establish (2.16), we prove that for $n\geq 91$,

$Q(u_{n})>0.$ (3.2)

The quadratic equation $Q(t)=0$ has two solutions, namely

$t_{0}=\frac{1-\sqrt{1-s(n)}}{s(n)},\ \ \text{and}\ \
t_{1}=\frac{1+\sqrt{1-s(n)}}{s(n)}.$

Thus $Q(t)>0$ when $t<t_{0}$ or $t>t_{1}$. From Theorem 2.3 and 2.7, we have
that for $n\geq 91$,

$u_{n}<\varphi(s(n))$ (3.3)

Set $t_{0}=\varphi(s(n))$ and conclude that (3.2) holds for $n\geq 91$. To
confirm (2.5) for $42\leq n\leq 91$, we can directly verify by checking
numerically in Mathematica.

##  4\. Conclusion

We conclude the paper by undertaking a brief study on totally positive
matrices with entries from sequences of overpartitions. Due to Engel [9], we
know that for $n\geq 2$,

$\text{det}\ M_{2}(\overline{p}(n)):=\text{det}\
\begin{pmatrix}\overline{p}(n)&\overline{p}(n+1)\\\
\overline{p}(n-1)&\overline{p}(n)\\\ \end{pmatrix}>0.$

Theorem 1.6 states that det $M_{3}(\overline{p}(n))>0$ for $n\geq 42$, more
specifically it is worthwhile to observe that for $n\geq 3$, the determinant
of each $2\times 2$ minor of the matrix $M_{3}(\overline{p}(n))$ is
nonnegative. This construction leads to the question whether one can always
construct a matrix $M_{k}(\overline{p}(n))$ of order $k$ with positive
determinant if we already know its all minors (of lower order) are totally
positive. More precisely,

######  Problem 4.1.

For a given $k\in\mathbb{Z}_{\geq 4}$, does there exists a
$n(k)\in\mathbb{Z}_{\geq 1}$ such that for $n>n(k)$,

$\det\ (\overline{p}(n-i+j))_{1\leq i,j\leq k}>0,$ (4.1)

and if (4.1) holds true, then what is the asymptotic growth of $n(k)?$

Acknowledgements

The author would like to express sincere gratitude to her advisor Prof. Manuel
Kauers for his valuable suggestions on the paper. The research was funded by
the Austrian Science Fund (FWF): W1214-N15, project DK13.

## References

  * [1] M. Aissen, A. Edrei, I. J. Schoenberg and A. M. Whitney.  On the generating functions of totally positive sequence.  Proc. Nat. Acad. Sci., 37: 303–307, 1951.
  * [2] F. Brenti.  Unimodal, $\log$-concave, and Pólya frequency sequences in combinatorics.  Mem. Amer. Math. Soc. 81, 1989.
  * [3] W. Y. C. Chen.  Recent developments on log-concavity and q-log-concavity of combinatorial polynomials.  22nd International Conference on Formal Power Series and Algebraic Combinatorics $($FPSAC $2010)$, 2010.http://www.billchen.org/talks/2010-FPSAC.pdf.
  * [4] W. Y. C. Chen, L. X. W. Wang and G. Y. B. Xie.  Finite differences of the logarithm of the partition function.  Math. Comput., 85:825–847, 2016.
  * [5] W. Y. C. Chen, D. X. Q. Jia and L. X. W. Wang.  Higher order Turán inequalities for the partition function.  Trans. Amer. Math. Soc., 372(3):2143–2165, 2019.
  * [6] S. Corteel and J. Lovejoy.  Overpartitions.  Trans. Am. Math. Soc. 356, 1623–1635, 2004.
  * [7] T. Craven and G. Csordas.  Jensen polynomials and the Turán and Laguerre inequalities.  Pacific J. Math., 136:241–260, 1989.
  * [8] S. DeSalvo and I. Pak.  Log-concavity of the partition function.  Ramanujan J., 38(1):61–73, 2015.
  * [9] B. Engel.  Log-concavity of the overpartition function.  Ramanujan J. 43(2), 229–241, 2017.
  * [10] M. Griffin, K. Ono, L. Rolen and D. Zagier.  Jensen polynomials for the Riemann zeta function and other sequences.  Proc. Natl. Acad. Sci. USA 116(23), 11103–11110, 2019.
  * [11] G. H. Hardy, J. E. Littlewood and G. Pólya.  Inequalities.  Cambridge University Press $($2d ed$)$, $1952$ Cambridge, England.
  * [12] G. H. Hardy, S. Ramanujan.  Asymptotic Formulae in Combinatory Analysis.  Proc. London Math. Soc. 17: 75–115, 1918.
  * [13] D. X. Jia and L. X. Wang.  Determinantal inequalities for the partition function.  Proceedings of the Royal Society of Edinburgh Section A: Mathematics, 150(3), 1451–1466, 2020.
  * [14] H. Larson and I. Wagner.  Hyperbolicity of the partition Jensen polynomials.  Res. Number Theory 5: 19, 2019.
  * [15] D. H. Lehmer.  On the remainders and convergence of the series for the partition function.  Trans. Amer. Math. Soc., 46:362–373, 1939.
  * [16] E. Y. S. Liu and H. W. J. Zhang.  Inequalities for the overpartition function.  Ramanujan J., 54(3): 485–509, 2021.
  * [17] J. Nicolas.  Sur les entiers N pour lesquels il y a beaucoup de groupes abéliens d’ordre N.  Ann. Inst. Fourier $(\text{Grenoble})$, 28(4), 1978.
  * [18] G. Pólya and J. Schur.  Über zwei Arten von Faktorenfolgen in der Theorie der algebraischen Gleichungen.  Reine Angew. Math., 144: 89–113, 1914.
  * [19] H. Rademacher.  A convergent series for the partition function $p(n)$.  Proc. Nat. Acad. Sci. 23, 78–84, 1937.
  * [20] L. X. W. Wang, G. Y. B. Xie, A. Q. Zhang.  Finite difference of the overpartition function.  Adv. Appl. Math. 92, 51–72, 2018.
  * [21] H. S. Zuckerman.  On the coefficients of certain modular forms belonging to subgroups of the modular group.  Trans. Am. Math. Soc. 45(2), 298–321, 1939.


# An Insight into the Dynamics and State Space Modelling of a 3-D Quadrotor

Rahul Vigneswaran K
Department of Mechanical Engineering
Amrita Vishwa Vidyapeetham
Amritapuri, India
<EMAIL_ADDRESS>&Soman KP
Center for Computational Engineering and Networking (CEN)
Amrita School of Engineering, Amrita Vishwa Vidyapeetham
Coimbatore, India
<EMAIL_ADDRESS>
https://rahulvigneswaran.github.ioWork done while interning at CEN, Amrita
Vishwa Vidyapeetham, India.

###### Abstract

Drones have gained popularity in a wide range of field ranging from aerial
photography, aerial mapping, and investigation of electric power lines. Every
drone that we know today is carrying out some kind of control algorithm at the
low level in order to manoeuvre itself around. For the quadrotor to either
control itself autonomously or to develop a high-level user interface for us
to control it, we need to understand the basic mathematics behind how it
functions. This paper aims to explain the mathematical modelling of the
dynamics of a 3 Dimensional quadrotor. As it may seem like a trivial task, it
plays a vital role in how we control the drone. Also, additional effort has
been taken to explain the transformations of the drone’s frame of reference to
the inertial frame of reference.

##  1 Introduction

The popularization of the drone technology has rendered the need for
understanding its dynamics and answering the question of _how to control
according it to our desirability?_ as vital. Drones are formally known as
Unmanned Aerial Vehicles (UAVs). It can either be controlled wirelessly by
pilots on base stations or can autonomously execute a pre-programmed task.
They are very in size and configurations. The drone studied in this paper
belongs to the category of micro-drones with 4 symmetrical rigid rotors. Even
though the quadrotors have become common, the underlying mathematics is still
complicated. This is so, because, quadrotors are _highly coupled_ , that is,
the fluid dynamics of propellers influence each other. Since its dynamics
involve multiple inputs and outputs, it is considered a multivariate problem
and call for the requirement of designing a non-linear mathematical model. It
must be noted that the system is actually under-actuated. The term under-
actuated means, when the number of actuators is less than the degree of
freedom. Here too, the drone can move along x,y,z axes and can rotate about
the x,y,z axes which makes up for a total of 6 DOF but the quadrotor, as the
name suggests is equipped with only 4 actuators, hence effectively making it
as an underactuated system. Further on how to overcome this and manoeuvre
through all 6 DOF will be discussed later on in this paper. Section 2 gives a
deeper understanding of rigid body transformations (both translational and
rotational). Section 3 introduces Euler angles from basics and move onto
derive a method to find the derivatives of them. Section 4 explains the
assumptions involved in the modelling and further proceeds to explaining the
kinematics model. Then it shows the design of a dynamics model from scratch by
discussing the various forces and aerodynamic moments that affects the model
and derives a state-space model for it.

##  2 Rigid body Transformations

Due to the reason that the drone is always in constant motion, it becomes
necessary to understand its current position with respect to its previous
position or how a set of mutually orthogonal vectors can be mapped to another
set of mutually orthogonal vectors. This can be achieved by using _Rigid body
transformations_ [1]. For the purpose of generalization, an arbitrary body is
taken into account instead of the quadrotor itself.

Figure 1: Arbitary points $p,q$ and their corresponding vectors _V,W_ are
being mapped on a rigid body.

For example, let us take the case that is shown in Figure 1. When a point in
$p$ is displaced, it is denoted by $g(p)$.111When a point is transformed, it
is denoted by $g$ We must note that displacement is essentially a
transformation of points. Similarly, $q$ is displaced to $g(q)$. Now consider
the vector V from $p$ to $q$. As the points $p$, $q$ moves to $g(p)$, $g(q)$,
the V is mapped to $g*(\textbf{\emph{V}})$.222When a vector is transformed, it
is denoted by $g_{*}$ A transformation is valid only when it satisfies the
following conditions.

  * •

A set of mutually orthogonal vectors after getting mapped must remain mutually
orthogonal.

  * •

The transformation $g$ must preserve lengths.

$\left\|{g(p)-g(q)}\right\|=\left\|{p-q}\right\|$

  * •

The transformation g* must preserve cross products.

${g_{*}}(\emph{V})\times{g_{*}}(\emph{W})={g_{*}}(\emph{V}\times\emph{W})$

Figure 2: Mapping of two rigid bodies $\\{A\\}$ and $\\{B\\}$ with different
frame of references.

Now let us take two rigid bodies $\\{A\\}$ and $\\{B\\}$ in a configuration as
shown in Figure 2 and try mapping them a point $P$ on ${B}$ to ${A}$. The body
${A}$ has $x-y-z$ frame of reference with origin $O$ attached to it while the
body $\\{B\\}$ has $x^{\prime}-y^{\prime}-z^{\prime}$ frame of reference with
origin $O^{\prime}$ attached to it. This is a two-stepped process. As we may
understand from Figure 2, the mapping involves a rotation mapping and
translation mapping. In order to map the $x-y-z$ frame to
$x^{\prime}-y^{\prime}-z^{\prime}$ frame, it must be rotated once and then
translated. The operations can be carried out the other way too.

${}^{A}{r^{P}}={}^{A}{R_{B}}{}^{B}{r^{P}}+{}^{A}{r^{O^{\prime}}}$

The equation (1) gives the transformation of point $P$ from frame of reference
of $\\{B\\}$ to $\\{A\\}$. ${}^{A}{\textbf{\emph{r}}^{P}}$ represents the
transformation of position vector of $P$ in $\\{B\\}$ to position vector of
$P$ in $\\{A\\}$ while ${}^{B}{\textbf{\emph{r}}^{P}}$ represents the relative
position of $\\{B\\}$ in reference to $\\{A\\}$.
${}^{A}{\textbf{\emph{r}}^{B}}$ represents the rotation of $\\{B\\}$ in
reference to $\\{A\\}$ and ${}^{A}{\textbf{\emph{r}}^{O^{\prime}}}$ represents
the translation of origin of $\\{B\\}$ in reference to $\\{A\\}$. Therefore,
the equation (1) first indicates the rotation transformation of $\\{B\\}$ in
reference to $\\{A\\}$ and then the translation transformation of the
resulting body.

###  2.1 Rotational Transformation

It is essential to understand the basic nature of the rotational
transformation to obtain intuitiveness for the rest of the paper. For the sake
of simplicity let’s take a 2-D configuration as shown in Figure 3 and later
extrapolate it to 3-D configuration. From it, the following components of the
$V$ and $V^{\prime}$ can be obtained.

Figure 3: A $x-y$ configuration in which a point from position $V$ is rotated
to position $V^{\prime}$.

${V_{x}}=r*\cos\alpha$ ${V_{y}}=r*\sin\alpha$
${{V^{\prime}}_{x}}=r*\cos(\alpha+\beta)$
${{V^{\prime}}_{y}}=r*\sin(\alpha+\beta)$

Now trigonometric identities are used to break down the equations (4) and (5)
to obtain the following,

${{V^{\prime}}_{x}}=r*\cos\alpha\cos\beta-r*\sin\alpha\sin\beta$
${{V^{\prime}}_{y}}=r*\sin\alpha\cos\beta-r*\cos\alpha\sin\beta$

Then by substituting equations (2) and (3) in equations (6) and (7), we get,

${{V^{\prime}}_{x}}={V_{x}}*\cos\beta-{V_{y}}*\sin\beta$
${{V^{\prime}}_{y}}={V_{y}}*\cos\beta-{V_{x}}*\sin\beta$

This when we convert into matrix notation, it gives the rotation of $V$ about
an arbitrary point $O$.

$\left[{\begin{array}[]{*{20}{c}}{{{V^{\prime}}_{x}}}&{{{V^{\prime}}_{y}}}\end{array}}\right]=\left[{\begin{array}[]{*{20}{c}}{{V_{x}}}&{{V_{y}}}\end{array}}\right]\left[{\begin{array}[]{*{20}{c}}{\cos\beta}&{\sin\beta}\\\
{-\sin\beta}&{\cos\beta}\end{array}}\right]$

Therefore, a 2-D rotation is defined by an angle and an origin. In 3-D, unlike
2-D, it is defined by an angle and a rotational axis. To obtain the rotational
transformation in 3-D, as the axis of rotation is restricted to one of the
three major axis, on component always remains constant. Then, by following the
same procedure as in 2-D, the rotational transformation about x-axis, y-axis
and z-axis are obtained.

The rotation around x-axis,

$\left[{\begin{array}[]{*{20}{c}}{{{V^{\prime}}_{x}}}&{{{V^{\prime}}_{y}}}&{{{V^{\prime}}_{z}}}\end{array}}\right]=\left[{\begin{array}[]{*{20}{c}}{{V_{x}}}&{{V_{y}}}&{{V_{z}}}\end{array}}\right]\left[{\begin{array}[]{*{20}{c}}1&0&0\\\
0&{\cos\beta}&{\sin\beta}\\\ 0&{-\sin\beta}&{\cos\beta}\end{array}}\right]$

The rotation around y-axis,

$\left[{\begin{array}[]{*{20}{c}}{{{V^{\prime}}_{x}}}&{{{V^{\prime}}_{y}}}&{{{V^{\prime}}_{z}}}\end{array}}\right]=\left[{\begin{array}[]{*{20}{c}}{{V_{x}}}&{{V_{y}}}&{{V_{z}}}\end{array}}\right]\left[{\begin{array}[]{*{20}{c}}{\cos\beta}&0&{-\sin\beta}\\\
0&1&0\\\ {\sin\beta}&0&{\cos\beta}\end{array}}\right]$

The rotation around z-axis,

$\left[{\begin{array}[]{*{20}{c}}{{{V^{\prime}}_{x}}}&{{{V^{\prime}}_{y}}}&{{{V^{\prime}}_{z}}}\end{array}}\right]=\left[{\begin{array}[]{*{20}{c}}{{V_{x}}}&{{V_{y}}}&{{V_{z}}}\end{array}}\right]\left[{\begin{array}[]{*{20}{c}}{\cos\beta}&{\sin\beta}&0\\\
{-\sin\beta}&{\cos\beta}&0\\\ 0&0&1\end{array}}\right]$

##  3 Euler angles

Leonhard Euler showed that any rotation can be described by three successive
rotations [2] about linearly independent axes and those three angles of
rotations are called as roll ($\phi$), pitch ($\theta$) and yaw ($\psi$).
Figure 4 illustrates these angles in reference to a glider. The Euler angles
$\phi$, $\theta$, $\psi$ can be determined from the rotational transformation
matrix. For the purpose of explanation let us take a $Z-Y-X$ rotational
transformation. This means that the body is first rotated angle about Z-axis,
then angle about Y-axis and then angle about X-axis. Equation (10) is the
resulting transformation matrix after $Z-Y-X$ rotation.

Figure 4: Roll($\phi$), pitch($\theta$), yaw($\psi$) in reference to a glider.

$\begin{gathered}R={R_{z}}(\psi){R_{y}}(\theta){R_{x}}(\phi)\hfill\\\
\quad=\left[{\begin{array}[]{*{20}{c}}{\cos\psi}&{\sin\psi}&0\\\
{-\sin\psi}&{\cos\psi}&0\\\
0&0&1\end{array}}\right]\left[{\begin{array}[]{*{20}{c}}{\cos\theta}&0&{-\sin\theta}\\\
0&1&0\\\
{\sin\theta}&0&{\cos\theta}\end{array}}\right]\left[{\begin{array}[]{*{20}{c}}1&0&0\\\
0&{\cos\phi}&{\sin\phi}\\\
0&{-\sin\phi}&{\cos\phi}\end{array}}\right]\hfill\\\
\quad=\left[{\begin{array}[]{*{20}{c}}{\cos\theta\cos\psi}&{\cos\psi\sin\theta\sin\phi-\cos\phi\sin\psi}&{\sin\phi\sin\psi+\cos\phi\cos\psi\sin\theta}\\\
{\cos\theta\sin\psi}&{\cos\phi\cos\psi+\sin\theta\sin\phi\sin\psi}&{\cos\phi\sin\theta\sin\psi-\cos\psi\sin\phi}\\\
{-\sin\theta}&{\cos\theta\sin\phi}&{\cos\theta\cos\phi}\end{array}}\right]\hfill\\\
\quad=\left[{\begin{array}[]{*{20}{c}}{{R_{11}}}&{{R_{12}}}&{{R_{13}}}\\\
{{R_{21}}}&{{R_{22}}}&{{R_{23}}}\\\
{{R_{31}}}&{{R_{32}}}&{{R_{33}}}\end{array}}\right]\hfill\\\ \end{gathered}$

Every element in the matrix (10) is a combination of two or more angles except
for $R_{31}$. By using it as starting equation, we get the following.

${R_{31}}=-\sin\theta$ $\theta=-{\sin^{-1}}({R_{31}})$

As we may understand from the equation (11), there are 2 possibilities of
$\theta$,

${\theta_{1}}=-{\sin^{-1}}({R_{31}})$
${\theta_{2}}=\pi-{\theta_{1}}=\pi+{\sin^{-1}}({R_{31}})$

From this step on, there are several cases that must be taken into account.

Case 1.1 : $(-1<\left|{{R_{31}}}\right|<1)$ and $(\cos\theta\neq 0)$

To find the value of $\phi$, we must observe that,

$\frac{{{R_{32}}}}{{{R_{33}}}}=\tan\phi$

Now we use this equation to solve for $\phi$ by using the
$\operatorname{atan2}$333An $\operatorname{atan2(y,x)}$ is an implementation
of $\operatorname{atan}(y,x)$ function that solves the sign ambiguity caused,
by taking the ratios and signs of $y$ and $x$ into consideration. Also,
$\operatorname{atan}(y,x)=\operatorname{atan}\left({\frac{y}{x}}\right)={\tan^{-1}}\left({\frac{y}{x}}\right)$.
function, as

$\phi=\operatorname{atan2}({R_{32}},{R_{33}})$

If $\cos\theta>0$, then the equation (14) holds. However, when $\cos\theta<0$,
then $\phi=\operatorname{atan2}(-{R_{32}},-{R_{33}})$. This can be bypassed by
using the equation as,

$\phi=\operatorname{atan2}\left({\frac{{{R_{32}}}}{{\cos\theta}},\frac{{{R_{33}}}}{{\cos\theta}}}\right)$

So, by substituting the equations (12) and (13) in (15), we get,

${\phi_{1}}=\operatorname{atan2}\left({\frac{{{R_{32}}}}{{\cos{\theta_{1}}}},\frac{{{R_{33}}}}{{\cos{\theta_{1}}}}}\right)$
${\phi_{2}}=\operatorname{atan2}\left({\frac{{{R_{32}}}}{{\cos{\theta_{2}}}},\frac{{{R_{33}}}}{{\cos{\theta_{2}}}}}\right)$

To find the corresponding values of $\psi$, similar analysis holds good. We
observe that,

$\frac{{{R_{21}}}}{{{R_{11}}}}=\tan\psi$

We solve for $\psi$ using an equation similar to (14),

$\psi=\operatorname{atan2}\left({\frac{{{R_{21}}}}{{\cos\theta}},\frac{{{R_{11}}}}{{\cos\theta}}}\right)$

Now, by substituting the equations (12) and (13) in (18), we get,

${\psi_{1}}=\operatorname{atan2}\left({\frac{{{R_{21}}}}{{\cos{\theta_{1}}}},\frac{{{R_{11}}}}{{\cos{\theta_{1}}}}}\right)$
${\psi_{2}}=\operatorname{atan2}\left({\frac{{{R_{21}}}}{{\cos{\theta_{2}}}},\frac{{{R_{11}}}}{{\cos{\theta_{2}}}}}\right)$

Therefore, for the case where $(-1<\left|{{R_{31}}}\right|<1$ and
$\cos\theta\neq 0)$, we get two sets of triplet Euler angles,

$\begin{gathered}({\phi_{1}},{\theta_{1}},{\psi_{1}})\hfill\\\
({\phi_{2}},{\theta_{2}},{\psi_{2}})\hfill\\\ \end{gathered}$

Both of these solutions are absolutely valid.

Case 1.2 : $(-1<\left|{{R_{31}}}\right|<1)$ and $(\cos\theta=0)$

In this case the elements ${R_{11}}$,${R_{21}}$,${R_{32}}$ and ${R_{33}}$ will
become zero causing problems444$\left({\frac{0}{0}}\right)$ is indeteminate as
it can hold any value. Therefore the equation cannot be solved to obtain the
values of $\phi$ and $\psi$. when solving the method used in case 1 and
doesn’t constrain the values $\phi$,$\psi$. As a result,equations (14) and
(18) becomes,

$\begin{gathered}\phi=\operatorname{atan}2\left({\frac{0}{0},\frac{0}{0}}\right)\hfill\\\
\psi=\operatorname{atan}2\left({\frac{0}{0},\frac{0}{0}}\right)\hfill\\\
\end{gathered}$

Therefore, we must use different elements of the rotational transformation
matrix (10) to find the solution.

Case 2.1 : $({R_{31}}=+1)$

Since we know that, $({R_{31}}=+1)$. From it, we can get
$\theta=\left({\frac{\pi}{2}}\right)$. Therefore, the rotational
transformation matrix (10) becomes,

$R=\left[{\begin{array}[]{*{20}{c}}0&{\cos\psi\sin\phi-\cos\phi\sin\psi}&{\sin\phi\sin\psi+\cos\phi\cos\psi}\\\
0&{\cos\phi\cos\psi+\sin\phi\sin\psi}&{\cos\phi\sin\psi-\cos\psi\sin\phi}\\\
{-1}&0&0\end{array}}\right]{\text{ }}$

As a result, $\psi$ and $\phi$ have an infinite number of valid solutions.

Case 2.2 : $({R_{31}}=-1)$

Similar to the case 2.1, since we know that, $({R_{31}}=-1)$. From it, we can
get $\theta=-\left({\frac{\pi}{2}}\right)$. Therefore, the rotational
transformation matrix (10) becomes,

$R=\left[{\begin{array}[]{*{20}{c}}0&{-\cos\psi\sin\phi-\cos\phi\sin\psi}&{\sin\phi\sin\psi-\cos\phi\cos\psi}\\\
0&{\cos\phi\cos\psi-\sin\phi\sin\psi}&{-\cos\phi\sin\psi-\cos\psi\sin\phi}\\\
1&0&0\end{array}}\right]{\text{ }}$

As a result, $\psi$ and $\phi$ have an infinite number of valid solutions. In
both the case 2.1 and 2.2, we have found that $\psi$ and $\phi$ are linked.
This phenomenon is called _Gimbal lock_. Although in this case, there is an
infinite number of solutions to the problem, in practice, one is often
interested in finding one solution. For this task, it is convenient to set
$\psi=0$ and compute $\phi$.

It is interesting to note that there is always more than one sequence of
rotations about the three principle axes that result in the same orientation
of an object. As an example, consider a book laying on a table face up in
front of you. Define the x-axis as to the right, the y-axis as away from you,
and the z-axis up. A rotation of $\pi$ radians about the y-axis will turn the
book so that the back cover is now facing up. Another way to achieve the same
orientation would be to rotate the book $\pi$ radians about the x-axis, and
then $\pi$ radians about the z-axis. Thus, there is more than one way to
achieve the desired rotation.

###  3.1 Derivative of Euler angles

In order to form the non-linear state space equation later on in this paper,
it is necessary that we derive the derivatives of Euler angles. This
subsection is dedicated to deriving an equation to obtain derivates of Euler
angles ($\dot{\phi},\dot{\theta},\dot{\psi}$) from angular velocities
($p,q,r$) while undergoing a $Z-Y-X$ rotational transformation. Angular
velocities in body frame can be expressed in terms of rate of change of the
Euler angles, as follows,

$\omega=p\hat{i}+q\hat{j}+r\hat{k}=\boldsymbol{\dot{\phi}}+\boldsymbol{\dot{\theta}}+\boldsymbol{\dot{\psi}}$

The first rotational transformation is
${X_{1}},{Y_{1}},{Z_{1}}\to{X_{2}},{Y_{2}},{Z_{2}}$ along Z-axis. In this
${Z_{1}}={Z_{2}}$ which implies that ${\hat{k}_{1}}={\hat{k}_{2}}$. Therefore,

$\boldsymbol{\dot{\psi}}=\dot{\psi}{\hat{k}_{1}}=\dot{\psi}{\hat{k}_{2}}$

The second rotational transformation is
${X_{2}},{Y_{2}},{Z_{2}}\to{X_{3}},{Y_{3}},{Z_{3}}$ along Y-axis. In this
${Y_{2}}={Y_{3}}$ which implies that ${\hat{j}_{2}}={\hat{j}_{3}}$. Therefore,

$\boldsymbol{\dot{\theta}}=\dot{\theta}{\hat{j}_{2}}=\dot{\theta}{\hat{j}_{3}}$

The third rotational transformation is
${X_{3}},{Y_{3}},{Z_{3}}\to{X},{Y},{Z}$555${X},{Y},{Z}$ refers to the
inertial/ground frame of reference. along Z-axis. In this ${X_{3}}={X}$ which
implies that ${\hat{i}_{3}}={\hat{i}}$. Therefore,

$\boldsymbol{\dot{\phi}}=\dot{\phi}{\hat{i}_{3}}=\dot{\phi}{\hat{i}}$

Now by combining equations (21), (22), (23) and (24), we get,

$\omega=p\hat{i}+q\hat{j}+r\hat{k}=\boldsymbol{\dot{\phi}}+\boldsymbol{\dot{\theta}}+\boldsymbol{\dot{\psi}}=\dot{\phi}\hat{i}+\dot{\theta}{\hat{j}_{3}}+\dot{\psi}{\hat{k}_{2}}$

In the transformation ${X_{2}},{Y_{2}},{Z_{2}}\to{X_{3}},{Y_{3}},{Z_{3}}$,

$\left[{\begin{array}[]{*{20}{c}}{{U_{2}}}\\\ {{V_{2}}}\\\
{{W_{2}}}\end{array}}\right]=\left[{\begin{array}[]{*{20}{c}}{\cos\theta}&0&{\sin\theta}\\\
0&1&0\\\
{-\sin\theta}&0&{\cos\theta}\end{array}}\right]\left[{\begin{array}[]{*{20}{c}}{{U_{3}}}\\\
{{V_{3}}}\\\ {{W_{3}}}\end{array}}\right]$

$\left[{\begin{array}[]{*{20}{c}}{{{\hat{i}}_{2}}}\\\ {{{\hat{j}}_{2}}}\\\
{{{\hat{k}}_{2}}}\end{array}}\right]=\left[{\begin{array}[]{*{20}{c}}{\cos\theta}&0&{\sin\theta}\\\
0&1&0\\\
{-\sin\theta}&0&{\cos\theta}\end{array}}\right]\left[{\begin{array}[]{*{20}{c}}{{{\hat{i}}_{3}}}\\\
{{{\hat{j}}_{3}}}\\\ {{{\hat{k}}_{3}}}\end{array}}\right]$

In the transformation ${X_{3}},{Y_{3}},{Z_{3}}\to X,Y,Z$,

$\left[{\begin{array}[]{*{20}{c}}{{U_{3}}}\\\ {{V_{3}}}\\\
{{W_{3}}}\end{array}}\right]=\left[{\begin{array}[]{*{20}{c}}1&0&0\\\
0&{\cos\phi}&{-\sin\phi}\\\
0&{\sin\phi}&{\cos\phi}\end{array}}\right]\left[{\begin{array}[]{*{20}{c}}U\\\
V\\\ W\end{array}}\right]$

$\left[{\begin{array}[]{*{20}{c}}{{{\hat{i}}_{3}}}\\\ {{{\hat{j}}_{3}}}\\\
{{{\hat{k}}_{3}}}\end{array}}\right]=\left[{\begin{array}[]{*{20}{c}}1&0&0\\\
0&{\cos\phi}&{-\sin\phi}\\\
0&{\sin\phi}&{\cos\phi}\end{array}}\right]\left[{\begin{array}[]{*{20}{c}}{\hat{i}}\\\
{\hat{j}}\\\ {\hat{k}}\end{array}}\right]$

Since, we already know that (26) and (28) are true, if we replace the velocity
vector with unit vectors like in (27) and (29), they still hold true.

From (27) we can find the expression for ${\hat{k}_{2}}$ in ${\hat{i}}$,
${\hat{j}}$, ${\hat{k}}$ and it is as follows,

${\hat{k}_{2}}=-\sin\theta{\hat{i}_{{}^{3}}}+\cos\theta{\hat{k}_{3}}$

But from (29) we know that, ${\hat{k}_{3}}=\sin\phi\hat{j}+\cos\phi\hat{k}$
and ${\hat{i}_{3}}=\hat{i}$, so,

${\hat{k}_{2}}=-\sin\theta\hat{i}+\cos\theta\sin\phi\hat{j}+\cos\theta\cos\phi\hat{k}$

Similarly, the expression for ${\hat{j}_{3}}$ in ${\hat{i}}$, ${\hat{j}}$,
${\hat{k}}$ can be derived from (29) and it is as follows,

${\hat{j}_{3}}=\cos\phi\hat{j}-\sin\phi\hat{k}$

Now by substituiting (30) and (31) in (25), we get,

$p\hat{i}+q\hat{j}+r\hat{k}=\dot{\phi}\hat{i}+\dot{\theta}(\cos\phi\hat{j}-\sin\phi\hat{k})+\dot{\psi}(-\sin\theta\hat{i}+\cos\theta\sin\phi\hat{j}+\cos\theta\cos\phi\hat{k})$

From (32), we can derive the following values,

$\begin{gathered}p=\dot{\phi}-(\sin\theta)\dot{\psi}\hfill\\\
q=(\cos\phi)\dot{\theta}+(\cos\theta\sin\phi)\dot{\psi}\hfill\\\
r=(\cos\theta\cos\phi)\dot{\psi}-(\sin\phi)\dot{\theta}\hfill\\\
\end{gathered}$

By rearranging the above in matrix form, we get,

$\left[{\begin{array}[]{*{20}{c}}p\\\ q\\\
r\end{array}}\right]=\left[{\begin{array}[]{*{20}{c}}1&0&{-\sin\theta}\\\
0&{\cos\phi}&{\cos\theta\sin\phi}\\\
0&{-\sin\phi}&{\cos\theta\cos\phi}\end{array}}\right]\left[{\begin{array}[]{*{20}{c}}{\dot{\phi}}\\\
{\dot{\theta}}\\\ {\dot{\psi}}\end{array}}\right]$

Or,

$\left[{\begin{array}[]{*{20}{c}}{\dot{\phi}}\\\ {\dot{\theta}}\\\
{\dot{\psi}}\end{array}}\right]=\left[{\begin{array}[]{*{20}{c}}1&{\sin\phi\tan\theta}&{\cos\phi\tan\theta}\\\
0&{\cos\phi}&{-\sin\phi}\\\ 0&{\sin\phi\sec\theta}&{\cos\phi
sec\theta}\end{array}}\right]\left[{\begin{array}[]{*{20}{c}}p\\\ q\\\
r\end{array}}\right]$

The equation (33) can be used to obtain the derivative of Euler angles when we
know the angular velocities.

##  4 System Modelling of Quadrotor

Let’s us take the configuration shown in Figure 5 for modelling the system. In
this, the propellers 1 and 3 rotate in the same direction while propellers 3
and 4 rotate in opposite direction. This is done so, to overcome the imbalance
caused moments of the rotating propellers. A quadrotor is an underactuated
system. It has only 4 rotors facing in the same direction in order to
manoeuvre through all 6 DOFs. This underactuation is tackled by executing two
or more operations (roll, pitch or yaw) at a time to move in the desired
direction. Let us take the task of a quadrotor moving forward. When a
quadrotor is hovering and is stationary, it has to produce enough overall
thrust force (say $F$) to equalize the gravitational force (say $G$) acting on
it, that is $F=G$. So, in order to move forward, it has to make a pitch
$\theta$. This creates a forward horizontal component $Fsin\theta$ as shown in
Figure 6 which enables the quadrotor to move forward. But doing so reduces the
vertical component to $Fcos\theta$. As we know from earlier, $Fcos\theta\neq
G$. So, the quadrotor will lose altitude and overcome it, we must increase the
thrust force F accordingly to both move forward and not lose altitude.
Similarly, all the other operations are executed with a combination of roll,
pitch and yaw rotations to overcome the hurdle of underactuation. There are a
few assumptions taken into account while modelling the configuration given in
Figure 5 to make the modelling easier and they are :

  * •

It is rigid, that is there are no other moving parts other than the rotor.

  * •

It has a symmetrical design.

  * •

The propellers are equally spaced from the centre of the body.

  * •

Various masses present in quadrotor (battery, circuit boards, etc.) are
distributed equally with a centre of mass coinciding with the frame’s centre.

Figure 5: Force diagram of a conventional quadrotor. Figure 6: Thrust force
diagram of a quadrotor moving in forward direction.

###  4.1 Kinematic Model

The first step to model a system is to define the coordinate frames and do the
necessary transformations. As shown in Figure 5, the earth inertial frame axes
are N, E, D and the body frame axes are X, Y, Z. The inertial frame is fixed
to a specific place and the notations N-E-D represent North, East and
Downwards respectively. For the body frame of reference, the origin is located
at the centre of the quadrotor’s body with X, Y and Z axes pointing towards
propellers 1, 2 and ground respectively. The orientation of the quadcopter is
defined by the three Euler angles ($\phi,\theta,\psi$).

The thrust force is defined from the body from of reference and we need a
transformation matrix (${}^{W}{R_{B}}$ \- World frame of reference to Body
frame of reference) to obtain that. First the rotation transformation is done
about Z-axis followed by Y-axis and X-axis. By doing so, we get ,

$\begin{gathered}{}^{W}{R_{B}}={R_{z}}(\psi){R_{y}}(\theta){R_{x}}(\phi)\hfill\\\
\quad=\left[{\begin{array}[]{*{20}{c}}{\cos\psi}&{\sin\psi}&0\\\
{-\sin\psi}&{\cos\psi}&0\\\
0&0&1\end{array}}\right]\left[{\begin{array}[]{*{20}{c}}{\cos\theta}&0&{-\sin\theta}\\\
0&1&0\\\
{\sin\theta}&0&{\cos\theta}\end{array}}\right]\left[{\begin{array}[]{*{20}{c}}1&0&0\\\
0&{\cos\phi}&{\sin\phi}\\\
0&{-\sin\phi}&{\cos\phi}\end{array}}\right]\hfill\\\
\quad=\left[{\begin{array}[]{*{20}{c}}{\cos\theta\cos\psi}&{\cos\psi\sin\theta\sin\phi-\cos\phi\sin\psi}&{\sin\phi\sin\psi+\cos\phi\cos\psi\sin\theta}\\\
{\cos\theta\sin\psi}&{\cos\phi\cos\psi+\sin\theta\sin\phi\sin\psi}&{\cos\phi\sin\theta\sin\psi-\cos\psi\sin\phi}\\\
{-\sin\theta}&{\cos\theta\sin\phi}&{\cos\theta\cos\phi}\end{array}}\right]\hfill\\\
\end{gathered}$

But we need a transformation from the body frame of reference to world frame
of reference which can be obtained by inverting ${}^{W}{R_{B}}$. By doing so,
we get,

${}^{B}{R_{W}}=\left[{\begin{array}[]{*{20}{c}}{\cos\theta\cos\psi}&{\cos\theta\sin\psi}&{-\sin\theta}\\\
{\cos\psi\sin\theta\sin\phi-\cos\phi\sin\psi}&{\cos\phi\cos\psi+\sin\theta\sin\phi\sin\psi}&{\cos\theta\sin\phi}\\\
{\sin\phi\sin\psi+\cos\phi\cos\psi\sin\theta}&{\cos\phi\sin\theta\sin\psi-\cos\psi\sin\phi}&{\cos\theta\cos\phi}\end{array}}\right]{\text{
}}$

###  4.2 Dynamic Model

The quadrotor’s motion can be divided into two subsystems - translational
motion (X, Y and Z) and rotational motion (roll, pitch and yaw). The forces
and moments acting on the quadrotor will be investigated using Newton-Euler
formalism. In Figure 5, $F_{1}$, $F_{2}$, $F_{3}$ and $F_{4}$ are the thrust
forces acting on the quadrotor while $M_{1}$, $M_{2}$, $M_{3}$ and $M_{4}$ are
the moments caused by the rotation of the propellers. Also, ${\tau_{x}}$,
${\tau_{y}}$ and ${\tau_{z}}$ are the torques caused by the thrust forces
about X, Y and Z axes respectively.

In some places, for the sake of alignment,
$\left[{\begin{array}[]{*{20}{c}}A\\\ B\\\ C\end{array}}\right]$ is written as
${\left[{\begin{array}[]{*{20}{c}}A&B&C\end{array}}\right]^{T}}$.

####  4.2.1 Translational Equations of Motion

According to Newton’s law,

$\begin{gathered}{{\vec{F}}_{\omega}}=\frac{d}{{dt}}(m.\vec{v})\\\
=m.\dot{v}\\\ \end{gathered}$

where, ${\vec{F}_{\omega}}=\vec{G}-\vec{T}$. $G$ is the gravity and $T$ is the
thrust generated by the propellers. As $\vec{T}$ is described in body frame of
reference, equation (35) is used to transform it to world frame of reference
and the expression (36) can be expressed as,

$\dot{v}=\frac{1}{m}\left[{\left[{\begin{array}[]{*{20}{c}}0\\\ 0\\\
{mg}\end{array}}\right]-{}^{B}{R_{W}}\left[{\begin{array}[]{*{20}{c}}0\\\ 0\\\
T\end{array}}\right]}\right]$

Now substitute (35) in (37), to get,

$\dot{v}=\frac{1}{m}\left[{\begin{array}[]{*{20}{c}}{-T(\cos\phi\sin\theta\cos\psi+\sin\phi\sin\psi)}\\\
{-T(\cos\phi\sin\theta\sin\psi-\sin\phi\cos\psi)}\\\
{mg-T(\cos\phi\cos\theta)}\end{array}}\right]$

This vector equation can be written in terms of its components as,

$\begin{gathered}\ddot{x}=-\frac{T}{m}(\cos\phi\sin\theta\cos\psi+\sin\phi\sin\psi)\\\
\ddot{y}=-\frac{T}{m}(\cos\phi\sin\theta\sin\psi-\sin\phi\cos\psi)\\\
\ddot{z}=g-\frac{T}{m}(\cos\phi\cos\theta)\\\ \end{gathered}$

####  4.2.2 Rotational Equations of Motion

In the inertial frame, the moment can be defined as the time derivative of the
angular momentum and can be expressed as follows,

${\vec{M}_{\omega}}=\frac{{d\vec{L}}}{{dt}}=\frac{d}{{dx}}(\vec{I}\bullet\vec{\omega})$

where, $\vec{I}$ is the body’s inertia tensor.

$\vec{I}=\left[{\begin{array}[]{*{20}{c}}{{I_{xx}}}&{{I_{xy}}}&{{I_{xz}}}\\\
{{I_{yx}}}&{{I_{yy}}}&{{I_{yz}}}\\\
{{I_{zx}}}&{{I_{zy}}}&{{I_{zz}}}\end{array}}\right]=\left[{\begin{array}[]{*{20}{c}}{{I_{xx}}}&0&0\\\
0&{{I_{yy}}}&0\\\ 0&0&{{I_{zz}}}\end{array}}\right]]$

In equation (41), the off-diagonal terms are zero due to the symmetric mass
distribution assumption as mentioned in Section 4.

According to Euler’s rotational equation,

${\vec{M}_{b}}=\left({\vec{I}\bullet\dot{\vec{\omega}}}\right)+\vec{\omega}\times\left({\vec{I}\bullet\vec{\omega}}\right)=\left[{\begin{array}[]{*{20}{c}}{{\tau_{x}}}\\\
{{\tau_{y}}}\\\ {{\tau_{z}}}\end{array}}\right]$

where, $\dot{\vec{\omega}}=\left[{\begin{array}[]{*{20}{c}}{\dot{p}}\\\
{\dot{q}}\\\ {\dot{r}}\end{array}}\right]$ and the right hand side is a vector
of applied torques.

Also, by simplifying
$\vec{\omega}\times\left({\vec{I}\bullet\vec{\omega}}\right)$, we get,

$\vec{\omega}\times\left({\vec{I}\bullet\vec{\omega}}\right)=\left[{\begin{array}[]{*{20}{c}}p\\\
q\\\ r\end{array}}\right]\times\left[{\begin{array}[]{*{20}{c}}{{I_{xx}}p}\\\
{{I_{yy}}q}\\\
{{I_{zz}}r}\end{array}}\right]=\left[{\begin{array}[]{*{20}{c}}{({I_{zz}}-{I_{yy}})qr}\\\
{({I_{yy}}-{I_{zz}})rp}\\\ {({I_{yy}}-{I_{xx}})pq}\end{array}}\right]$

Now, by substituting (43) in (42), we get,

$\begin{gathered}{\tau_{x}}={I_{xx}}\dot{p}+({I_{zz}}-{I_{yy}})qr\\\
{\tau_{y}}={I_{yy}}\dot{q}+({I_{yy}}-{I_{zz}})rp\\\
{\tau_{z}}={I_{zz}}\dot{r}+({I_{yy}}-{I_{xx}})pq\\\ \end{gathered}$

Here$p,q,r$ are actually equal to $\dot{\phi},\dot{\theta},\dot{\psi}$ since
by definition, angular velocity of the body are same as time derivatives of
Euler angles of it .When the eqution (44) is rearranged, we get,

$\begin{gathered}\dot{p}=\ddot{\phi}=({\tau_{x}}+{I_{yy}}qr-{I_{zz}}qr)/{I_{xx}}\\\
\dot{q}=\ddot{\theta}=({\tau_{y}}-{I_{xx}}pr+{I_{xx}}pr)/{I_{yy}}\\\
\dot{r}=\ddot{\psi}=({\tau_{z}}+{I_{xx}}pq-{I_{yy}}pq)/{I_{zz}}\\\
\end{gathered}$

Using equation (33), we can get the derivatives of the Euler angles as
follows,

$\begin{gathered}\dot{\phi}=p+{\text{
}}(r\cos(\phi)\sin(\theta))/\cos(\theta){\text{ }}+{\text{
}}(q\sin(\theta)\sin(\phi))/\cos(\theta)\hfill\\\
\dot{\theta}=q\cos(\phi)-r\sin(\phi)\hfill\\\ \dot{\psi}={\text{
}}(r\cos(\phi))/\cos(\theta){\text{ }}+{\text{
}}(q\sin(\phi))/\cos(\theta)\hfill\\\ \end{gathered}$

####  4.2.3 Aerodynamic Moment and thrust Force

As the propeller rotates, it generates a thrust force. This thrust force is
directly proportional [4] to the square of the rotor speed and can be
expressed as,

${F_{i}}=-{K_{a}}\omega_{i}^{2}$

Apart from thrust force, the rotating propellers generate aerodynamic moment
which is also directly proportional to the square of the rotor speed and can
be expressed,

${M_{i}}=-{K_{m}}\omega_{i}^{2}$

For the case of low altitude flight, the values of $K_{a}$ and $K_{m}$ in
expression (47) and (48) can be considered as constants. The directions of the
thrust force and aerodynamic moments generated by the propellers is
illustrated in the Figure 5. Now, we can take the moment caused by the thrust
force generated into account. Let’s say that the length of the arms reaching
out to the propellers from the centre of the body is $l$. Now, the total
moment caused by thrust force about the X-axis is,

$\begin{gathered}{\tau_{x}}=-{F_{2}}l+{F_{4}}l\hfill\\\
=-({K_{a}}\omega_{2}^{2})d+({K_{a}}\omega_{4}^{2})l\hfill\\\
={K_{a}}l(\omega_{4}^{2}-\omega_{2}^{2})\hfill\\\ \end{gathered}$

Total moment caused by the generated thrust about Y-axis is,

$\begin{gathered}{\tau_{y}}={F_{1}}l-{F_{3}}l\\\
=({K_{a}}\omega_{1}^{2})l-({K_{a}}\omega_{3}^{2})l\\\
={K_{a}}l(\omega_{1}^{2}-\omega_{3}^{2})\\\ \end{gathered}$

The total moment about Z-axis is influenced only by the aerodynamic moment
caused by the propellers (48) but not affected by the thrust force and it can
be expressed as,

$\begin{gathered}{\tau_{z}}={M_{1}}-{M_{2}}+{M_{3}}-{M_{4}}\\\ ={\text{
}}({K_{m}}\omega_{1}^{2})-({K_{m}}\omega_{2}^{2}){\text{
}}+({K_{m}}\omega_{3}^{2})-({K_{m}}\omega_{4}^{2})\\\
={K_{m}}(\omega_{1}^{2}-\omega_{2}^{2}+\omega_{3}^{2}-\omega_{4}^{2})\\\
\end{gathered}$

The total thrust force can be expressed as,

$T=-{K_{a}}(\omega_{1}^{2}+\omega_{2}^{2}+\omega_{3}^{2}+\omega_{4}^{2})$

The negative sign in expression (52) indicates that the total thrust force is
acting in the negative direction of Z-axis. Also, in this quadrotor
configuration, thrust force exists only in the direction of Z-axis.

###  4.3 State Space Equation

In general, the evolution of a dynamical system’s states are governed by set
of ordinary differential equations. In control problems, these ordinary
differential equations are often rearranged into state-space form,

$\dot{X}=f(X,U)$

where, $X$ is a matrix of states and $U$ is a matrix of inputs. The formation
of state space equation involves 3 steps :

  * •

Identify the governing equations and and find the order ($n$) of it.

  * •

Define the state and input vector.

  * •

Define the system of first-order differential equations.

####  4.3.1 Governing Equations

Equations (39) and (45) are the governing equations and are as follows,

$\begin{gathered}\ddot{x}=-\frac{T}{m}(\cos\phi\sin\theta\cos\psi+\sin\phi\sin\psi)\\\
\ddot{y}=-\frac{T}{m}(\cos\phi\sin\theta\sin\psi-\sin\phi\cos\psi)\\\
\ddot{z}=g-\frac{T}{m}(\cos\phi\cos\theta)\\\ \end{gathered}$

$\begin{gathered}\dot{p}=\ddot{\phi}=({\tau_{x}}+{I_{yy}}qr-{I_{zz}}qr)/{I_{xx}}\\\
\dot{q}=\ddot{\theta}=({\tau_{y}}-{I_{xx}}pr+{I_{xx}}pr)/{I_{yy}}\\\
\dot{r}=\ddot{\psi}=({\tau_{z}}+{I_{xx}}pq-{I_{yy}}pq)/{I_{zz}}\\\
\end{gathered}$

The order ($n$) of the all the governing equation is 2.

####  4.3.2 State and Input Vector

The order of each governing equation is 2. Therefore, each governing equation
needs $(n-1)$ states, that is 2 states each. Therefore, in order to define the
motion of a quadrotor, practically, we will need $(2\times 6)$ states, that is
12 states -
$x,y,z,\dot{x},\dot{y},\dot{z},\phi,\theta,\psi,\dot{\phi},\dot{\theta},\dot{\psi}$.
Therefore the matrix of the states of the system is,

$\begin{gathered}X={\left[{\begin{array}[]{*{20}{c}}{{x_{1}}}&{{x_{2}}}&{{x_{3}}}&{{x_{4}}}&{{x_{5}}}&{{x_{6}}}&{{x_{7}}}&{{x_{8}}}&{{x_{9}}}&{{x_{10}}}&{{x_{11}}}&{{x_{12}}}\end{array}}\right]^{T}}\hfill\\\
={\left[{\begin{array}[]{*{20}{c}}x&y&z&\phi&\theta&\psi&{\dot{x}}&{\dot{y}}&{\dot{z}}&{\dot{\phi}}&{\dot{\theta}}&{\dot{\psi}}\end{array}}\right]^{T}}\hfill\\\
\end{gathered}$

The matrix of the inputs to the system is,

$\begin{gathered}U=\left[{\begin{array}[]{*{20}{c}}{{U_{1}}}&{{U_{2}}}&{{U_{3}}}&{{U_{4}}}\end{array}}\right]\hfill\\\
=\left[{\begin{array}[]{*{20}{c}}{\omega_{1}^{2}}&{\omega_{2}^{2}}&{\omega_{3}^{2}}&{\omega_{4}^{2}}\end{array}}\right]\hfill\\\
\end{gathered}$

####  4.3.3 System of first-order equations

The first 6 equations simply related states to each other by finding their
derivatives like,

$\begin{gathered}\frac{{d{x_{1}}}}{{dt}}=\frac{{dx}}{{dt}}=\dot{x}={x_{7}}\\\
\frac{{d{x_{2}}}}{{dt}}=\frac{{dy}}{{dt}}=\dot{y}={x_{8}}\\\
\frac{{d{x_{3}}}}{{dt}}=\frac{{dz}}{{dt}}=\dot{z}={x_{9}}\\\
\frac{{d{x_{4}}}}{{dt}}=\frac{{d\phi}}{{dt}}=\dot{\phi}={x_{10}}\\\
\frac{{d{x_{5}}}}{{dt}}=\frac{{d\theta}}{{dt}}=\dot{\theta}={x_{11}}\\\
\frac{{d{x_{6}}}}{{dt}}=\frac{{d\psi}}{{dt}}=\dot{\psi}={x_{12}}\\\
\end{gathered}$

The next 6 equations can be derived by differentiating and rearranging the
above 6 first-order differential equations like,

$\begin{gathered}\frac{{d{x_{7}}}}{{dt}}=\frac{{d\dot{x}}}{{dt}}=\ddot{x}\\\
\frac{{d{x_{8}}}}{{dt}}=\frac{{d\dot{y}}}{{dt}}=\ddot{y}\\\
\frac{{d{x_{9}}}}{{dt}}=\frac{{d\dot{z}}}{{dt}}=\ddot{z}\\\
\frac{{d{x_{10}}}}{{dt}}=\frac{{d\dot{\phi}}}{{dt}}=\ddot{\phi}\\\
\frac{{d{x_{11}}}}{{dt}}=\frac{{d\dot{\theta}}}{{dt}}=\ddot{\theta}\\\
\frac{{d{x_{12}}}}{{dt}}=\frac{{d\dot{\psi}}}{{dt}}=\ddot{\psi}\\\
\end{gathered}$

Now, by substituting the values of values for
$\ddot{x},\ddot{y},\ddot{z},\dot{\phi},\dot{\theta},\dot{\psi},\ddot{\phi},\ddot{\theta},\ddot{\psi}$,
from the the equations (39), (45) and (46) onto the above 12 first-order
differential equationa=s, we get our state space equation in the form
$\\{\dot{X}=f(X,U)\\}$ as,

$\left[{\begin{array}[]{*{20}{c}}{{\dot{X}_{1}}}\\\ {{\dot{X}_{2}}}\\\
{{\dot{X}_{3}}}\\\ {{\dot{X}_{4}}}\\\ {{\dot{X}_{5}}}\\\ {{\dot{X}_{6}}}\\\
{{\dot{X}_{7}}}\\\ {{\dot{X}_{8}}}\\\ {{\dot{X}_{9}}}\\\ {{\dot{X}_{10}}}\\\
{{\dot{X}_{11}}}\\\
{{\dot{X}_{12}}}\end{array}}\right]=\left[{\begin{array}[]{*{20}{c}}{{f_{1}}}\\\
{{f_{2}}}\\\ {{f_{3}}}\\\ {{f_{4}}}\\\ {{f_{5}}}\\\ {{f_{6}}}\\\ {{f_{7}}}\\\
{{f_{8}}}\\\ {{f_{9}}}\\\ {{f_{10}}}\\\ {{f_{11}}}\\\
{{f_{12}}}\end{array}}\right]=\left[{\begin{array}[]{*{20}{c}}{\dot{x}}\\\
{\dot{y}}\\\ {\dot{z}}\\\ {\dot{\phi}}\\\ {\dot{\theta}}\\\ {\dot{\psi}}\\\
{\ddot{x}}\\\ {\ddot{y}}\\\ {\ddot{z}}\\\ {\ddot{\phi}}\\\ {\ddot{\theta}}\\\
{\ddot{\psi}}\end{array}}\right]=\left[{\begin{array}[]{*{20}{c}}{{x_{7}}}\\\
{{x_{8}}}\\\ {{x_{9}}}\\\ {p+{\text{
}}(r\cos(\phi)\sin(\theta))/\cos(\theta){\text{ }}+{\text{
}}(q\sin(\theta)\sin(\phi))/\cos(\theta)}\\\ {q\cos(\phi)-r\sin(\phi)}\\\
{(r\cos(\phi))/\cos(\theta){\text{ }}+{\text{ }}(q\sin(\phi))/\cos(\theta)}\\\
{-\frac{T}{m}(\cos\phi\sin\theta\cos\psi+\sin\phi\sin\psi)}\\\
{-\frac{T}{m}(\cos\phi\sin\theta\sin\psi-\sin\phi\cos\psi)}\\\
{g-\frac{T}{m}(\cos\phi\cos\theta)}\\\
{({\tau_{x}}+{I_{yy}}qr-{I_{zz}}qr)/{I_{xx}}}\\\
{({\tau_{y}}-{I_{xx}}pr+{I_{xx}}pr)/{I_{yy}}}\\\
{({\tau_{z}}+{I_{xx}}pq-{I_{yy}}pq)/{I_{zz}}}\end{array}}\right]$

Therefore, the non-linear state-space equation for a 3-D quadrotor has been
derived. It is non-linear due to the presence $sine$ and $cosine$ terms in it.

##  5 Conclusion

Dynamics of a 3-D quadrotor and its state-space equation has been derived from
the grass root level with detailed explanation and the logical reasoning
behind every step have been stated clearly. As a future work, the dynamical
system that is modelled in this paper can be used in conjunction with concepts
of control system theory to make filters and controllers that can make the
quadrotor follow a trajectory or tackle any kind of task with precision and
robustness using MATLAB. Apart from that, the concepts of Quaternions [5] can
be used in the modelling of the system to make the computations easier.

#### Acknowledgments

This paper includes a collective knowledge majorly from the lectures of
Dr.Vijay Kumar and Sarah Tang from University of Pennsylvania and the master
thesis work [3] of Haifeng Lyu from University of Rhode Island. Also, this
paper would have not been possible without the continued guidance of Dr. KP.
Soman of CEN, Amrita Vishwa Vidyapeetham, India.

## References

  * [1] Roger W Brockett.  Robotic manipulators and the product of exponentials formula.  In Mathematical theory of networks and systems, pages 120–129. Springer, 1984.
  * [2] Leonhard Euler.  Formulae generales pro translatione quacunque corporum rigidorum.  Novi Commentarii academiae scientiarum Petropolitanae, pages 189–207, 1776.
  * [3] Haifeng Lyu.  Multivariable control of a rolling spider drone.  2017\.
  * [4] Paul Pounds, Robert Mahony, Joel Gresham, Peter Corke, and Jonathan M Roberts.  Towards dynamically-favourable quad-rotor aerial robots.  In Proceedings of the 2004 Australasian Conference on Robotics & Automation. Australian Robotics & Automation Association, 2004.
  * [5] Stanley Wayne Shepperd.  Quaternion from rotation matrix.[four-parameter representation of coordinate transformation matrix].  1978\.