Disclaimer: I am not actually a math guy. The main reason I wanted to understand barycentric coordinates is because I needed to use them in my software renderer. I've implemented them in code, but now I really want to understand what happens „under the hood“ of this math concept. So my math might be a bit naive but I try to do my best.
The definition of the barycentric coordinates is the following:
Definition. Let $\mathbf{A}$ be an Euclidean space, a flat, or an affine space and the points $V_0, \ldots, V_n \in \mathbf{A}$. Then for any given point $P \in \mathbf{A}$ there are scalars $k_0, \ldots, k_n$ such that $\exists \; k_i \neq 0$ and
$$ \sum k_i \; \overrightarrow{OP} = \sum k_i \; \overrightarrow{OV_i} $$
for any $O \in \mathbf{A}$. These scalaras $k_0, \ldots, k_n$ are called the barycentric coordinates of the point $P$ with respect to the point $V_0, \ldots, V_n$.
But I also have seen that the barycentric coordinates of a point $P$ inside a triangle $ABC$ are expressed as a tuple $(\alpha, \beta, \gamma)$ such as
$$ \alpha = \frac{Area \; \triangle BCP}{Area \; \triangle ABC}, \;\;\; \beta = \frac{Area \; \triangle ACP}{Area \; \triangle ABC}, \;\;\; \gamma = \frac{Area \; \triangle ABP}{Area \; \triangle ABC}. $$
and
$$ P = \alpha A + \beta B + \gamma C. $$
So I want to prove that
Theorem. Let $\mathbf{A}$ be an Euclidean space, a flat, or an affine space; points $A, B, C \in \mathbf{A}$ form a triangle, and point $P \in \mathbf{A}$. Then scalars $\alpha, \beta, \gamma$ such that
$$ \alpha = \frac{Area \; \triangle BCP}{Area \; \triangle ABC}, \;\;\; \beta = \frac{Area \; \triangle ACP}{Area \; \triangle ABC}, \;\;\; \gamma = \frac{Area \; \triangle ABP}{Area \; \triangle ABC} $$
are the barycentric coordinates of the point $P$ with respect to the $\triangle ABC$.
As I see it in order to prove this theorem I need to find a general formula of the barycentric coordinates of a point and then express the areas of the triangles through their vertices ($A, B, C$) and if these results will be equal than the theorem will be proved.
I also think that the solution will not depend on the dimension of $\mathbf{A}$ but it would be easier to start with a 2D flat and then generalize the prove to a $n$-dimensional space.
Proof. First we need to express the triangles areas using their vertices coordinates. It is easy to do using the shoelace formula:
$$ Area \; \triangle ABC = \frac{1}{2} \; \left | \; det \begin{pmatrix} A_x & B_x & C_x \\ A_y & B_y & C_y \\ 1 & 1 & 1 \end{pmatrix} \right | = \frac{1}{2} (A_x B_y + A_y C_x + B_x C_y - B_y C_x - C_y A_x - A_y B_x), \\ $$
$$ Area \; \triangle BCP = \frac{1}{2} \; \left | \; det \begin{pmatrix} B_x & C_x & P_x \\ B_y & C_y & P_y \\ 1 & 1 & 1 \end{pmatrix} \right | = \frac{1}{2} (B_x C_y + B_y P_x + C_x P_y - C_y P_x - P_y B_x - B_y C_x), \\ $$
$$ Area \; \triangle ACP = \frac{1}{2} \; \left | \; det \begin{pmatrix} A_x & C_x & P_x \\ A_y & C_y & P_y \\ 1 & 1 & 1 \end{pmatrix} \right | = \frac{1}{2} (A_x C_y + A_y P_x + C_x P_y - C_y P_x - P_y A_x - A_y C_x), \\ $$
$$ Area \; \triangle ABP = \frac{1}{2} \; \left | \; det \begin{pmatrix} A_x & B_x & P_x \\ A_y & B_y & P_y \\ 1 & 1 & 1 \end{pmatrix} \right | = \frac{1}{2} (A_x B_y + A_y P_x + B_x P_y - B_y P_x - P_y A_x - A_y B_x). $$
Then we can express the scalars $\alpha, \beta, \gamma$ as follows:
$$ \alpha = \frac{Area \; \triangle BCP}{Area \; \triangle ABC} = \frac{B_x C_y + B_y P_x + C_x P_y - C_y P_x - P_y B_x - B_y C_x}{A_x B_y + A_y C_x + B_x C_y - B_y C_x - C_y A_x - A_y B_x}, $$ $$ \beta = \frac{Area \; \triangle ACP}{Area \; \triangle ABC} = \frac{A_x C_y + A_y P_x + C_x P_y - C_y P_x - P_y A_x - A_y C_x}{A_x B_y + A_y C_x + B_x C_y - B_y C_x - C_y A_x - A_y B_x}, $$ $$ \gamma = \frac{Area \; \triangle ABP}{Area \; \triangle ABC} = \frac{A_x B_y + A_y P_x + B_x P_y - B_y P_x - P_y A_x - A_y B_x}{A_x B_y + A_y C_x + B_x C_y - B_y C_x - C_y A_x - A_y B_x}. $$
Then we need to find the „real“ barycentric coordinates of the point $P$ respect to the points $A, B, C$. Let's call them $u, v, w$:
$$ P = uA + vB + wC. $$
Since we are in 2D now we can rewrite this equation as a system of two linear equations:
$$ \begin{cases} A_x u + B_x v + C_x w = P_x \\ A_y u + B_y v + C_y w = P_y \end{cases} $$
And here I am stuck. I can't to solve this system using Gauss elimination neither Cramer's method or inverse matrix (because in this system the number of the variables is not equal to the number of the equations).
Thanks @daw. I managed to solve the problem using his hint.
Since $v + u + w = 1$:
$$ \begin{cases} A_x (1 - v - w) + B_x v + C_x w = P_x \\ A_y (1 - v - w) + B_y v + C_y w = P_y \end{cases} $$
What leads us to a system with two equations and two variables:
$$ \begin{cases} \overrightarrow{AB}_x v + \overrightarrow{AC}_x w = \overrightarrow{AP}_x \\ \overrightarrow{AB}_y v + \overrightarrow{AC}_y w = \overrightarrow{AP}_y \end{cases} $$
Which can be solved using Cramer's rule:
$$ v = \frac{\det \begin{pmatrix} \overrightarrow{AP}_x & \overrightarrow{AC}_x \\ \overrightarrow{AP}_y & \overrightarrow{AC}_y \end{pmatrix}} {\det \begin{pmatrix} \overrightarrow{AB}_x & \overrightarrow{AC}_x \\ \overrightarrow{AB}_y & \overrightarrow{AC}_y \end{pmatrix}} = \frac{\overrightarrow{AP}_x \overrightarrow{AC}_y - \overrightarrow{AP}_y \overrightarrow{AC}_x}{ \overrightarrow{AB}_x \overrightarrow{AC}_y - \overrightarrow{AB}_y \overrightarrow{AC}_x }, \\ w = \frac{\det \begin{pmatrix} \overrightarrow{AB}_x & \overrightarrow{AP}_x \\ \overrightarrow{AB}_y & \overrightarrow{AP}_y \end{pmatrix}} {\det \begin{pmatrix} \overrightarrow{AB}_x & \overrightarrow{AC}_x \\ \overrightarrow{AB}_y & \overrightarrow{AC}_y \end{pmatrix}} = \frac{\overrightarrow{AB}_x \overrightarrow{AP}_y - \overrightarrow{AB}_y \overrightarrow{AP}_x}{ \overrightarrow{AB}_x \overrightarrow{AC}_y - \overrightarrow{AB}_y \overrightarrow{AC}_x }. $$
Then we can go from vectors to scalars:
$$ \begin{align*} v &= \frac{(P_x - A_x)(C_y - A_y) - (P_y - A_y)(C_x - A_x)} {(B_x - A_x)(C_y - A_y) - (B_y - A_y)(C_x - A_x)} \\ &= \frac{ C_y P_x + P_y A_x + A_y C_x - A_x C_y - A_y P_x - C_x P_y }{ A_x B_y + A_y C_x + B_x C_y - B_y C_x - C_y A_x - A_y B_x }, \\ w &= \frac{(B_x - A_x)(P_y - A_y) - (B_y - A_y)(P_x - A_x)} {(B_x - A_x)(C_y - A_y) - (B_y - A_y)(C_x - A_x)} \\ &= \frac{ B_x P_y + B_y A_x + A_y P_x - B_x A_y - A_x P_y - B_y P_x }{ A_x B_y + A_y C_x + B_x C_y - B_y C_x - C_y A_x - A_y B_x }. \end{align*} $$
And then we can find $u$:
$$ u = 1 - v - w = \frac{B_x C_y + B_y P_x + C_x P_y - C_y P_x - P_y B_x - B_y C_x}{A_x B_y + A_y C_x + B_x C_y - B_y C_x - C_y A_x - A_y B_x}. $$
As we can see
$$ u = \alpha,\;\; v = \beta,\;\; w = \gamma. $$