For my computer science project I am working on a calculator that will find(or approximate) the point of equidistance along the surface of the Earth from N number of points. Using basic vector math I can quite easily find this equidistant point (which I will be calling the centre from now on) for N=2 and N=3 points, however for higher values of N I will need to approximate.
Imagining the Earth as a sphere with a constant radius in a 3 dimensional plane, we know that the arclength along the surface of a sphere can be calculated using:
$$l = r\arctan\frac{|n_{1}\times n_{2}|}{n_{1} \cdot n_{2}}$$
where
- $l$ is arclength across the surface of the Earth
- $n_{1}, n_{2}$ are two points on the surface of the Earth represented as 3 dimensional vectors
- $r \approx 6371$ as radius of the Earth
This is gathered from https://en.wikipedia.org/wiki/Great-circle_distance#Vector_version. Points are taken in cartesian form (calculated from polar form).
From here, I imagined the attributes of what each point on the surface of the Earth (vertex) and their distance from the unknown centre may have in relation to each other if they were to be equidistant. This can best be shown by the variance distances between each vertex and the centre. Variance is defined as:
$$\sigma^2 = \frac{\sum_{i=1}^{n}(x_i - \mu)^2} {n}$$
But, an easier form to work with would be:
$$\sigma^2 = \frac{\sum_{i=1}^{n}x_{i}^2} {n} - \Bigl(\frac{\sum_{i=1}^{n}x_i}{n}\Bigr)^2$$
Then we can see that arclength between a general vertex (V) and the centre (C) is written as:
$$\arctan\frac{|V\times C|}{V \cdot C}$$
The radius above is ignored as it acts as a constant that doesn't affect the maths going on behind the method and only the end result (now the problem is effectively minimizing the azimuth angle between vectors but that's a less catchy title). Now I substituted that into the variance formula to give this beautiful jumble of letters:
$$\sigma^2 = \frac{\sum_{i=1}^{n}arctan^2\frac{|V_i\times C|}{V_i \cdot C}} {n} - \Bigl(\frac{\sum_{i=1}^{n}arctan\frac{|V_i\times C|}{V_i \cdot C}}{n}\Bigr)^2$$
With this expression then, you can take the derivative of variance with respect to the centre, and by setting the derivative equal to zero, we can minimize the equation. This means now doing vector calculus, and this is where I have likely gone horribly wrong (Note my multivariable calculus is self-taught as I am only a year 13 student, so apologies if this is poorly done).
First setting each vertex to be defined as a 3d vector:
$$V_i = [V_{i1}, V_{i2}, V_{i3}]^T$$
And similarly to the centre:
$$C = [x, y, z]^T$$
Then I set about finding the derivative of just arclength:
$$l = \arctan(u) \\ u = \frac{v}{w} \\ v = |V_i\times C| \\ w = V_i \cdot C$$
From here, I worked out derivatives back in terms of $V_i$ and $C$ which gave me these:
$$\frac{\partial w}{\partial C} = V_i$$ $$\frac{\partial v}{\partial C} = \frac{1}{|V_i\times C|}K_iC, \\ K_i = \begin{pmatrix} V_{i2}^2 + V_{i3}^2 & -V_{i1}V_{i2} & -V_{i1}V_{i3}\\ -V_{i1}V_{i2} & V_{i1}^2 + V_{i3}^2 & -V_{i2}V_{i3}\\ -V_{i1}V_{i3} & -V_{i2}V_{i3} & V_{i1}^2 + V_{i2}^2 \end{pmatrix}$$ $$\frac{\partial u}{\partial C} = \frac{1}{(V_i \cdot C)^2}\Bigl[\frac{(V_i \cdot C)}{|V_i\times C|}K_iC - V_i|V_i\times C|\Bigr]$$ $$\frac{\partial l}{\partial C} = \frac{1}{|V_i|^2|C|^2}\Bigl[\frac{(V_i \cdot C)}{|V_i\times C|}K_iC - V_i|V_i\times C|\Bigr]$$
This was done with an application of the chain and product rule, as well as using the basics of multivariable calculus, and the use of the identity: $|V_i|^2|C|^2 = (V_i \cdot C)^2 + |V_i\times C|^2$ From here, I replaced every arctan in the variance formulae with $l$ and took the derivative of variance with respect to $l$, and multiplied expectedly.
$$\frac{\partial \sigma^2}{\partial C} = \frac{2}{n}\sum_{i=1}^{n}\frac{1}{|V_i|^2|C|^2}\Bigl[\frac{(V_i \cdot C)}{|V_i\times C|}K_iC - V_i|V_i\times C|\Bigr]\arctan\frac{|V_i\times C|}{V_i \cdot C} \\ -\frac{2}{n^2} \Bigl(\sum_{i=1}^{n}\arctan\frac{|V_i\times C|}{V_i \cdot C}\Bigr)\Bigl(\sum_{i=1}^{n}\frac{1}{|V_i|^2|C|^2}\Bigl[\frac{(V_i \cdot C)}{|V_i\times C|}K_iC - V_i|V_i\times C|\Bigr]\Bigr)$$
Which as you can see is an ugly mess. I thought at this point I was finished, as I could rearrange this to separate out a C by using inverse matrices, setting the derivative equal to zero, and define an iterative formula from that; however when programming that in, I end up resulting in a 0 = 0 scenario.
I thought that was it for the project, but then only recently in my lessons I learnt about the Newton Raphson method of approximation (yes only recently still a secondary school student), and that may serve as the solution to my problem. The Newton Raphson method can be applied in n dimensions with the formula:
$$x_{n+1} = x_n - J^{-1}(x_n)f(x_n)$$
Where $J^{-1}(x_n)$ is the inverse Jacobian of $f(x_n)$. Calculating the inverse of a matrix is quite easy through numpy in the code, which means all I need is a way to write the inverse Jacobian of $f(x_n)$. This effectively means finding the second derivative of my variance formula, and this is where I have hit a block. I understand the layout of a Jacobian matrix, and how to calculate them if they were just given in terms of the direction vectors of C, but when given in such a vector format, I don't know how to differentiate the vector with respect to a vector.
So, I would really appreciate any help anyone can offer on the solution; and in truth the answer would just be as nice (since whoever marks my computing project will not understand this part of the maths anyway).
Oh and just to add, I have tried doing it through Wolfram Cloud/Mathematica and the format it gives it in is unusable as I really need it in terms of the vectors.
Thank you for any responses in advance, sorry for any LaTeX or maths errors.