Why is it that by working backwards from the euclidean algorithm one can find the modular inverse of a number?
Further, there is also another method for finding inverses discussed here which seems similar to the extended euclidean algorithm method but is much shorter. How is this method related to the extended euclidean algorithm one, and why does this method work?
Using the Euclidean algorithm for two coprime numbers gives you $$a=q_1b+r_1\ ,\quad b=q_2r_1+r_2\ ,\ldots,\quad r_{n-2}=q_nr_{n-1}+r_n\ ,$$ where $r_n=1$. Therefore $$1=r_{n-2}-q_nr_{n-1}\ ,$$ and then by substituting for $r_{n-1}$ in terms of $r_{n-2}$ and $r_{n-3}$, and so on, you eventually get $$1=ax+by$$ for some $x,y$. Doing a few examples should convince you that this always works; if you want a more formal proof, you can do it by induction on $n$, the number of steps in the algorithm. We now have $$1\equiv ax\pmod b\ ,$$ which means by definition that $x$ is the inverse of $a$ modulo $b$.
With regard to the other algorithm that you linked, if you do it the Euclidean algorithm way starting with the same numbers $811$ and $216$ you will see that the remainders you get are just the same as in the other method: this is why they are related, and it works for essentially the same reason as above. Once again, you could prove it formally by induction.
By the way, I'm not convinced that the other method is very much shorter than the Euclidean method - maybe a little. Note that in the page you linked they just wrote down the values of $f,e,d,c,b,a$, but there is actually some calculation to be done here which they didn't show. If you do the same example using the "reverse Euclidean" method I think you will find that the arithmetic work is fairly similar.
The induction step simply lifts up the Bezout identity for the gcd by $\rm\color{#c00}{eliminating}$ the remainder
$\quad \color{#0a0}{a = q a'}\color{#c00}{ + a''} \Rightarrow\, \gcd(a,a') = \gcd(a',a'')$
$\quad\ \ {\rm induction}\ \Rightarrow\,\ \gcd(a',a'') = j\, a' + k\, \color{#c00}{a''} = {\rm linear\ combination\ of}\ a',a''\ {\rm (Bezout)}$
$\quad\begin{eqnarray}{\rm elmination}\,\Rightarrow\ \gcd(a,a') &\,=\,& j\, a' + k (\color{#0a0}{a\!-\!qa'})\ \ {\rm by}\ \ \rm\color{#c00}{eliminating\,\ a''} =\, \rm remainder\\ \\ &=& k\, a + (j\! -\!kq)\, a' = {\rm linear\ combination\ of}\ a,a'\ {\rm (Bezout)}\end{eqnarray}$
Therefore, when the gcd $\, = 1 = k\, a + n\, a'\,$ then $\, 1 \equiv k a\pmod{a'},\,$ so $\ a^{-1}\equiv k.$
The method in the linked article employs the same Euclidean remainder sequence as the standard back-substitution extended Euclidean algorithm. The only difference is that it explicitly writes out the equations connecting the steps, which apparently makes the back substituion simpler for you.
The back-substution process is notoriously error-prone. A better way to eliminate the complex back-substitution process is to forward propagate the Bezout identities by using $ $ row operations, see this answer. For further optimization one can employ least magnitude remainders, i.e. allow negative remainders, which generally halves the length of the remainder sequence, e.g. here.
