Markov Chains - A Mathematical Treatise

Markov Chains

We have a set of stated, $\mathcal{S}$ = {${s_1, s_2, ..., s_m}$}. A Markov process starts in one of these states and moves successively from one state to another. Each move is called a step . If the chain is currently in state $s_i$, then it moves to state $s_j$ at the next step with a probability denoted by $p_{ij}$, and this probability does not depend upon which states the chain was in before the current state.These probabilities $p_{ij}$ are called transition probabilities. The process can remain in the state it is in, and this occurs with probability $p_{ii}$ . Initially a probability is defined by specifying a particular state as the starting state.

Say, we stay in Patna (Capital of the State of Bihar, India), and at my home we usually eat either Roti ($\text{Ro}$), Rice ($\text{Ri}$), Parantha ($\text{Pa}$)) or Bread ($\text{Br}$) for dinner. There are also some days when we eat one particular dinner on two or three or more consecutive days.

Now, we define transition probabilities of using a Transition Matrix which is given as,
$$\text{T}=\begin{array}{c c} & \begin{array}{c c c} \text{Ri} & \text{Ro} & \text{Pa} & \text{Br} \\ \end{array} \\ \begin{array}{c c c} \text{Ri} \\ \text{Ro} \\ \text{Pa} \\ \text{Br}\\ \end{array} & \left[ \begin{array}{c c c} \frac{1}{2} & \frac{1}{4} & \frac{1}{4} & 0 \\ 0 & \frac{1}{4} & \frac{1}{4} & \frac{1}{2} \\ 0 & \frac{1}{2} & \frac{1}{4} & \frac{1}{4} \\ \frac{1}{2} & 0 & \frac{1}{4} & \frac{1}{4} \\ \end{array} \right] \end{array}$$
From the above matrix, it can be seen that $P_{(Rice \implies Rice)}=\frac{1}{2}$ and similarly, $P_{(Roti \implies Parantha)}=\frac{1}{4}$.

In the similar lines of the above concept, a Transition Diagram can also be drawn.

Let us dig further into the concept now. Imagine today we eat Rice at home and two days from today we eat Bread. This is written as $P^{2}_{(Rice \implies Bread)}$. We see that if we ate Rice today then the event of us eating Bread two days from today is a disjoint union of the following,

1. We ate Rice tomorrow and Bread day after.
2. We ate something else tomorrow and Bread day after.
3. We ate Bread tomorrow and Bread day after.

Therefore, we can write $P^{2}_{(Rice \implies Bread)}$ as ,
$$P^{2}_{(Rice \implies Bread)} = [P_{(Rice \implies Rice)}]\cdot[P_{(Rice \implies Bread)}]+[P_{(Rice \implies Roti)}]\cdot[P_{(Roti \implies Bread)}]+ \\ [P_{(Rice \implies Parantha)}]\cdot[P_{(Parantha \implies Bread)}]+[P_{(Rice \implies Bread)}]\cdot[P_{(Bread \implies Bread)}]$$
The has pointed towards a more generalized concept of Dot product of two vectors.

Considering there are $r$ states in the Markov chain,
$$P^{2}_{ij}=\sum_{k=1}^{r}P_{ik}P_{kj}$$
This study was for $n=2$ , meaning we are yet to define this in a more generalized manner. Obviously, looking at the present scenario of high-end computing software packages available, this attempt seems to be a waste of efforts, but then knowing the math is always good fun!

Before doing that, we will again go into the basic definition as defined above and try to make an iterative approach into the problem.

We have assumed that our meal starts at state Rice, so, for us,
$$\text{p}_0=\begin{bmatrix} 1 \\ 0 \\ 0 \\ 0\\ \end{bmatrix}$$
$\text{p}_0$ is the state of our system at $0^{th}$ state or the beginning of the system.

Similarly, if we want to find $\text{p}_1$ we will do it as,
$$\text{p}_1 = \text{T} \cdot \text{p}_0$$
In the similar lines,
$$\boxed{\text{p}_5 = \text{T} \cdot \text{p}_4 = \text{T}^{2} \cdot \text{p}_3 = \text{T}^3 \cdot \text{p}_2 = \text{T}^4 \cdot \text{p}_1 = \text{T}^{5} \cdot \text{p}_0}$$
This means, for generalization’s sake,
$$\text{p}_n = \text{T}^{n}\cdot \text{p}_0$$
Using some very basic Linear Algebra methods to compute $\text{T}^n$,

Recall : $\text{T}=\mathcal{u}\cdot D \cdot\mathcal{u}^{-1}$

Here, $D$ is a Diagonal matrix and $\mathcal{u}$ is a matrix whose columns correspond to the Eigen-vectors of $\text{T}$.

I leave the computation to your own practice.

There for,
$$\text{T}^n = \mathcal{u} \cdot D^n . \mathcal{u}^{-1}$$
Since, $D$ is a diagonal matrix, which is obviously of form,
$$D = \begin{bmatrix} a & 0 & 0 & 0 \\ 0 & b & 0 & 0 \\ 0 & 0 & c & 0 \\ 0 & 0 & 0 & d \end{bmatrix}$$
Therefore, $D^n$ will be,
$$D = \begin{bmatrix} a^n & 0 & 0 & 0 \\ 0 & b^n & 0 & 0 \\ 0 & 0 & c^n & 0 \\ 0 & 0 & 0 & d^n \end{bmatrix}$$
Hence, we can easily compute $T^n$.

The Python implementation of this concept is relatively simple once we understand the math behind it.

All the best with that.

Cheers!

Resources:

1. Matrix Diagonalization

2. Eigen Values and Eigen Vectors

3. Markov Chains

4. Markov Chains

Catalan's Conjecture - A learning exercise for the bored mind - Contd.

Before we move to Mihailescu’s proof, let us cover some more mathematical concepts (so that we feel superior to the people around us, just kidding :smile:),

There is an interesting theorem that Mihailescu uses in his proof, called the Stickelberger’s theorem,

---

Stickelberger’s theorem:

This is a result of algebraic number theory, which gives more information about Galois Module structure of class groups of Cyclotomic Fields.
This theorem consists of Stickelberger’s element and Stickelberger’s ideal.
I will now state the complete definition of the theorem and then visit its corners as we move along,

Let $K_m$ denote the $m$-th cyclotomic field${[1]}$ . It is a Galois extension of $\mathbb{Q}$ with Galois Group $G_m$ isomorphic to the multiplicative group of integers modulo $m$ $(\mathbb{Z}/m\mathbb{Z})^{\times}$.

The Stickelberger element (of level $m$ or of $K_m$) is an element in thr group ring $\mathbb{Q}[G_m]$ and the Stickelberger ideal is an ideal in the group ring $\mathbb{Z}[G_m]$. (Note: $\mathbb{Z} \subset \mathbb{Q}$).

The definition of both the Stickelberger element and ideal are, let $\zeta_m$ denote a primitive $m$-th root of unity ${[2]}$. The isomorphism from $(\mathbb{Z/mZ})^{\times}$ to $G_m$ is given by sending $a$ to $\sigma_a$ by the relation $$\sigma_a(\zeta_m)=\zeta^{a}_m.$$
The Stickelberger element of level $m$ is given by,
$$\theta(K_m)=\frac{1}{m}\sum_{a=1 \ (a,m)=1}^ma\cdot \sigma_a^{-1} \in \mathbb{Q}[G_m].$$
The Stickelberger ideal of level $m$ is given by,
$$I(K_m)=\theta(K_m) \mathbb{Z}[G_m] \cap \mathbb{Z}[G_m].$$

---

Inkeri$[3]$, used the concept of Wiefrich pair (explain in the previous blogpost of this series) in the context of Catalan’s equation as follows:

A Wieferich pair is a pair $(p,q)$ of primes such that $$p^{q-1} \equiv 1( \text{mod} \ q^2)$$ and $$q^{p-1} \equiv 1(\text{mod} \ p^2).$$

He showed that if the Catalan’s equation holds, then either $(p,q)$ is a Wieferich pair, or $q$ divides $h_p$, the class number$[4]$ of cyclotomic field $\mathbb{Q}(\zeta_p)$, or $p$ divides $h_q$, the class number of cyclotomic field $\mathbb{Q}(\zeta_q)$, there were other developments in this direction too.

Bugeaud and Hanrot $[5]$proved a class number criterion concerning Catalan’s equation, which implies that the Catalan’s Equation ( $x^p - y^q =1$) has no solutions in non-zero integers $x$ and $y$ if $p$ and $q$ are primes such that one of them is smaller than 43. This was a huge achievement , I recommend you to have a look at the paper at [5].

Mihailescu proved that the Catalan equation has no solutions if $p$ and $q$ are odd and $q$ does not divide $p - 1$. By this result the Catalan conjecture became a theorem. And later Mihailescu succeeded in finding a more elegant proof of Catalan’s conjecture in the case where $q$ does divide $p-1$. Thus, Catalan’s conjecture is a theorem with an algebraic proof in which no computer calculations
are needed.

---

---

In this section e wll discuss some breakthrough results by Mihailescu. The most important one is that $q^2$ divides $x$.

The following lemma will be used for that, an element $a$ of a ring $R$ is called nilpotent if $\exists$ and integer $n$ such that $a^n=0$.

Lemma 1: The ring $\mathbb{O}_k/(q)$ does not contain nilpotent elements, if $\alpha$, $\beta$ $\in$ $\mathbb{O}_k = \mathbb{Z}[\zeta]$ satisfy the congruence $\alpha^p \equiv \beta^q (\text{mod} \ q^2)$

Theorem 1: For $\theta \in {I_s}^{-}$, the element $(x-\zeta)^\theta$ is a $q^{th}$ in $K=\mathbb{Q}(\zeta)$. We also have that $q^2$ divides $x$ and $p^2$ divides $y$.

The proofs of the above Lemma and Theorem are beyond the scope of this blog; regardless to say that the pre-requisites are already covered in detail. For the interested, you can refer Catalan’s Conjecture - A cyclotomic field.

Cheers!

Catalan's Conjecture - A learning exercise for the bored mind.

I was going through a video by Numberphile where they were talking about the Catalan’s conjecture.
This is another conjecture which is simple to state; however extremely difficult to prove; just like Collatz conjecture, about which I already have a blog in place.
At the very outset of the blog; let me bring it to your knowledge that this blog is just for learning and understanding and contains “very little”-to-“null” original work on the subject. However, this blog will be exhaustive and will contain a lot of references that will help a newbie (like myself) get into the depths of the conjecture and fully understand concepts used in its proof.
Let us now define the conjecture;
This statement was conjectured by Eugene Catalan (1814–1894) and was sent to the editor of Journal fur die Reine und Angewandte Mathematik$^{[1]-\text{Introduction}}$.
The conjecture is as follows,
$2^3$ and $3^2$are two powers of natural numbers whose values are consecutive (i.e., 8 and 9); the conjecture is for a mathematical statement such as,
$$\boxed{x^a-y^b=1} \tag*{(1)}$$
The only solution of $(1)$ is for $a$, $b$ $>$ $1$ and $x$, $y$ $>$ $0$ is $x=3$, $a=2$, $y=2$, $b=3$.
Catalan’s conjecture was proven true by Preda Mihăilescu in 2002; the proof involves the theory of cyclotomic fields and Galois Modules$^{[2]-\text{History}}$.
So, as you see now, the breadth of the subject blew up! On second were talking about squares and cubes and now we are talking of Cyclotomic fields and Galois Modules!
Nevertheless, I will try and cover each topic in brief and quench our mathematical thirst!
Let us consider $(1)$ (with $a=p$ and $b=q$; for the sake of eliminating any confusion between $a$ and an english “a”), unless otherwise stated, $x$ and $y$ can be negative integers as well. Now, we re-write $(1)$ as follows,
$$(x-1)\frac{x^p-1}{x-1}=y^q \tag*{(2)}$$
The GCD of the two factors on the left hand side of the equation (after considering $x^p = ((x-1)+1)^p$) is either $p$ or $1$ (How did this happen?).

---

Some concepts before we move ahead:

The Wieferich pairs
$^{[3]-\text{Wieferich pair}}$In mathematics,
a Wieferich pair is a pair of prime number $p$ and $q$ that satisfy,
$$p^{q-1} \equiv 1( \mod p^2)$$
Let us write $(2)$ as, $$(x-1)\frac{x^p-1}{x-1}=pb^q$$ for, $p \nmid b$
This suggests a traditional approach of factorizing the left hand side
in $\mathbb{Z}[\zeta]$, the ring of integers in the $p$th cyclotomic
field $\mathbb{Q}(\zeta)(\zeta=e^{2\pi i/p})$

---

Ring$^{[4]-\text{Ring-Wolfram}}$:
A ring in the mathematical sense is a set $S$ together with two binary
operators $+$ and $.$ (Addition and multiplication), satisfying the
following conditions:
1. Additive associativity: For all $a,b,c \in S$, $(a+b)+c=(a+b)+c$,
2. Additive commutativity: For all $a,b\in S$, $a+b=b+a$,
3. Additive identity: There exists an element $0$ in $S$ such that for all a in $S$, $0+a=a+0=a$.
4. Additive inverse: For every a in $S$ there exists $-a \in S$ such that $a+(-a)=(-a)+a=0$,
5. Left and right distributivity: For all $a,b,c \in S$, $a.(b+c)=(a.b)+(a.c)$ and $(b+c).a=(b.a)+(c.a)$,
6. Multiplicative associativity: For all $a,b,c \in S$, $(a.b).c=a.(b.c)$ ($a$ ring satisfying this property is sometimes
explicitly termed an associative ring). Conditions 1-5 are always
required. Though non-associative rings exist, virtually all texts also
require condition 6.
7. Multiplicative commutativity: For all $a,b \in S$, $a.b=b.a$ ($a$ ring satisfying this property is termed a commutative ring),
8. Multiplicative identity: There exists an element $1 \in S$ such that for all $a \neq 0 \in S$, $1.a=a.1=a$ (a ring satisfying this
property is termed a unit ring, or sometimes a “ring with identity”),
9. Multiplicative inverse: For each $a \neq 0 \in S$, there exists an element $a^{-1} \in S$ such that $\forall a \neq 0 \in S, a.a^{-1}=a^{-1}.a=1$, where $1$ is the identity element.

---

A brief history of the past developments on this conjecture is a must to be read and understood; the significance comes due to the period of 150 years for which it remained an open problem$^{[5]-\text{Catalan's Conjecture: Are 8 and 9 the Only Consecutive Powers?}}$,

• Only after six years after Catalan formally defined the conjecture, a result was proposed by French mathematician, Victor Lebesgue. He stated that, for the equation, $x^p-y^2=1$; where $p$ is a prime; has no solutions for positive values of $x$ and $y$. A proof of the same will be discussed in brief in the later part of the blog.
• After Lebesgue’s work, all development solely consisted of small exponents, and then Naggel showed in 1921 that the difference between a third power and an other perfect power never is equal to 1.
• In 1932, Selberg proved that, $x^4-y^n=1$ has no solution in positive integers when $n>1$. A stronger result to this was proved by Ko Cho in 1965, that stated that the equation $x^2-y^q=1$ has no solutions for positive integers when $q \geq 5$.
• Cassels made some observations for $x^p-y^q=1$ where $p$ and $q$ are odd-primes. He proved that is this equality holds for positive integers $x$ and $y$, then $p$ divides $y$ and $q$ divides $x$. For the case $p=2$, this had already been shown by Naggel
• Inkeri defined the concept of a Wieferich pair [the definition and explanation of the same is given above] in the concept of Catalan equation as follows:
  If the Catalan’s equation $(1)$ holds, then either $(p,q)$ is a Wieferich Pair, or $q$ divides $h_p$, the class number of the cyclotomic field $\mathbb{Q}(\zeta_p)$, or $p$ divides $h_q$, the class number of the cyclotomic field $\mathbb{Q}(\zeta_q)$

---

Cyclotomic Field$^{[5]-\text{Cyclotomic field}}$:
In number theory, a cyclotomic field is a number field obtained by adjoining a complex primitive root of unity to $\mathbb{Q}$, the field of rational numbers. The $n$-th cyclotomic field $\mathbb{Q}(\zeta_n)$ (where $n > 2$) is obtained by adjoining a primitive $n$-th root of
Primitive root of unity$^{[6]-\text{Root of unity}}$
In mathematics, a root of unity, occasionally called a de Moivre number, is any complex number that gives $1$ when raised to some positive integer power $n$.

---

• Some time later Mihailescu proved that the Catalan equation has no solutions if $p$ and $q$ are odd and $q$ does not divide $p-1$. By this result the Catalan conjecture became a theorem

---

---

More on Cyclotomic Fields:
Let $p$ be an odd-prime number. Let $\Phi_p$ be the $p$-th cyclotomic polynomial in $\mathbb{Q}[X]$ i.e., $\Phi_p = \frac{X^p-1}{X-1}$. Consider the field extension $\mathbb{Q}[X]/(\Phi_p) \cong \mathbb{Q}(\zeta)$ of $\mathbb{Q}$, where $\zeta$ denotes a primitive $p$-th root of unity. This is a field extension of degree $p-1$ and it is reducible in $\mathbb{Q}[X]$. We denote $\mathbb{Q}(\zeta)$ by $k$ from now on.
This field extension is Galois with Galois Group,
$$G=\text{Gal}(\mathbb{Q}(\zeta)/\mathbb{Q}) \cong (\mathbb{Z}/p\mathbb{Z})^{*},$$
Since the map,
$$(\mathbb{Z}/p\mathbb{Z})^{*} {{\sim} \atop {\to}}\text{Gal}(\mathbb{Q}(\zeta)/\mathbb{Q})$$
$$a (\mod p) \mapsto (\sigma_a : \zeta \mapsto \zeta^a)$$
is an isomorphism
The automorphism $\sigma_{p-1}$ acts in all embeddings as complex conjugation. Therefore, we call $\sigma_{p-1}$ complex conjugation.
The fixed field of complex conjugation is $\mathbb{Q}(\zeta+\zeta^{-1})$, which is called the maximal real subfield of $\mathbb{Q}(\zeta)$. We denote $\mathbb{Q}(\zeta + \zeta^{-1})$ by $K^{+}$. The field extension $\mathbb{Q}(\zeta+\zeta^{-1})$ of $\mathbb{Q}$ has degree $\frac{p-1}{2}$ and it is Galois with Galois theory.
$$G^{+}=\text{Gal}(\mathbb{Q}(\zeta+\zeta^{-1})/\mathbb{Q}) \cong (\mathbb{Z}/p\mathbb{Z})^{*}/(\pm 1)$$

Another important concept that is

Mihailescu’s proof

To be contd…in the next blogpost!

i!

Define the Factorial of a Complex number.

In usual sense, factorial is defined as,
$$n! = \prod_{k=1}^{n} k=1 \cdot2\cdot3\ldots(n-2)\cdot(n-1)\cdot n \tag*{(1)}$$
Now, the not so usual definition is based on the famous Gamma Function,
$$\Gamma(t) = \int_0^\infty e^{-x}x^{t-1}\mathrm{d}x \tag*{(2)}$$
There is an unique and very useful property,
$$\Gamma(n+1)=n! \tag*{(3)}$$
To extend $(3)$ into the complex domain, we will first have to go through Analytic Continuation, please read about Analytic Continuation here.
Therefore, after analytic continuation, we can write it as,
$$z! \overset{\text{def}}{=} \Gamma(z+1) \tag*{(4)}$$
For, $z \in \mathcal C-{-1,-2...}$
So, now, clearly,
$$i!= \Gamma(i+1) \tag*{(5)}$$
By $(3)$
$$\Gamma(i+1)=\int_0^\infty e^{-x} x^{(i+1)-1}\mathrm{d}x \tag*{(6)}$$
Clearly,
$$\Gamma(i+1)=\int_0^\infty e^{-x} x^{i}\mathrm{d}x \tag*{(6)}$$
For easier computation, please catch that,
$$x^i = e^{i \ln x} = \cos ( \ln x)+i \sin( \ln x)$$
Let’s break it down,
$$i! = \int_0^\infty e^{-x} (\cos ( \ln x)+i \sin( \ln x))\mathrm{d}x$$
$$\implies i! = \int_0^\infty e^{-x} \cos ( \ln x)\mathrm{d}x+i\int_0^\infty e^{-x} \sin ( \ln x)\mathrm{d}x$$
If you have reached this far, you obviously know how to solve the above integral.
Cheers!

Study of Brachistochrone Problem.

This blog will deal with one of the most elegant topics in mathematics. The famous Brachistochrone curve. Vsuace did a great job at explaining it in this video. This looks good to me, but it doesn’t feel good; because of course the video does not contain tons of fancy mathematical vocabulary and millions of lines of $\LaTeX$.

To get started, we will just put forward a small fact that, this topic is from a larger field of mathematics called the Calculus of Variations.

Before we get into the Analysis of Brachistochrone Curve (which happens to be the heart of this post), let us brush-up on some concepts.

---

Partial Differentiation (advanced preliminary):

Let us look at how we define the derivative of a given function, say we have $g(x)$ and we have,
$$\boxed{\frac{\mathrm{d}}{\mathrm{d}x}g(x)=g'(x)=\lim_{c \to 0}\frac{f(x+ c)-f(x)}{c}}$$
Provided this limit exists.

Now, for a function of many variables it is not easy to compute the total derivative (the usual derivative). Therefore, we calculate the Partial Derivative. You might have already guessed it by now, we define partial derivative as, say we have a function $h(x,y,z)$, then,
$$\boxed{\frac {\partial}{\partial x}h(x,y,z) = \lim_{c \to 0} \frac{h(x+c, y,z)-f(x,y,z)}{c}}$$
Let us look into a small example now, say we have,
$$f(x,y,z) = x^2y^3+x+y+z+\sin(x) \cos(y)$$
Now,
$$\frac{\partial}{\partial x}f(x,y,z) = 2xy^3+1+0+0+\cos(x) \cos(y)$$
similarly,
$$\frac{\partial}{\partial y}f(x,y,z) = 3x^2y^2+0+1+0-\sin(x) \sin(y)$$
and,
$$\frac{\partial}{\partial z}f(x,y,z) = 0 +0+0+1+0$$

---

Idea of Speed and Time (basic preliminary):

Let us consider, the following scenario, (Made using Geogebra)

We have, in the figure we have denoted the path ACB as $\pi_1$ and the path ADB as $\pi_2$.
Consider that, The time taken for a body to move from A to B using path $\pi_1$ is $t_{AB_{\pi_1}}$ similarly, the time taken for a body to move from A to B using path $\pi_2$ is $t_{AB_{\pi_2}}$.
Also, consider that the distance for $\pi_1$ be $l_1$ and the distance for $\pi_2$ be $l_2$.

Then,
Speed ($\pi_1$) =$v_1$= $\frac{l_1}{t_{AB_{\pi_1}}}$ and Speed ($\pi_2$) =$v_2$ =$\frac{l_2}{t_{AB_{\pi_2}}}$.

3. Idea of Distance(basic preliminary):

Distance between two points $(a,b)$ and $(c,d)$ is given by,
$$l = \sqrt{(a-c)^2+(b-d)^2}$$
In case, the distance is measured from origin ($O(0,0)$) to some point $A(x,y)$, we get,
$$l = \sqrt{(x-0)^2+(y-0)^2} = \sqrt{x^2 + y^2}$$

---

Now that the notion is clear, we shall proceed into the analysis.
[NOTE : Some topics which are new to me (or you) will be explained as and when it is required].

---

Time needed for a body to travel from $A$ to $B$ is given by,
$t_{AB}=\frac{l}{v}$ for a linear path.

For a curve or a non-linear path, we consider piece-wise linear distance ($\mathrm{d}l$) and speed ($v$). We can define, $t_{AB}$ as,
$$\boxed{t_{AB} = \int_{A}^{B}\frac{\mathrm{d}l}{v}} \tag*{(1)}$$
Now, we must understand a fact that, all distance is composed of $x$ and $y$ coordinates. We consider the point $A$ as the origin ($O(0,0)$) and $B$ as $(x,y)$. Hence,
$$l = \sqrt{x^2 + y^2}\\ \implies \mathrm{d}l = \sqrt{\mathrm{d}x^2 + \mathrm{d}y^2}\\ \implies \mathrm{d}l = \sqrt{(1 + \frac{\mathrm{d}y^2}{\mathrm{d}x^2})\mathrm{d}x^2}$$
$$\boxed{\implies \mathrm{d}l = \sqrt{(1 + \frac{\mathrm{d}y^2}{\mathrm{d}x^2})}\mathrm{d}x} \tag*{(2)}$$

---

Also, since the translation is both in the $x$ and $y$ axes, we can say that the speed is gained by $KE$ and $PE$ equality,
$$KE = PE\\ \implies \frac{1}{2}mv^2 = mgy \\ $$
$$\implies \boxed{v = \sqrt{2gy}} \tag*{(3)}$$

---

Use $(2)$ and $(3)$ in $(1)$, we have,
$$\boxed{t_{AB} = \int_{A}^{B}\frac{\sqrt{(1 + \frac{\mathrm{d}y^2}{\mathrm{d}x^2})}\mathrm{d}x}{\sqrt{2gy}}} \tag*{(4)}$$
We write
$$\frac{\mathrm{d}y}{\mathrm{d}x}=y'$$
Now, $(4)$ becomes,
$$t_{AB} = \int_{A}^{B}\frac{\sqrt{(1 +y'^2 )}}{\sqrt{2gy}}\mathrm{d}x \tag*{(4.1)}$$
Therefore, the function is,
$$\theta = \frac{\sqrt{(1 +y'^2 )}}{\sqrt{2gy}}\\ \boxed{\theta = (1 +y'^2 )^{1/2}(2gy)^{-1/2}} \tag*{(5)}$$

---

Please notice that, the integral $t_{AB}$ can be written as,
$$t_{AB} = \int_{A}^{B} f(y, y', x) \mathrm{d}x \tag*{(6)}$$
In our case, $f(y, y', x)$ is $\theta$.

Equation$(6)$ closely follows the The Euler-Lagrange equation.

It is pretty simple to define,

---

The Euler-Lagrange equation

For any $I$ such that,
$$I = \int_{a}^{b}f(y,y',x)\mathrm{d}x$$
Where, $y' = \frac{\mathrm{d}y}{\mathrm{d}x}$, then $I$ has an stationary point, if the Euler-Lagrange Equation, given by,
$$\frac{\partial f}{\partial y}-\frac{\mathrm{d}}{\mathrm{d}x}\frac{\partial f}{\partial y'}=0$$
is satisfied.

Therefore, for our analysis,
$$\boxed{\frac{\partial \theta}{\partial y} - \frac{\mathrm{d}}{\mathrm{d}x}\frac{\partial \theta}{\partial y'}=0} \tag*{(7)}$$

---

Stationary value

This is a value at the stationary point.
A stationary point $a$ is the point at which the first derivative of a function becomes zero,
$$f'(x)|_{x=a} = 0 ; a = \text{Stationary point}$$
Side note :
Find the stationary points of $f(x) = x^3-6x^2+9x$.

---

Beltrami Identity

Our $\theta$ (from $(5)$) is such a cool expression, because $x$ does not appear in that; therefore, of course $\boxed{\frac{\partial \theta}{\partial x} = 0}$; this again leads us to another beautiful form, Beltrami Identity, which is given by,
$$\boxed{\theta - y' \frac{\partial \theta}{\partial y'}=C ; C \text{ is a constant}}.\tag*{(8)}$$

---

We now find, $\frac{\partial \theta}{\partial y'}$,
$$\frac{\partial \theta}{\partial y'}=\frac{\partial}{\partial y'}[(1 +y'^2 )^{1/2}(2gy)^{-1/2}] \\ \implies \frac{\partial \theta}{\partial y'} =(2gy)^{-1/2} \frac{1}{2}(1 +y'^2 )^{-1/2}\cdot 2y'$$
$$\therefore \boxed{\frac{\partial \theta}{\partial y'} = (2gy)^{-1/2} (1 +y'^2 )^{-1/2}\cdot y'} \tag*{(9)}$$

---

Use $\theta$ and $\frac{\partial \theta}{\partial y'}$ in equation $(8)$, we get,
$$(1 +y'^2 )^{1/2}(2gy)^{-1/2}-y' \cdot (2gy)^{-1/2} (1 +y'^2 )^{-1/2}\cdot y' = C$$
$$\implies \frac{\sqrt{(1 +y'^2 )}}{\sqrt{2gy}}-\frac{y'^2}{\sqrt{2gy}\sqrt{1+y'^2}} = C$$
$$\implies \frac{1+y'^2-y'^2}{\sqrt{2gy}\sqrt{1+y'^2}}=C$$
$$\implies \boxed{\frac{1}{\sqrt{2gy}\sqrt{1+y'^2}}=C} \tag*{(10)}$$
Rearranging a little gives,
$$\sqrt{1+y'^2}=\frac{1}{\sqrt{2gy}C}$$
$$\implies \boxed{(1+\frac{\mathrm{d}y}{\mathrm{d}x})y = \frac{1}{2gC^2} = M^2} \tag*{(11)}$$
This is the equation of the cycloid as per [4], I am yet to figure out how this happened.
The solution of $(11)$ can be found using a parametric equation,
$$\boxed{x(t) = \frac{1}{2}M^2(t-\sin t)\\ y(t)=\frac{1}{2}M^2(1-\cos t)} \tag*{(12)}$$
To derive the above equation, please refer Math-stackexchange|Derive the parametric equations of a cycloid.

I have even plotted the solutions at $(12)$, we can see that it is a cycloid,

If you wish to play with the visualisation please go to Cycloid | Parametric @ Desmos by Pragyaditya

Cheers!

---

References:
1. Stationary Points
2. Euler Lagrange Equation
3. Derivation of Beltrami Identity
4. Introduction to calculus of variations
5. Brachistochrone @ Wolfram

A Nice Problem Involving Gamma and Beta Function.

[Playing with Gamma and Beta functions].
Solve :
$$\int_{0}^{\pi/2}(\cos^{2}x-\cos^{3}x) \mathrm{d}x$$

Let,
$$I = \int_{0}^{\pi/2}(\cos^{2}x-\cos^{3}x) \mathrm{d}x\tag*{(1)}$$

Beta function is defined as:
$$\boxed{B(x,y) = \int_{0}^{1}t^{x-1}(1-t)^{(y-1)}\mathrm{d}t \hspace{0.2cm}\forall x, y >0}$$

Now, there are many properties of this beautiful mathematical tool, some of them are,

1. $$B(x,y) = B(y,x)$$ They are symmetric.
2. We can write Beta Functions in terms of Gamma Functions as follows, $$B(x,y) = \frac{\Gamma(x) \cdot \Gamma(y)}{\Gamma(x+y)}.$$
3. We can obviously add factorial into the mix too, $$B(x,y)=\frac{(x-1)!(y-1)!}{(x+y-1)}$$
4. We can define it in terms of trigonometric angles too,
   $$B(x,y)=2 \int_{0}^{\pi/2}(\sin \theta)^{2x-1}(\cos \theta)^{2y-1}\mathrm{d}\theta$$

Now,
$$I = \int_{0}^{\infty}(\cos^{2}x)\mathrm{d}x-\int_{0}^{\infty}(\cos^{3}x) \mathrm{d}x \\ \implies I = I_1 + I_2$$

We will solve both separately,
$$I_1 = \int_{0}^{\infty}(\cos^{2}x)\mathrm{d}x \\ \implies I_1 = \int_{0}^{\infty}(\sin x)^0(\cos x)^2\mathrm{d}x \\ \implies I_1 = \int_{0}^{\infty} (\sin x)^{2.\frac{1}{2}-1}(\cos x)^{2.\frac{3}{2}-1}\mathrm{d}x \\ \implies I_1 = \frac{B(\frac{1}{2},\frac{3}{2})}{2}\\ \implies \boxed{ I_1 = \frac{\Gamma(\frac{1}{2}). \Gamma(\frac{3}{2})}{2\Gamma(\frac{1}{2}+\frac{3}{2})}}$$
Similarly, for the other integral,
$$I_2 = \int_{0}^{\infty}(\cos^{3}x)\mathrm{d}x\\ \implies I_2 = \int_{0}^{\infty}(\sin x)^0 (\cos x)^3 \mathrm{d}x\\ \implies I_2 = \int_{0}^{\infty}(\sin x)^{2.\frac{1}{2}-1}(\cos x)^{2.2-1}\mathrm{d}x\\ \implies I_2 = \frac{B(\frac{1}{2},2)}{2}\\ \implies \boxed{ I_2 = \frac{\Gamma(\frac{1}{2}). \Gamma(2)}{2\Gamma(\frac{1}{2}+2)}}$$
Therefore,
$$\boxed{\int_{0}^{\infty}(\cos^{2}x-\cos^{3}x) \mathrm{d}x = \frac{\Gamma(\frac{1}{2}). \Gamma(\frac{3}{2})}{2\Gamma(\frac{1}{2}+\frac{3}{2})}+\frac{\Gamma(\frac{1}{2}). \Gamma(2)}{2\Gamma(\frac{1}{2}+2)}}$$
Cheers $$\bowtie$$

Integration : Problem 2

$$\text{I have been inactive for quite a long time.}\\ \text{As usual, I will start this blogging streak}\\ \text{with a small warm-up exercise.}\\ \text{I will solve some mildly complex Integrals, }$$

$\underline{\text{Question 1.}} \\$ $\LARGE \int _{0}^{\frac{1}{2}}\frac{x^2}{(1+x^2)^{3/2}} \mathrm{d}x$
$\underline{\text{Solution}}$:
$$\text{Let, } I_{1} = \int _{0}^{\frac{1}{2}}\frac{x^2}{(1+x^2)^{3/2}} \mathrm{d}x \\ \text{We can write, } I_1 = \int _{0}^{\frac{1}{2}} \frac{x^2 + 1 - 1}{(1+x^2)^{3/2}} \mathrm{d}x \\ \implies I_1 = \int _{0}^{\frac{1}{2}} \frac{x^2+1}{(1+x^2)^{3/2}}\mathrm{d}x - \int_{0}^{\frac{1}{2}}\frac{1}{(1+x^2)^{3/2}}\mathrm{d}x\\ \text{Which with some algebraic manipulation comes to, }\\ I_1 = \int_{0}^{\frac{1}{2}} \frac{1}{\sqrt{1+x^2}}\mathrm{d}x - \int_{0}^{\frac{1}{2}}\frac{1}{(1+x^2)^{3/2}}\mathrm{d}x\\ \text{Let,} I_2 = \int \frac{1}{\sqrt{1+x^2}}\mathrm{d}x \\ \text{and, Let , } I_3 = \int \frac{1}{(1+x^2)^{3/2}}\mathrm{d}x\\ \text{We know that, } \boxed{\boxed{\int \frac{\mathrm{d}x}{\sqrt{x^2 + a^2}} = \log|x +\sqrt{x^2 + a^2}|+c }} \\ \therefore I_2 = \int \frac{1}{\sqrt{1+x^2}}\mathrm{d}x = \log|x+\sqrt{1+x^2}| \\ \text{ neglecting } c \text{ because it will ultimately} \\ \text{ be a definite integral }.\\ \\ \text{Now the rigor starts, } I_3 = \int \frac{1}{(1+x^2)^{3/2}}\mathrm{d}x\\ \text{Take, } x = \tan {t} \rightarrow \frac{\mathrm{d}x}{\mathrm{d}t} = \sec^2 {t} \\ \therefore I_3 = \int \frac{\sec^2{t}}{(1+\tan^2{t})^{3/2}} \mathrm{d}t\\ \text{We know, } \boxed{\boxed{1 + \tan^2 \theta = \sec^2 \theta}} \\ \text{Using this } I_3 \text{ changes to, }\\ I_3 = \int \frac{1}{\sec t} \mathrm{d} t \implies \int \cos {t} \mathrm{d}t = \sin t \\ \text{We get the substitution in place ,}\\ t = \tan^{-1} x \rightarrow \sin t \rightarrow \sin{(\tan^{-1} x)} = \frac{x}{\sqrt{x^2+1}} \\ \therefore I_3 = \sin t = \frac{x}{\sqrt{x^2+1}} \text{Now,} \\ I_1 = I_2 - I_3 \text{ with no limits. }\\ \therefore I_1 = \int \frac{1}{\sqrt{1+x^2}}\mathrm{d}x = \log|x+\sqrt{1+x^2}| - \frac{x}{\sqrt{x^2+1}}\\ \text{Put limits } \frac{1}{2} \text{ to } 0 \text{ and get answer as, }\\ I_1 = {I_1}|_{0}^{1/2}=\boxed{0.033}.\\ \text{Cheers! } \bowtie$$

Life of \pi

Is $\pi = \frac{22}{7}$?

There was an interesting problem that had come up during one of the JEE exams.

It asked us to integrate,
$$I = \int_{0}^{1} \frac{x^4 (1-x)^4}{1+x^2} dx$$
Use Polynomial Long Division,
$$I=={\displaystyle\int}\left(-\dfrac{4}{x^2+1}+x^6-4x^5+5x^4-4x^2+4\right)\mathrm{d}x$$

Use Linearity,
$$=-\class{steps-node}{\cssId{steps-node-1}{4}}{\displaystyle\int}\dfrac{1}{x^2+1}\,\mathrm{d}x+$$

$${\displaystyle\int}x^6\,\mathrm{d}x-\class{steps-node}{\cssId{steps-node-2}{4}}{\displaystyle\int}x^5\,\mathrm{d}x+\class{steps-node}{\cssId{steps-node-3}{5}}{\displaystyle\int}x^4\,\mathrm{d}x-\class{steps-node}{\cssId{steps-node-4}{4}}{\displaystyle\int}x^2\,\mathrm{d}x+\class{steps-node}{\cssId{steps-node-5}{4}}{\displaystyle\int}1\,\mathrm{d}x$$

$$I=-4\arctan\left(x\right)+\dfrac{x^7}{7}-\dfrac{2x^6}{3}+x^5-\dfrac{4x^3}{3}+4x$$

Substitute the values and we will get,
$$I = -\dfrac{7{\pi}-22}{7} = -(\pi - \frac{22}{7})$$
Which comes to,
$$\displaystyle I = \frac{22}{7} - \pi \tag*{(1)}$$

Since, Area cannot be negative, it is evident from $(1)$ that,
$$\frac{22}{7} \geq \pi$$


Also, $\pi \neq \frac{22}{7}$ is evident from,
$$\boxed{ e^{\pi i}+1 \neq e^{\frac{22}{7} i}+1}$$
Cheers!

Preparation of exam.

Preparation Notes - Hell lot of Math and Algorithms; very helpful.

Note : These notes are made during my preparation for my exams. They might be crude and may not be in their best form of presentation. Read at your own risk.

---

Subject : IC 021 - Computational Techniques in Control Engineering

Professor : Dr. A Ramakalyan.

Difficulty : Runs to Heavens and goes around it $\infty$ times.

Books:
1. B N Datta’s Numerical Methods for Linear Control Systems. - referenced as [1]
2. Gilbert Strang’s Linear ALgebra and it’s applications. - referenced as [2]
3. S Kumaresen’s Linear Algebra - A Geometric Approach. - referenced as [3]

---

1. Computing $e^\mathcal{A}$[1]

1. There are many ways of doing this, one of the very popular one is the Taylor Series method.
   We simply expand $e^\mathcal{A}$ as per the Taylor series.
   Hence, $e^\mathcal{A}$ can be written as,
   $$\boxed{e^\mathcal{A} = I + \mathcal{A} + \frac{\mathcal{A}^{2}}{2!} + \frac{\mathcal{A^3}}{3!} + ...}$$

   • The drawbacks to this method are that a large number of terms is needed for convergence, and even when convergence occurs, the answer can be totally wrong.

2. Hence, I will now get into the Pade’s Scaling and Squaring algorithm for this task.
   Suppose that, $f(x)$ is a power series presented as,
   $$f(x) = f_0 + f_1 x+ f_2 x^2+...$$
   Then the Pade’s Scaling and Squaring Algorithm, given an approximation of $f(x)$ as ,
   $$\boxed{f(x) \approx \frac{c(x)}{d(x)}= \frac{\sum_{k=0}^p c_{k} x^{k}}{\sum_{k=0}^q d_{k} x^{k}}}$$
   The $c_k$ and $d_k$ are chosen such that the terms containing $x^0, x^1, x^2,...,x^{p+q}$ are cancelled out in $f(x)d(x) - c(x)$. The order of Pade’s approximation is $p+q$. The ratio $(\frac{c(x)}{d(x)})$ is unique and the coefficients $c_0, c_1,...,c_p$ and $d_0, d_1,...,d_q$ always exist.
   
   Similarly, Pade’s approximation for $e^\mathcal{A}$ is given by,
   $$\boxed{R_{pq}(\mathcal{A}) = [D_{pq}(\mathcal{A})]^{-1}N_{pq}(\mathcal{A})}$$
   Where,
   $$\boxed{D_{pq}(\mathcal{A})=\sum_{k=0}^q\frac{(p+q-j)!q!}{(p+q)!j!(q-j)!}(-\mathcal{A})^{j}}$$
   and
   $$\boxed{N_{pq}(\mathcal{A})=\sum_{k=0}^p\frac{(p+q-j)!p!}{(p+q)!j!(p-j)!}(\mathcal{A})^{j}}$$

   • It can be shown (Exercise 5.16) that if p and q are large, or if A is a matrix having all its eigenvalues with negative real parts, then $D_{pq}(\mathcal{A})$ is nonsingular.
   • The difficulty of Round off errors in Pade’s method can be controlled by Scaling and squaring.
   • Since $e^\mathcal{A}$ is same as $(e^\frac{\mathcal{A}}{m})^m$, the idea is it to find out $m$ in terms of powers of $2$ so that $e^\frac{\mathcal{A}}{m}$ can be computed accurately and then find $(e^\frac{\mathcal{A}}{m})^m$ with repetitive squaring.
   • This method has favorable numerical properties when $p=q$.
     
     Let $m = 2^j$ be chosen such that $||A||_{\infty} \leq 2^{j-1}$, then Moler and Van Loan(1978) have shown that there exists an $E$ such that ,
     $$\boxed{[R_{pp}(A/2^j)]^{2^j} = e^{\mathcal{A} + E}},$$
     Where $||E||_{\infty} \leq \epsilon ||\mathcal{A}||_{\infty}$, with
     $$\epsilon = 2^{3-2p} (\frac{(p!)^2}{(2p)!(2p+1)!)}) \\ \boxed{\epsilon = \frac{8}{2^{2p}} (\frac{(p!)^2}{(2p)!(2p+1)!)})}$$
     Given an error tolerance of $\delta$ we should choose $p$ such that $||E||_{\infty} \leq \delta ||\mathcal{A}||_{\infty}$
     
     The Algorithm is as follows:
     Input : $A \in \mathcal{R}^{n \times n}$, $\delta > 0$
     Output : $F = e^{\mathcal{A} + E}$ with $||E||_{\infty} \leq \delta ||\mathcal{A}||_{\infty}$
     
     Step 1: Choose $j$ such that $||\mathcal{A}||_{\infty} \leq 2^{j-1}$, set $\mathcal{A} \equiv \mathcal{A}/2^j$.
     Step 2: Choose $p>0$ and smallest non-negative integer satisfying,
     $$\boxed{(\frac{8}{2^{2p}}) (\frac{(p!)^2}{(2p)!(2p+1)!)}) \leq \delta}.$$
     Step 3: Set $\boxed{D \equiv I, N \equiv I, Y \equiv I, c = 1}$
     Step 4: For $k = 1, 2, ..., p$ do
     $$\boxed{c \equiv c(p - k + 1)/[(2p - k + 1)k] \\ Y \equiv AY,\\ N \equiv N + cY , \\D = D+(-1)^{k}cY.}$$
     Step 5: Solve for $F$ : $\boxed{DF = N}$.
     Step 6: For $k = 1, 2, ..., j$ do,
     $$\boxed{F \equiv F^{2}}$$
   • The algorithm requires about $2(p + j + (1/3)n^3$ flops.
   • The algorithm is numerically stable when A is normal. When A is non-normal, an analysis of the stability property becomes difficult, because, in this case e a may grow before it decays during the squaring process; which is known as “hump” phenomenon.
     Example at Page 143 of [1]

3. We have another effective method called the Matrix decomposition method uses the Real Schur form to do so.
   
   Let $\mathcal{A}$ be transformed to a real Schur matric $R$ using an orthogonal similarity transformation:
   $$\boxed{P^{T}\mathcal{A}P = R}$$
   Then,
   $$\boxed{e^{\mathcal{A}}=Pe^{R}P^{T}}$$
   Note that $R$ is upper-triangular.
   
   The algorithm is:
   Step 1: Transform $\mathcal{A}$ to $R$ using the $QR$ iteration algorithm.
   $$\boxed{P^{T}\mathcal{A}P = R}$$
   Step 2: Compute $e^R = G = g_{ij}$:
   for $i = 1,...,n$ do,
   $$g_{ii} = e^{r_{ii}}$$
   end.
   
   for $k = 1, 2, ..., n-1$ do
   for $i = 1, 2, ..., n-k$ do
   set $j = i+k$
   $$\boxed{g_{ij}=\frac{1}{(r_{ii}-r_{jj})}\Big[r_{ij}(g_{ii}-g_{jj})+\sum_{p=i+1}^{j-1}(g_{ip}r_{pj}-r_{ip}g_{pj})\Big]}$$
   end
   end
   Step 3: Now we have $e^R$ then we compute $e^\mathcal{A}$

---

2.Lyapunov Stability theorem: [1]

Lyapunov Stability theorem,
The system :
$$\boxed{\dot{x}(t)=\mathcal{A}x(t)}$$
is asymptotically stable $iff$, for any symmetric positive definite matrix $M$, there exists a unique symmetric positive definite matrix $X$ satisfying the equation:
$$\boxed{X\mathcal{A}+\mathcal{A}^{T}X=-M}$$

Kronecker Product and Lyapunov:
$$\boxed{\mathcal{A} \bigotimes B = \big[ a_{ij} [B] \big]}$$
For example,
if
$$\mathcal{A} = \begin{bmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \\ \end{bmatrix}$$
and,
$$B = \begin{bmatrix} b_{11} &b_{12} &b_{13}\\ b_{21} &b_{22} &b_{23}\\ \end{bmatrix}$$
Then, the Kronecker product will look like,
$$\mathcal{A} \bigotimes B = \begin{bmatrix} a_{11}b_{11} &a_{11}b_{12} &a_{11}b_{13} &a_{12}b_{11} &a_{12}b_{12} &a_{12}b_{13} \\ a_{21}b_{21} &a_{21}b_{22} &a_{21}b_{23}&a_{22}b_{11} &a_{22}b_{12} &a_{22}b_{13} \\ \end{bmatrix}$$
It is evident that,
$$\boxed{\mathcal{A} \bigotimes B \neq B\bigotimes\mathcal{A}}$$
Now, we formulate the $X\mathcal{A}+\mathcal{A}^{T}X=M$ using the Kronecker product. Note that $M$ is negative semi-definite. $X$ is positive definite and symmetric.

Vectorisation of a matrix $\large \triangleq$ Stack all elements of
the matrix into one column.

Going directly to the final answer,
$$\boxed{X\mathcal{A}+\mathcal{A}^{T}X=M = [(A^{T} \bigotimes I )+( I\bigotimes A^{T})]\overrightarrow{X} = \overrightarrow{M}}$$
Verify it once.

If $\mathcal{A}$ is stable then $X$ has unique solution.


Schur method of Lyapunov:
Method proposed by Bartels and Stewart. It is based on the reduction of $A^{T}$ to $RSF$.

Step 1: Reduce the problem.
For,
$$X\mathcal{A}+\mathcal{A}^{T}X = C$$
is reduced to,
$$\boxed{YR^{T}+RY=\hat{C}}$$
Where,
$R = U^{T}A^{T}U$, $\hat{C}=U^{T}CU$ and $Y=U^{T}XU$.

Step 2: Solution of the reduced form.
Let,
$$Y = (y_1, y_2,...,y_n),$$
$$\hat{C}=(c_1,...,c_n)$$
and $$R=(r_{ij})$$
Assume that the columns $y_{k+1}$ through $y_{n}$ have been computed, and consider the following two cases,

Case 1: $r_{kk-1} = 0$. Then $y_k$ is determined by solving the quasi-triangular system,
$$\boxed{(R+r_{kk}I)y_k=c_k-\sum_{j=k+1}^nr_{kj}y_j}$$
If, in particular, R is upper triangular, that is there are no “Schur bumps” on the diagonal, then each $y_i$, $i=n, n - 1 ..... 2, 1$ can be obtained
by solving an $n \times n$ upper triangular system as follows:
$$\boxed{(R+r_{nn}I)y_n=c_n} \\ \boxed{(R+r_{n-1,n-1}I)y_{n-1}=c_{n-1}-r_{n-1,n}y_n}$$

Case 2: $r_{kk-1} \neq 0$ this shows a Schur Bump.
$$\boxed{R(y_{k-1}, y_k)+(y_{k-1}, y_k) \begin{bmatrix} r_{k-1,k-1} &r_{k,k-1}\\ r_{k-1,k} &r_{k,k}\\ \end{bmatrix}\\ =(c_{k-1},c_k)-\sum_{j=k+1}^n(r_{k-1,j}y_j, r_{kj}y_j)=(d_{k-1},d_k)}$$
Example at Page 274


If $\mathcal{A}$ doesn’t have Eigen Values in the Limit circle , we modify $\mathcal{A}$ as,
$$A_{1} = (I-A)^{-1}(I+A) \text{ and } Q_{1} = (I+A^{T}_{1})(I+A_{1})\\ \text { solve for A } \implies (I-A)A_{1}=(I+A)\\ \implies A_{1}-AA_{1}=I+A\\ \implies A_{1} - I = AA_{1}+A\\ \implies A(A_{1}+I) = A_{1}-I\\ \implies A=(A_{1}-I)(A_{1}+I)^{-1}$$
Use this $A$ in the Discrete Lyapunov which is given by,
$$\mathcal{L}_{D} : A^{T}PA-P=Q$$
The resulting equation is,
$$A^{T}_1P+PA_{1}=Q_{1} \text{ assume, } Q = -2I$$
You will get a similar equation to that of the normal Lyapunov equation.

---

3. Modified QR algorithms:[4]

1. Modified Gram-Schmidt method:
   for
   $$k = 1:n\\ q^{T}_k = a^{T}_k$$
   end
   
   for
   $$k = 1:n\\ r_{kk} = (q^{T}_{k}q_k)^{1/2}\\ q^{T}_k= \frac{q^{T}_k}{r_{kk}}$$
   -for
   $$j = k+1:n\\ r_{kj}=a^{T}_{j}q_{k}\\ q^{T}_{j}=q^{T}_{j}-r_{kj}q^{T}_{k}$$
   -end
   end
   It requires $\boxed{O(\frac{8}{3}n^{3})}$ operations.
   Practice example from the above link

2. Original algorithm of QR decomposition:
   Aim is to write matrix $\mathcal{A}$ in terms of $\boxed{Q\times R}$.
   Here, $Q$ is an orthogonal matrix with orthonormal columns; and $R$ is an upper-triangular matrix.
   
   Now, we take an example,
   $$\mathcal{A}=\begin{bmatrix} 1 & 2 & 3\\ -1 & 1 & 1\\ 1 & 1 & 1\\ 1 & 1 & 1 \end{bmatrix}$$
   We will first find $q_1$, $q_2$ and $q_3$ as the three orthonormal vectors using the GS-strategy.
   Now, we will define $a$, $b$ and $c$ as,
   $$a = \begin{bmatrix} 1\\ -1\\ 1\\ 1 \end{bmatrix}, b=\begin{bmatrix} 2\\ 1\\ 1\\ 1 \end{bmatrix} \text{and}, c = \begin{bmatrix} 3\\ 1\\ 1\\ 1 \end{bmatrix}$$
   We define a partial basis as $\mathcal{X}$ which will be updated at each step.
   
   For the first step,
   we find $v_1$ as,
   $$\boxed{v_1 = \frac{a}{1} = \begin{bmatrix} 1\\ -1\\ 1\\ 1 \end{bmatrix}} \text{now the partial basis is,} \mathcal{X} = \big[\begin{bmatrix} 1\\ -1\\ 1\\ 1 \end{bmatrix} \big]$$
   similarly, we find $v_2$ as follows,
   $$v_2 = b - proj_{\mathcal{X}}\begin{bmatrix} 2\\ 1\\ 1\\ 1 \end{bmatrix}\\ = \begin{bmatrix} 2\\ 1\\ 1\\ 1 \end{bmatrix}-\frac{\begin{bmatrix} 2\\ 1\\ 1\\ 1 \end{bmatrix}.\begin{bmatrix} 1\\ -1\\ 1\\ 1 \end{bmatrix}}{\begin{bmatrix} 1\\ -1\\ 1\\ 1 \end{bmatrix}.\begin{bmatrix} 1\\ -1\\ 1\\ 1 \end{bmatrix}}\begin{bmatrix} 1\\ -1\\ 1\\ 1 \end{bmatrix} = \begin{bmatrix} 2\\ 1\\ 1\\ 1 \end{bmatrix}-\frac{3}{4}\begin{bmatrix} 1\\ -1\\ 1\\ 1 \end{bmatrix}\\ \text{Now,} \boxed{4v_2=\begin{bmatrix} 8\\ 4\\ 4\\ 4 \end{bmatrix}-\begin{bmatrix} 3\\ -3\\ 3\\ 3 \end{bmatrix}\\ =\begin{bmatrix} 5\\ 7\\ 1\\ 1 \end{bmatrix}}$$
   Now, we update $\mathcal{X}$,
   $$\mathcal{X} = \Big[\begin{bmatrix} 1\\ -1\\ 1\\ 1 \end{bmatrix}, \begin{bmatrix} 5\\ 7\\ 1\\ 1 \end{bmatrix} \Big]$$
   Now, we find $v_3$ as follows,
   $$v_3 = \begin{bmatrix} 3\\ 1\\ 1\\ 1 \end{bmatrix} - proj_{\mathcal{X}}\begin{bmatrix} 3\\ 1\\ 1\\ 1 \end{bmatrix}\\ v_3 = \begin{bmatrix} 3\\ 1\\ 1\\ 1 \end{bmatrix}-\frac{\begin{bmatrix} 3\\ 1\\ 1\\ 1 \end{bmatrix}.\begin{bmatrix} 1\\ -1\\ 1\\ 1 \end{bmatrix}}{\begin{bmatrix} 1\\ -1\\ 1\\ 1 \end{bmatrix}.\begin{bmatrix} 1\\ -1\\ 1\\ 1 \end{bmatrix}}\begin{bmatrix} 1\\ -1\\ 1\\ 1 \end{bmatrix}-\frac{\begin{bmatrix} 3\\ 1\\ 1\\ 1 \end{bmatrix}.\begin{bmatrix} 5\\ 7\\ 1\\ 1 \end{bmatrix}}{\begin{bmatrix} 5\\ 7\\ 1\\ 1 \end{bmatrix}.\begin{bmatrix} 5\\ 7\\ 1\\ 1 \end{bmatrix}}\begin{bmatrix} 5\\ 7\\ 1\\ 1 \end{bmatrix}\\ v_3 = \begin{bmatrix} 3\\ 1\\ 1\\ 1 \end{bmatrix}-\frac{4}{4}\begin{bmatrix} 1\\ -1\\ 1\\ 1 \end{bmatrix}-\frac{24}{76}\begin{bmatrix} 5\\ 7\\ 1\\ 1 \end{bmatrix}\\ \text{now},\\ \boxed{19v_3 = \begin{bmatrix} 8\\ -4\\ -6\\ -6 \end{bmatrix}}$$
   Hence, the partial basis is,
   $$\boxed{\mathcal{X} = \Big[\begin{bmatrix} 1\\ -1\\ 1\\ 1 \end{bmatrix}, \begin{bmatrix} 5\\ 7\\ 1\\ 1 \end{bmatrix}, \begin{bmatrix} 8\\ -4\\ -6\\ -6 \end{bmatrix}\Big]}$$

$$proj_{\overrightarrow{a}}{\overrightarrow{b}} = \frac{{\overrightarrow{a}}\cdot {\overrightarrow{b}}}{{\overrightarrow{a}}^2}{\overrightarrow{a}}$$

This partial basis is the orthogonal set, hence we now make it orthonoamal,
$$Q_{orthogonal} = \begin{bmatrix} 1 & 5 & 8\\ -1 & 7 & -2\\ 1 & 1 & -3\\ 1 & 1 & -3\\ \end{bmatrix}$$
Hence,
norm of first column, $$||n_1|| = \sqrt{1^2+(-1)^2+1^2+1^2} = 2$$
Similarly,
$$||n_2|| = \sqrt{76}$$
and
$$||n_2|| = \sqrt{38}$$
Hence, the orthonormal $Q$ is,
$$Q_{orthonormal}=\begin{bmatrix} 1/2 & 5/\sqrt{76} & 8/\sqrt{38} \\ -1/2 & 7/\sqrt{76} & -2/\sqrt{38}\\ 1/2 & 1/\sqrt{76} & -3/\sqrt{38}\\ 1/2 & 1/\sqrt{76} & -3/\sqrt{38}\\ \end{bmatrix}$$
Now, we have, $\mathcal{A}$ and $Q$; we now find $R$; which by theory has to be upper-triangular,
$$\mathcal{A}= Q \cdot R \\ Q^{T} \cdot \mathcal{A} = Q^{T} \cdot Q \cdot R\\ \text{NOTE:}\text{because } Q \text{ is orthonormal } \boxed{Q^{T} \cdot Q = I}.\\ \therefore \boxed{R = Q^{T} \cdot \mathcal{A}}$$
Make up the algorithm from the above steps.

---

4. Condition Number:[5]

As a rule of thumb, if the condition number $\kappa (A)=10^{k}$, then you may lose up to $k$ digits of accuracy on top of what would be lost to the numerical method due to loss of precision from arithmetic methods.

A problem with a low condition number is said to be well-conditioned, while a problem with a high condition number is said to be ill-conditioned. The condition number is a property of the problem.

Condition number is mathematically defined as,
$$\kappa(A) = ||A||_{1}||A^{-1}||_{1}$$
We define the $1$-norm as follows,
For a matrix $M$
$$M = \begin{bmatrix} 2 & 3\\ 1 & -1 \end{bmatrix}$$
The $1$-norm is,
$$\text{Maximum of }(2+1,3+1) = 4.\\ \text{Hence, }||M||_{1} = 4.$$

So, for a given matrix,
$$P = \begin{bmatrix} 1 & 10^{4}\\ -1 & 2 \end{bmatrix}$$
We get,
$$P^{-1} = \frac{1}{2+10^{4}}\begin{bmatrix} 2 & -10^{4}\\ 1 & 1 \end{bmatrix}$$
Hence,
$$||P||_{1} = \text{maximum of }(1+1, 2+10^{2}) = 2+10^{4}$$
and,
$$||P^{-1}||_{1} = \frac{1}{2+10^{4}}\text{maximum of }(2+1, 10^{4}+1) = \\ \frac{1}{2+10^{4}}(1+10^{4})$$
Therefore,
$$\text{Condition number},\\ \kappa(P) = ||P||_{1}||P^{-1}||_{1} = 1+10^{4}$$
Thus, the matrix (or system) is ill-conditioned.

---

The question paper that came for this exam was pretty sadistic.
Give it a try, if you can and tell me how awesome it feels. I will score around 30/50.

a FAILED mathematical approach into the search for the E.T.

Are we alone?

Scientists have been searching for another form of intelligent life for decades now. We have searched and have failed, we have the evidence and yet do not understand what they are.
Numerous Physicists, Astronomers, Conspiracy theorists and even common man ( :p ) have spent most or a significant time of their lives looking into the big blue sky and asking the question, “Am I alone?”
This sounds highly unlikely that we are the ONLY intelligent form of life in this huge universe.
I don’t care if the probability of any life in this universe, other than us, is $\rightarrow 0$, until and unless it is not $0$, I want to keep looking.
There have been so many things, let us take the Nazka lines of Peru. There are usually called Geoglyphs. These took a lot of time to come in light, Pedro Cieza de León mistook them as trail marks; some people mistook them as irrigation lines.
It was not until the days of flights and aeroplanes, when people saw these clearly; that is when they understood that these were not trail marks or anything. These were signs and drawings.
The one I like the most is as follows:



Now, what I ask is,
- How did they make these?
- Who did they make these for?
Unfortunately, these answers are not yet answered.

Some archeologists and geologists studied these lines and the results are shocking. These are not just lines with some nonsensical scribbling; these are made to last for a long time. Let alone the design; the methodology of design is scientific too. An excerpt from wikipedia is says the following,

On the ground, most of the lines are formed by a shallow trench with a
depth between 10 and 15 cm (4 and 6 in). Such trenches were made by
removing the reddish-brown iron oxide-coated pebbles that cover the
surface of the Nazca Desert. When this gravel is removed, the
light-colored clay earth which is exposed in the bottom of the trench
produces lines which contrast sharply in color and tone with the
surrounding land surface. This sublayer contains high amounts of lime
which, with the morning mist, hardens to form a protective layer that
shields the lines from winds, thereby preventing erosion.

Are these just drawings after all?
These are maps, some theorists say. Don’t believe them? have a look at the image below.




One of my favorite conspiracy theory is by Erich von Däniken, it is called the Cargo Cult

Erich von Däniken’s theory is the most famous approach to solve the
mystery of Nazca. He had the idea that long time ago visitors from
other stars visited the earth and naturally Nazca. At this place they
landed, during the landing stones was blown away by the power of
rocket propulsion. By approaching more the power was increasing and
the cleaned band broader. In this way the first trapezes emerged.
Later the Aliens disappeared and left confused people. Like in the
modern cargo cults they tried to call the Gods back by drawing lines,
figures and trapezes. Never Däniken said the formations was made by
Aliens. He discovered the GGF/Mandala/Zodiac and the mirror -
Formation and compares them with modern VASIS or PAPI-Signs.

To read more please visit : Conspiracy theories of Nazca Lines.

---

Mathematics:

There is a mathematical formula by American Astronomer and Astrophysicist Frank Drake [3]. He is hugely regarded as the father of modern SETI(Search for extraterrestrial intelligence).
The criticism related to Drake’s equation is not about the correctness of the equation itself, but about the highly ambiguous estimations of the various variable used in the equation.
The equation is as follows,
$$N = R^* \times f_p \times n_e \times f_l \times f_i \times f_c \times L$$
Here,
$N$ = the number of civilizations in our galaxy with which communication might be possible (i.e. which are on our current past light cone);
$R*$ = the average rate of star formation in our galaxy
$f_p$ = the fraction of those stars that have planets
$n_e$ = the average number of planets that can potentially support life per star that has planets
$f_l$ = the fraction of planets that could support life that actually develop life at some point
$f_i$ = the fraction of planets with life that actually go on to develop intelligent life (civilizations)
$f_c$ = the fraction of civilizations that develop a technology that releases detectable signs of their existence into space
$L$ = the length of time for which such civilizations release detectable signals into space

For people wondering what is a light cone, here is a picture. This is
a gloried way of looking into the dynamics of time.


In special and general relativity, a light cone is the path that a flash of light, emanating from a single event (localized to a single point in space and a single moment in time) and traveling in all directions, would take through spacetime

---

Criticism:

1. One major setback of this equation is that it takes no account of the cosmological developmental phases and time, the value of $N$ would heavily rely on the speed of cosmological development[1].
2. There is no derivation or methodology of how this equation came into being; which suggests that it is just a stab in the dark.
3. The factor $n_e$ determines habitability; how in the fuck’s world do you find that? Oxygen and Water makes Earth habitable for humans; and probably in some planet $\mathcal{X}$, Nitrogen and Alcohol makes it habitable [2].

---

Here is a beautiful image about Drake’s equation,

---

Conclusion:

1. This is a very badly written equation; not because it is wrong. But, because it has no methodology involved. Looks like couple of fresh college graduates wrote this while having some really potent weed.
   2.Drake’s equation is therefore a Fermi Paradox.

The Fermi paradox or Fermi’s paradox, named after physicist Enrico
Fermi, is the apparent contradiction between the lack of evidence and
high probability estimates, e.g., those given by the Drake equation,
for the existence of extraterrestrial civilizations.

---

An interesting coming from xkcd[3], it added another term $B_6$ which is the amount of bull shit you are gonna buy from the the drake equation.

---

Cheers!

Euler - Riemann Zeta Function.

“Madam, I have just come from a country where people are hanged if
they talk.” ― Leonhard Euler

The Riemann Zeta function, denoted as $\zeta{(s)}$ is a function of a complex variable $s$ that analytically continues the sum of the Dirichlet Series, [1]
$$\LARGE \zeta{(s)}= \sum_{n=1}^{\infty}\frac{1}{n^{s}}$$
The entire blog will be divided into the following parts:

1. What is Analytic Continuation?
2. What is Dirichlet Series?
3. Recurrence relation between Bernoulli Numbers.
4. Relationship between Bernoulli, Riemann and Euler.

---

1.

Before we talk about Analytic continuation, we will have to know what an Analytic function is.

There are multiple ways of defining the Analytic functions[2]. The one that I find most comfortable is given below:
$\underline{\text{Definition 1}}$: A function $f(z)$ is said to be analytic in a region $\mathcal{R}$ of the complex plane if $f(z)$ has a derivative at each point of $\mathcal{R}$ and if $f(z)$ is single valued.

For better understanding, I will extend an example,
$\underline{\text{Example 1}}$: Consider $f(z)$ = $2xy + i(x^{2}-y^{2})$, determine if the given function is analytic or not.
$\underline{\text{Solution}}$: We will use the Cauchy-Riemann equations.

For a given, $f(z) = f(x+iy) = u(x,y) + i v(x,y)$
Now we have the following definitions,
$$u_{x} = \frac{\partial{u}}{\partial{x}}, \\ u_{y} = \frac{\partial{u}}{\partial{y}}, \\ v_{x} = \frac{\partial{v}}{\partial{x}}, \\ v_{y} = \frac{\partial{v}}{\partial{y}}.$$
Now, according to the Cauchy-Riemann Equation, any $f(z)$ is an analytic function, iff,
$$u_{x}=v_{y}, \\ u_{y} = -v_{x}.$$

Now, for the given $f(z)$,
$$u(x,y) = 2xy, \\ v(x,y) = x^{2}-y^{2}.$$
Hence, the relations are,
$$u_{x} = \frac{\partial{(2xy)}}{\partial{x}} = 2y, \\ u_{y} = \frac{\partial{(2xy)}}{\partial{y}}=2x, \\ v_{x} = \frac{\partial{(x^{2}-y^{2})}}{\partial{x}}=2x, \\ v_{y} = \frac{\partial{(x^{2}-y^{2})}}{\partial{y}}=-2y.$$
It is clear that,
$$u_{x} \neq v_{y}, \\ u_{y} \neq -v_{x}.$$
Therefore, the given $f(z)$ is not an analytic function.

Analytic continuation is a pretty simple concept to understand really.
$\underline{\text{Definition 2}}$: Say $f(z)$ is an analytic function defined over a non-empty open subset $W$ of the complex plane $\mathcal{C}$. If $V$ is a larger open subset of $\mathcal{C}$, containing $W$, and $\mathcal{F}$ is an analytic function defined on $V$ such that,
$$\mathcal{F}(z) = f(z), \forall z \in V$$
In other words, Analytic continuation is the method of extending the domain of an analytic function. [3]

For some cool examples on Analytic continuation refer Virginia Tech’s paper here [4].

---

2.

Dirichlet Series[5] is any series of general form,
$$\sum_{n=1}^{\infty}\frac{a_n}{n^{s}}$$
Now, with a slight modification,
$$a_n = 1 , \forall n \in \mathcal{R}$$
We get,
$$\sum_{n=1}^{\infty}\frac{1}{n^{s}}$$
For a function defined as,
$$\zeta{(s)} = \sum_{n=1}^{\infty}\frac{1}{n^{s}}$$
is the Riemann Zeta function.

---

3.

Definitions,[6]
$$s_{0}(n)=n,\\ s_{1}(n) = \frac{n(n-1)}{2},\\ s_{2}(n) = \frac{n(n-1)(2n-1)}{6}...$$
Therefore, $s_{p}(n)$ is given as,
$$s_{p}(n) = \frac{1}{p+1}n^{p+1}-\frac{1}{2}n^{p}+\frac{p}{12}n^{(p-1)}+0 \times n^{(p-2)}+...$$
More generally,
$$s_{p}(n)=\sum_{k=0}^{p} \frac{B_{k}}{k!} \frac{p!}{(p+1-k)!}n^{(p+1-k)}$$
for $n=1$,
$$0=\sum_{k=0}^{p} \frac{B_{k}}{k!} \frac{p!}{(p+1-k)!}$$
or equivalently,
$$B_p=-\frac{1}{p+1}\sum_{k=0}^{p-1} \binom{p+1}{k}$$
The above relation is symbolically written as,
$$(B+1)^{(p+1)}-B^{(p+1)}=0$$
On expansion, all $k$-th powers of $B$, $B^{k}$ must be written as $B_k$ and treated as Bernoulli Numbers.
The expansion is done on the basis of Binomial Theorem, which statest that,
$$\forall a,b \in \mathcal{Z} \ \text{and} \ n \in \mathcal{N}\\ \large (a+b)^{n}= \sum_{k=0}^{n}\binom{n}{k}a^{n-k}b^{k}$$
Therefore, for $n=5$,
$$(a+b)^{5} = \sum_{k=0}^{5}\binom{5}{k}a^{5-k}b^{k}\\ \implies(a+b)^{5} = \binom{5}{0}a^{5-0}b^{0}+\binom{5}{1}a^{5-1}b^{1}+\binom{5}{2}a^{5-2}b^{2}+\\ \binom{5}{3}a^{5-3}b^{3}+\binom{5}{4}a^{5-4}b^{4}++\binom{5}{5}a^{5-5}b^{5}\\ \implies(a+b)^{5} = 1 \times a^{5}b^{0}+5\times a^{5-1}b^{1}+10\times a^{5-2}b^{2}+\\ 10\times a^{5-3}b^{3}+5\times a^{5-4}b^{4}+1\times a^{5-5}b^{5}\\ $$
Therefore, for $(B+1)^{(p+1)}-B^{(p+1)}=0$, after expanding for $p=4$ we have,
$$(B+1)^{(4+1)}-B^{(4+1)}=0 \\ \implies (B+1)^{5}-B^{5}=0 \\ \implies 1\times B_5+5\times B_4+10 \times B_3 + 10 \times B_2 + \\ 5 \times B_1+1\times B_0 =0\\ \implies 1\times 0+5\times B_4+10 \times 0+ 10 \times (\frac{1}{6}) + \\ 5 \times (\frac{-1}{2})+1\times 1 =0\\ \implies 5\times B_4+10 \times (\frac{1}{6}) + 5 \times (\frac{-1}{2})+1 =0\\ \implies B_4=-\frac{1}{30}\\ $$

---

4.

Euler found a formula which easily defined the **even-numbered zeta**functions as follows[7]:
$$\zeta(2n)=(-1)^{(n+1)}\frac{{B_{2n}}{(2\pi)^{2n}}}{2(2n)!}$$
Interestingly, $B_{2n+1}=0 \ \ \forall n \geq 1$, so there is no counter-part for $\zeta(2n+1)$ and yet a value of $\zeta(3)$ exists. Well, that is beyond my scope of this blog.
Cheers!