Finding the Shortest Way to Connect Three Points
Imagine you need to lay pipes, wires, or roads connecting three fixed locations to a single junction box, and you want to use as little material as possible. Where should that junction go? This is the classic Fermat–Torricelli point problem: given three points in a plane, find the point that minimizes the sum of distances to all three.
What makes this problem especially beautiful is that it has a direct mechanical analogy. If you drill three holes at the vertex positions on a horizontal board, thread a string through each hole, tie all three strings together in a single knot above the board, and hang equal weights from the other end of each string over the edge, the knot will settle exactly at the Fermat point. Physics finds the minimum for you.
The Math Behind It
Given three points $P_1, P_2, P_3$, we want to find point $P = (x, y)$ that minimizes:
$$
L(P) = \sum_{i=1}^{3} \lVert P - P_i \rVert
$$
Taking the gradient of $L$ with respect to $P$:
$$
\nabla L(P) = \sum_{i=1}^{3} \frac{P - P_i}{\lVert P - P_i \rVert}
$$
Each term is a unit vector pointing away from $P_i$ — physically, this is exactly the tension force of a string under equal load. At the minimum, these forces must cancel:
$$
\sum_{i=1}^{3} \hat{u}_i = 0
$$
Since all three vectors have the same magnitude (1) and must sum to zero, they can only be arranged one way: at 120° angles to each other. This is the famous geometric property of the Fermat point — as long as no angle of the triangle exceeds 120°, the point that satisfies this force-balance condition lies inside the triangle.
Setting $\nabla L(P) = 0$ and rearranging algebraically gives a fixed-point equation:
$$
P = \frac{\displaystyle\sum_{i=1}^{3} \frac{P_i}{\lVert P - P_i \rVert}}{\displaystyle\sum_{i=1}^{3} \frac{1}{\lVert P - P_i \rVert}}
$$
Iterating this equation is known as Weiszfeld’s algorithm, a classic and very fast way to solve this kind of “geometric median” problem.
Example Setup
Let’s use three concrete points, imagining them as three factories that need to be connected to a shared distribution hub:
- $A = (0, 0)$
- $B = (8, 0)$
- $C = (3, 7)$
We’ll solve for the Fermat point two different ways:
- Naive relaxation — a direct simulation of the physical system: the “knot” is dragged step by step in the direction of the net string tension (an overdamped, no-inertia model of the real strings-and-weights setup).
- Weiszfeld’s algorithm — the fast, mathematically derived fixed-point solver.
Comparing the two shows both why the mechanical picture works and how much faster a properly derived numerical method is than literally simulating physics step by step.
Full Source Code
1 | import numpy as np |
Code Walkthrough
Section 1–2 (setup and objective): We define the three points as a NumPy array and write total_length, which computes $L(P)$ using np.linalg.norm with axis=1 — this subtracts the candidate point from all three points at once and returns their norms in a single vectorized call, rather than looping in Python.
Section 3 (net_force): This function computes $\nabla L(P)$. diff holds the three vectors from each fixed point to the candidate point $P$; dividing each by its own length turns it into a unit vector — exactly the tension direction of a string under load. Summing them gives the net physical force acting on the imaginary knot. We clip the distance with np.maximum(dist, eps) purely to avoid a division-by-zero if $P$ ever lands exactly on one of the vertices.
Section 4 (naive relaxation): This is the literal simulation of the physical system: at each time step, the knot moves a small amount (lr) in the direction of the net force, just like an overdamped object being pulled by three strings with no inertia. It takes hundreds of small steps to settle down, which is why we call it “naive” — it mirrors real, gradual physical motion.
Section 5 (Weiszfeld’s algorithm): Instead of nudging the point step by step, this directly applies the fixed-point formula derived earlier from $\nabla L(P) = 0$. Each iteration recomputes $P$ as a distance-weighted average of the three fixed points — points that are closer pull harder, points that are farther pull more gently, exactly balancing at the solution. This converges in a tiny fraction of the iterations that naive relaxation needs.
Section 6–7 (running and verifying): We run both methods starting from the centroid of the triangle and confirm that they land on essentially the same point. We then measure the angles between the three lines connecting the Fermat point to $A$, $B$, and $C$ — mathematically, these should each equal 120°, confirming the force-balance condition we derived earlier.
Section 9 (2D plot): This draws the triangle in dashed gray, marks the three fixed points, draws solid crimson lines from the Fermat point to each vertex (the “strings”), and overlays both convergence paths — the naive relaxation (dotted blue, many points) and Weiszfeld’s path (orange, very few points) — so you can visually see how much faster the fixed-point method converges.
Section 10 (3D plot): Here we build a grid of $(x, y)$ values covering the area around the triangle and compute $L(x, y)$ at every grid point in a fully vectorized way (looping only 3 times — once per fixed point — never once per grid cell). The resulting surface is a bowl-shaped function with a single global minimum, which we mark in red at the Fermat point. This is a great visual way to see that the optimization problem is well-behaved: one smooth valley, one clear minimum, no risk of getting stuck in a wrong answer.
About Performance
Both methods here are extremely lightweight — even a few hundred iterations on two-dimensional vectors run in well under a millisecond, so no further speed optimization is strictly necessary for this three-point example. That said, the code already demonstrates the general principle for scaling up: total_length, net_force, and the surface computation in Section 10 are all vectorized with NumPy rather than using per-point Python loops, and Weiszfeld’s algorithm converges in a small constant number of steps regardless of the geometry, so the same code would scale comfortably even if you needed to repeat this calculation for thousands of different triangles (for example, optimizing hub placement across many separate clusters of three locations at once).
Visualizing the Result
Run the code above in a fresh cell. It will print a summary to the console and generate two figures.
2D result — the triangle, the strings, and both convergence paths:
3D result — the total-length surface with the minimum marked:
Console output — point coordinates, iteration counts, and the 120° angle check:
======================================================= Fermat-Torricelli point ======================================================= A=[0. 0.], B=[8. 0.], C=[3. 7.] Naive relaxation -> point=[3.34728069 2.2650171 ], iterations=400, total length=13.964065 Weiszfeld (fast) -> point=[3.34017112 2.26202749], iterations=43, total length=13.964055 Angle A-F-B = 120.000 deg Angle B-F-C = 120.000 deg Angle C-F-A = 120.000 deg Sum of angles = 360.000 deg (should be ~360)
Once you drop in your own screenshots and console text above, look closely at the 2D plot: you should see the orange Weiszfeld path collapse onto the Fermat point in just a handful of steps, while the blue dotted naive-relaxation path spirals in much more gradually — a nice visual confirmation that solving the force-balance equation directly is far more efficient than simulating the physical settling process step by step, even though both arrive at exactly the same 120°-balanced point.
