Solution

Answer: 23.386029052

Let the silo have radius $R$, and let the grain be poured from the point $(x,0)$ on the top of the cylinder.

If the angle of repose is $\alpha$, then the free surface of the grain makes slope

$$\tan\alpha$$

with the horizontal.

For any point $(u,v)$ in the circular cross-section, the vertical empty height above the grain is therefore

$$h(u,v)=\tan\alpha;\sqrt{(u-x)^2+v^2}.$$

Hence the wasted volume is

$$V(x)=\tan\alpha\iint_D \sqrt{(u-x)^2+v^2},dA,$$

where $D$ is the disk $u^2+v^2\le R^2$.

1. Converting the integral

Use polar coordinates centered at the pouring point:

$$u=x+r\cos\theta,\qquad v=r\sin\theta.$$

Along direction $\theta$, the ray exits the silo when

$$(x+r\cos\theta)^2+(r\sin\theta)^2=R^2.$$

Solving for $r$,

$$r_{\max}(\theta) = -x\cos\theta+\sqrt{R^2-x^2\sin^2\theta}.$$

Since the integrand is simply $r$,

$$V(x) = \tan\alpha \int_0^{2\pi}\int_0^{r_{\max}(\theta)} r\cdot r,dr,d\theta.$$

Therefore

$$\boxed{ V(x)= \frac{\tan\alpha}{3} \int_0^{2\pi} \left( -x\cos\theta+\sqrt{R^2-x^2\sin^2\theta} \right)^3 d\theta }.$$

2. Checking the example

The statement’s first example uses:

diameter $6\Rightarrow R=3$,
$\alpha=30^\circ$,
$x=0$.

Then

$$V(0) = \frac{\tan30^\circ}{3} \int_0^{2\pi} 3^3,d\theta = 18\pi\tan30^\circ = 6\pi\sqrt3 \approx 32.648388556,$$

which matches the problem statement exactly.

So the formula is correct.

3. Actual problem data

Now use:

$$R=6,\qquad \alpha=40^\circ.$$

We seek all $x\in[0,6]$ such that $V(x)$ is a perfect square.

Since $V(x)$ is strictly increasing with $x$, there is at most one solution for each square.

Compute the range:

$$V(0)\approx 379.60, \qquad V(6)\approx 644.43.$$

Thus the only possible perfect squares are

$$20^2,21^2,22^2,23^2,24^2,25^2.$$

Numerically solving $V(x)=n^2$ for each $n=20,\dots,25$ gives:

$$\begin{aligned} x_{20}&\approx1.6097585009208716,\ x_{21}&\approx2.8060119465435142,\ x_{22}&\approx3.6783279154320040,\ x_{23}&\approx4.4264032917821366,\ x_{24}&\approx5.1094202074447333,\ x_{25}&\approx5.7561071894994169. \end{aligned}$$

Summing:

$$\sum x = 23.386029051622676\ldots$$

Rounded to $9$ decimal places:

$$\boxed{23.386029052}$$

4. Python implementation

import mpmath as mp

mp.mp.dps = 50

R = mp.mpf(6)
alpha = mp.radians(40)
T = mp.tan(alpha)

# Volume function
def V(x):
    def integrand(theta):
        return (
            -x * mp.cos(theta)
            + mp.sqrt(R*R - x*x * mp.sin(theta)**2)
        ) ** 3

    return (T / 3) * mp.quad(integrand, [0, 2 * mp.pi])

# Solve V(x) = n^2 by bisection
def solve_for_square(n):
    target = n * n

    lo = mp.mpf(0)
    hi = R

    for _ in range(100):
        mid = (lo + hi) / 2

        if V(mid) < target:
            lo = mid
        else:
            hi = mid

    return (lo + hi) / 2

roots = []

for n in range(20, 26):
    x = solve_for_square(n)
    roots.append(x)
    print(n, x)

answer = mp.fsum(roots)

print("Sum =", answer)
print("Rounded =", mp.nstr(answer, 12))

5. Code walkthrough

V(x) implements the exact integral formula derived above.
For each angle theta, the ray length inside the disk is

-x*cos(theta) + sqrt(R^2 - x^2*sin(theta)^2)

The inner radial integral contributes a factor $r^3/3$.
We numerically integrate over $[0,2\pi]$.
Since $V(x)$ is monotonic, bisection safely solves

$$V(x)=n^2.$$

Only squares from $20^2$ to $25^2$ lie in the attainable range.

6. Final verification

The answer is positive and less than $6\times6=36$, which is reasonable since each $x_i\in[0,6]$.
The computed roots are strictly increasing.
Substituting each root back into $V(x)$ reproduces the required squares to high precision.

Answer: 23.386029052

Project Euler Problem 431