Problem Description#
Two random points, one red and one blue, are chosen uniformly and independently from the interior of a square. To ten decimal places, what is the probability that there exists a point on the side of the square closest to the blue point that is equidistant to both the blue point and the red point?
You can check out the problem on their website, however if you want to check their solution be aware that as of writing there is a mistake in their proposed integral, so please refer to the one I have provided. I have made them aware of this issue.
Solution#
Mathematical representation of the problem#
For this problem, we need to consider which parts of the selected wall are closest to blue, closest to red, and an equal distance.
Let’s plot some examples, highlighting the part of the relevant side which is closest to red or blue:
As seen above, for any two points, there are two possibilites:
- one point is closest to the whole side
- each point has a share of the side
If one point is closest to the whole side, then there cannot be a point on that side which is equidistant from both. However, in the 2nd case, there must be a point on the side equidistant of the two, located where the two highlighted areas meet.
To check these cases, we check the two corners. If the closest point to each corner is the same, then that point is closest to the whole side, but if one corner is closest to blue and the other closest to red, then each point has a share of the side and there will be a solution.
Let’s fix the blue point, and consider its distance from each corner:
The red point needs to closer to one corner, and not the other. We can find the boundary for being closer to a corner by drawing a circle from it:
To get a solution the red point must be inside one circle, but not both:
Recap#
We have found an alternative way to represent the original problem as the probability of the red point being in the highlighted red area. From now on we will define the total area of the square to be 1, so if the red area is 0.63 then the probability of there being a valid solution is 63%.
The question asks about the average probability across all possible blue point positions, so the answer to the question is the average red area.
Creating an equation to solve#
simplifying the representation#
Due to the symmetry of the square, we can choose to view the blue point from the perspective where its closest side is along the x-axis, and its closest corner is along the y axis.
Probability for a single blue position#
To find the area, we can add the areas of the individual quarter-circles, and then subtract their intersection twice (one time to account for double counting, and the second to actually remove the area from the valid space).
Let \(r_c,r_f\) be the distance from the blue point to the closest and furthest corner of its side respectively.
Let \(p(x,y)\) be the probability of a valid solution given the blue point is at position \((x,y)\).
$$ p(x,y) = \frac{\pi}{4}(r_c^2+r_f^2) - 2*\text{intersection} $$
To find the area of intersection, we add together the segments from each circle that makes it, and subtract the total area formed by the triangle between the two corners and \((x,y)\). $$ \text{intersection} = \frac{\alpha}{2 \pi} \pi r_c^2 + \frac{\beta}{2 \pi} \pi r_f^2 - \frac{y}{2} $$
where $$ \alpha = \arctan(\frac{y}{x}) $$ $$ \beta = \arctan(\frac{y}{1-x}) $$ $$ r_c = \sqrt{x^2+y^2} $$ $$ r_f = \sqrt{(1-x)^2+y^2} $$
mean probability for all blue positions#
To get the average of our probability function for all blue points, we integrate over all blue points and divide by the area of the bounds.
The bounds on the blue point’s position are represented by the cyan triangle seen earlier. This triangle has a slope that follows the gradient \(y=x\), so we represent the triangle using the following inequalities: $$ 0 \le y \le x $$ $$ 0 \le x \le 0.5 $$
The area of these bounds is \(\frac{1}{8}\), so the mean probability (and the answer to the question) is given by: $$ 8\int_{0}^{\frac{1}{2}} \int_{0}^{x} p(x,y) \ dx \ dy $$
Solving the equation#
Given that we don’t need to find the exact solution to the problem, we can use numerical methods to approximate the integral:
from scipy.integrate import dblquad
from math import atan, pi
def p(y,x):
if x==0 or y==0:
return 0
r1_sq = x**2 + y**2
r2_sq = (1-x)**2 + y**2
area = pi / 4 * (r1_sq+r2_sq)
intersection = 1/2*atan(y/x)*r1_sq + 1/2*atan(y/(1-x))*r2_sq - y/2
return area - 2*intersection
integral_val= dblquad(p, 0.0, 0.5, lambda x: 0.0, lambda x: x)[0]
result = integral_val*8
This outputs 0.49140757883830793
However, we can do better! Using python’s sympy library, we can calculate the exact result of the integral:
import sympy as sp
x,y = sp.symbols('x y')
r1_sq = x**2 + y**2
r2_sq = (1-x)**2 + y**2
area = sp.pi / 4 * (r1_sq+r2_sq)
intersection = sp.Rational(1,2)*(sp.atan(y/x)*r1_sq + sp.atan(y/(1-x))*r2_sq - y)
integrand = area - 2*intersection
integral_val = sp.integrate(integrand, (y,0,x), (x,0,sp.Rational(1,2))
result = integral_val*8
This outputs -17*log(2)/6 + 1/12 + pi/6 + 4*log(4)/3
We can simplify this further:
simplified_log = sp.simplify(-17 * sp.log(2) / 6 + 4*sp.log(4)/3)
# -log(2)/6
Therefore the final answer is: $$ \frac{1+2\pi - \ln(4)}{12} $$
Bonus: An easy mistake when using sympy#
When I first wrote the previous code, I used 0.5 instead of sp.Rational(1,2) as the bound
for the integral:
integral_val = sp.integrate(integrand, (y,0,x), (x,0,0.5))
This is a mistake because it prevented sympy from giving an exact answer, however it still tried, leading to a potentially misleading output that mixes exact terms with numerical approximations.
result # -2*log(2)/3 + 0.125*pi + 0.560806617512881
However, you can still use evaluate the output to get the correct numerical approximation:
result.evalf()*8 # 0.491407578838308