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
Assume $r_1 = 0, r_2, r_3 \ne 0$. Let $t_1 \tan\frac{\psi}{2}$, $t_2 = \tan\frac{\phi}{2}$ and $t_3 = \tan\frac{\theta}{2}$.
With help of a CAS, I get two set of solutions: $$t_2 = \frac{s_2 - \epsilon \sqrt{s_2^2 - 4r_2^2}}{2r_2},\quad t_3 = \frac{s_3 + \epsilon \sqrt{s_3^2 + 4r_3^2}}{2r_3}\tag{*1}$$ where $\epsilon = \pm 1$ and $s_2 = r_2^2+r_3^2+1$, $s_3 = r_2^2+r_3^2-1$.
For $(r_1,r_2,r_3) = (0,−0.101233861621737,0.365119069777688)$, the approximation solution you get corresponds to the $\epsilon = +1$ solution. Numerically, the corresponding angles evaluated to
$$(\psi, \phi, \theta ) \sim (-0.0657294323,-0.1779886273,0.7060269791)$$
The formula in $(*1)$ is relatively simple, it should be possible to derive that by hand.
However, I can't figure that out yet. You can plug those expression to a CAS and verify that yourself.Update
I sort of figure out how to derive the two solutions in $(*1)$.
When $r_1 = 0$, first equation give us $t_1 = t_2t_3$. Substitute this into second and third equation, we have
$$\frac{t_2(1+t_3^2)}{1 + t_2^2t_3^2} = r_2\quad\text{ and }\quad \frac{t_3(1-t_2^2)}{1 + t_2^2t_3^2} = r_3\tag{*2}$$
By brute force, we have $$s_2 = r_2^2 + r_3^2 + 1 = \frac{(1+t_2^2)(1+t_3^2)}{1+t_2^2t_3^2} \implies r_2(1+t_2^2) - s_2t_2 = 0 $$ Similarly, $$s_3 = r_2^2 + r_3^2 - 1 = -\frac{(1-t_2^2)(1-t_3^2)}{ 1 + t_2^2t_3^2} \implies r_3(t_3^2-1) - s_3t_3 = 0$$
Solving these two quadratic equations give us $4$ possible combinations of $(t_2,t_3)$. Throwing them into an CAS, one find only two of them, those appear in $(*1)$, satisfy $(*2)$.