1) How does assuming any point as an initial approximation narrow down to the correct solutions?
2) why should the coefficient matrix be diagonally dominant?
3) if I were to think of such a method, what would my path of reasoning be?
Here's a similar question, but I don't quite get the answers. I specifically don't understand what the first two answers, assume to be evident,
which "naturally" defines a fixed-point iteration
You can do this in general for some system $F(X)=0$, $F=(f_1,...,f_n)$, $x=(x_1,...,x_n)$. The idea is to examine each equation and find its dominant variable, if there is one. If one can find a one-to-one dominance relationship of variables and equations, and if the dominance is strong enough, then you can solve each equation for its dominant variable. Assume that $x_i$ is the dominant variable in $f_i$.
Then one can loop through the equations in Jacobi-fashion, that is all-at-once in parallel (conceptually), $$ x_i^{new}=\text{ Solve }( 0=f_i(x^{old}_{1..i-1},u,x^{old}_{i+1..n}) \text{ for } u ) $$ or in Gauss-Seidel fashion, one-equation-at-a-time sequentially, replacing the computed variable value for the next equations $$ x_i^{new}=\text{ Solve }(0=f_i(x^{new}_{1..i-1},u,x^{old}_{i+1..n}) \text{ for } u) $$
This gives a fixed-point method $x^{new}=G(x^{old})$ and all the theorems about them apply, esp. the Banach fixed-point theorem. Thus it is a sensible demand that $G$ be contracting locally around the solution (which poses a problem in practical application as you do not know the solution).
In the non-linear case, there is no easy criterion of convergence, especially if you do not start close to a root (but one might transfer the convergence of the linearization). In the linear situation, contractivity of the method is a global property, so if it is present at all, then you get convergence from any initial point. As is known, one obtains this contractivity from diagonal dominance.