From $A(x_1,y_1,z_1),B(x_2,y_2,z_2),C(x_3,y_3,z_3)$ we can get their position vectors.
$\vec{AB}=(x_2-x_1,y_2-y_1,z_2-z_1)$ and $\vec{AC}=(x_3-x_1,y_3-y_1,z_3-z_1)$.
Then $||\vec{AB}\times\vec{AC}||=0\implies A,B,C$ collinear.
The regular Newton-Raphson method is initialized with a starting point $x_0$ and then you iterate $x_{n+1}=x_n-\dfrac{f(x_n)}{f'(x_n)}$.
In higher dimensions, there is an exact analog. We define:
$$F\begin{bmatrix}x\\y\\z\end{bmatrix} = \begin{bmatrix}f_1(x,y,z) \\ f_2(x,y,z) \\ f_3(x,y,z)\end{bmatrix} = \begin{bmatrix}\sqrt{(x-20000)^2+(y-19400)^2+(z-19740)^2}-20000 = 0 \\ \sqrt{(x-18700)^2+(y-1800)^2+(z-18500)^2}-19000 = 0 \\ \sqrt{(x-18900)^2+(y-17980)^2+(z-20000)^2}-19500= 0 \end{bmatrix}$$
The derivative of this system is the $3x3$ Jacobian given by:
$J(x,y,z) = \begin{bmatrix} \dfrac{\partial f_1}{\partial x} & \dfrac{\partial f_1}{\partial y} & \dfrac{\partial f_1}{\partial z}\\ \dfrac{\partial f_2}{\partial x} & \dfrac{\partial f_2}{\partial y} & \dfrac{\partial f_2}{\partial z} \\ \dfrac{\partial f_3}{\partial x} & \dfrac{\partial f_3}{\partial y} & \dfrac{\partial f_3}{\partial z}\end{bmatrix} =\\
\begin{bmatrix}
\frac{x-20000}{\sqrt{(x-20000)^2+(y-19400)^2+(z-19740)^2}} & \frac{y-19400}{\sqrt{(x-20000)^2+(y-19400)^2+(z-19740)^2}} & \frac{z-19740}{\sqrt{(x-20000)^2+(y-19400)^2+(z-19740)^2}} \\
\frac{x-18700}{\sqrt{(x-18700)^2+(y-1800)^2+(z-18500)^2}} & \frac{y-1800}{\sqrt{(x-18700)^2+(y-1800)^2+(z-18500)^2}} & \frac{z-18500}{\sqrt{(x-18700)^2+(y-1800)^2+(z-18500)^2}} \\
\frac{x-18900}{\sqrt{(x-18900)^2+(y-17980)^2+(z-20000)^2}} & \frac{y-17980}{\sqrt{(x-18900)^2+(y-17980)^2+(z-20000)^2}} & \frac{z-20000}{\sqrt{(x-18900)^2+(y-17980)^2+(z-20000)^2}} \\
\end{bmatrix}$
The function $G$ is defined as:
$$G(x) = x - J(x)^{-1}F(x)$$
and the functional Newton-Raphson method for nonlinear systems is given by the iteration procedure that evolves from selecting an initial $x^{(0)}$ and generating for $k \ge 1$,
$$x^{(k)} = G(x^{(k-1)}) = x^{(k-1)} - J(x^{(k-1)})^{-1}F(x^{(k-1)}).$$
We can write this as:
$$\begin{bmatrix}x^{(k)}\\y^{(k)}\\z^{(k)}\end{bmatrix} = \begin{bmatrix}x^{(k-1)}\\y^{(k-1)}\\z^{(k-1)}\end{bmatrix} + \begin{bmatrix}y_1^{(k-1)}\\y_2^{(k-1)}\\y_3^{(k-1)}\end{bmatrix}$$
where:
$$\begin{bmatrix}y_1^{(k-1)}\\y_2^{(k-1)}\\y_3^{(k-1)}\end{bmatrix}= -\left(J\left(x^{(k-1)},y^{(k-1)},z^{(k-1)}\right)\right)^{-1}F\left(x^{(k-1)},y^{(k-1)},z^{(k-1)}\right)$$
Here is a starting vector:
$$x^{(0)} = \begin{bmatrix}x^{(0)}\\y^{(0)}\\z^{(0)}\end{bmatrix} = \begin{bmatrix}1\\1\\1\end{bmatrix}$$
The iterates will be:
- $x_0 = (1,1,1)$
- $x_1 = (17289.9,15734.3,-8508.52)$
- $x_2 = (13614.9,9566.68,1369.85)$
- $x_3 = (16370.9,11019.8,1760.91)$
- $x_4 = (16274.,10923.6,2010.87)$
- $x_5 = (16276.2,10924.5,2011.53)$
- $x_6 = (16276.2,10924.5,2011.53)$
You should end up with a solution that looks something like:
$$\begin{bmatrix}x\\y\\z\end{bmatrix} = \begin{bmatrix} 16276.2 \\ 10924.5 \\ 2011.53 \end{bmatrix}$$
Of course, for a different starting vector you could get a different solution and perhaps no solution at all.
Best Answer
Your reference triangle lives naturally in 2D. So I'd define it as $$T_{ref} = \{(0,0),(1,0),(0,1)\}.$$ Now your transformation is $$ B\begin{bmatrix}\hat{x}\\\hat{y}\end{bmatrix} + c = \begin{bmatrix} x_1-x_3 & x_2-x_3 \\ y_1-y_3 & y_2-y_3 \\ z_1-z_3 & z_2-z_3 \end{bmatrix} \begin{bmatrix}\hat{x}\\\hat{y}\end{bmatrix} + \begin{bmatrix}x_3\\y_3\\z_3\end{bmatrix}. $$ To translate an integral over the reference triangle to an integral over the 3D triangle, you've got to account for the degree to which $B$ distorts area. To do so, note that $B$ maps the square in 2D to a parallelogram spanned by the columns of $B$ in 3D. Thus, the area distortion is simply the area of this parallelogram which can be computed as the magnitude of the cross-product of those column vectors.
Note that, more generally, if $\vec{p}:U\to\mathbb R^3$ maps a region $U$ in a $uv$-plane to $\mathbb R^3$, then the local area distortion of $\vec{p}$ is $$||\vec{p}_u \times \vec{p}_v||.$$ This is the standard distortion factor that's used when computing surface integrals, as described on Wikipedia.