grazeint.py

import math
from scipy.optimize import brentq
import matplotlib.pyplot as plt

# Solve the 10/12/18 Riddler Classic problem on fivethirtyeight.com.
# The following is the result of several steps...
#
# 1. WLOG, we take the field to be a unit circle, and arrange the coordinate
#    system so the stake is at (1,0).
# 2. The goat-grazeable region is bounded on the left by the arc corresponding
#    to what the goat can reach and on the right by the field.
# 3. Given the coordinate system, the x axis divides the goat-grazeable region
#    precisely in half.
# 4. The two circles cross (above the x axis) where sqrt(r^2 - (x - 1)^2) =
#    sqrt(1 - x^2), i.e. at x = (2 - r^2)/2. Let's call that the "zenith".
# 5. The area of the top half of the goat-grazeable region is thus the sum of
#    @  the integral of sqrt(r^2 - (x - 1)^2) dx from 1-r to the zenith
#    @  the integral of sqrt(1 - x^2) from the zenith to 1
# 6. Thanks to https://www.integral-calculator.com/, we know that the
#    antiderivative of sqrt(1 - x^2) is (arcsin(x) + x*sqrt(1 - x^2))/2
#    (we ignore the constant term because it will subtract out)
# 7. The antiderivative of sqrt(r^2 - (x - 1)^2) was hairy, so we reduced
#    it to the unit circle problem by moving the center back to the origin,
#    scaling the coordinates down by a factor of r, doing the unit circle
#    integral, and multiplying the result by r^2. (Whew.)
# 8. That gets us out of the numerical integration business, but not to the
#    point that there's a closed-form solution, so thank goodness for
#    scipy.optimize.brentq.

# For a given choice of r, return the difference between the goat-grazeable
# region and half the field.
def delta(r):
    zenith = (2 - r*r) / 2
    antideriv = lambda x: (math.asin(x) + x*math.sqrt(1 - x*x))/2
    lhs = r*r * (antideriv((zenith - 1)/r) - antideriv(-1))
    rhs = antideriv(1) - antideriv(zenith)
    return 2 * (lhs + rhs) - math.pi/2

# Find r that works for R = 1. We know 1 < r < 2, so...oops. delta() assumes
# the circles intersect above the x axis, so we'll go slightly below 2.
def goatRopeLength():
    return brentq(delta, 1, 1.8)

# Plot the solution.
def plotsoln():
    r = goatRopeLength()
    fig, ax = plt.subplots()
    ax.add_artist(plt.Circle((0, 0), 1, color='#f0000060'))
    ax.add_artist(plt.Circle((1, 0), r, color='#0000f060'))
    plt.grid(True)
    plt.xlim(min(-1, 1-r) -0.25, 1+r + 0.25)
    plt.ylim(-r - 0.25, r + 0.25)
    plt.title('solution to 10/12/18 Riddler Classic: r={0:.5f}'.format(r))
    plt.show()

Added plot routine. Clearly I need to make my way through matplotlib in an organized fashion, or cut to the chase and do the same with plotnine.

Got out of the numerical integration business, which should get rid of a source of error. Also reduced the problem of the area of the intersection bounded by the radius r circle to the unit circle, getting r out of the integration as much as possible.