Problem 525: Rolling Ellipse

Project Euler Problem 525 Solution

EulerSolve provides an optimized solution for Project Euler Problem 525, Rolling Ellipse, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary An ellipse with semiaxes $a$ and $b$ rolls on a straight line without slipping. Let $C(a,b)$ denote the arc length of the trajectory traced by its center during one complete revolution. The problem asks for $$C(1,4)+C(3,4).$$ The implementations do not search for an elementary antiderivative. Instead, they derive an exact speed formula for the center, reduce the geometry to a smooth one-dimensional integral, and evaluate that integral numerically to high precision. Mathematical Approach Work in the ellipse's own coordinate system and parameterize the boundary point that is instantaneously in contact with the ground by $$r(t)=(x(t),y(t))=(a\cos t,b\sin t),\qquad 0\le t\le 2\pi.$$ As $t$ runs once around the ellipse, the contact point visits the whole boundary once, so one full revolution of the rolling motion is captured by integrating over one full period of $t$. Step 1: Parameterize the ellipse and its local geometry The first and second derivatives are $$r'(t)=(-a\sin t,b\cos t),\qquad r''(t)=(-a\cos t,-b\sin t).$$ The distance from the center of the ellipse to the contact point is $$\lVert r(t)\rVert=\sqrt{a^2\cos^2 t+b^2\sin^2 t}.$$ This distance varies with $t$ unless $a=b$, which is exactly why the center of a rolling ellipse does not move as simply as the center of a rolling circle....

Detailed mathematical approach

Problem Summary

An ellipse with semiaxes $a$ and $b$ rolls on a straight line without slipping. Let $C(a,b)$ denote the arc length of the trajectory traced by its center during one complete revolution. The problem asks for

$$C(1,4)+C(3,4).$$

The implementations do not search for an elementary antiderivative. Instead, they derive an exact speed formula for the center, reduce the geometry to a smooth one-dimensional integral, and evaluate that integral numerically to high precision.

Mathematical Approach

Work in the ellipse's own coordinate system and parameterize the boundary point that is instantaneously in contact with the ground by

$$r(t)=(x(t),y(t))=(a\cos t,b\sin t),\qquad 0\le t\le 2\pi.$$

As $t$ runs once around the ellipse, the contact point visits the whole boundary once, so one full revolution of the rolling motion is captured by integrating over one full period of $t$.

Step 1: Parameterize the ellipse and its local geometry

The first and second derivatives are

$$r'(t)=(-a\sin t,b\cos t),\qquad r''(t)=(-a\cos t,-b\sin t).$$

The distance from the center of the ellipse to the contact point is

$$\lVert r(t)\rVert=\sqrt{a^2\cos^2 t+b^2\sin^2 t}.$$

This distance varies with $t$ unless $a=b$, which is exactly why the center of a rolling ellipse does not move as simply as the center of a rolling circle.

Step 2: Turn rolling without slipping into a speed law

For pure rolling, the contact point has zero velocity relative to the ground at the instant of contact. Therefore that point is the instantaneous center of rotation of the rigid body.

If $\theta(t)$ is the current orientation angle of the ellipse, then the center moves around the contact point with speed

$$v(t)=\left|\frac{d\theta}{dt}\right|\lVert r(t)\rVert.$$

So the geometric problem is reduced to finding the angular speed $\theta'(t)$ associated with the changing tangent direction of the ellipse.

Step 3: Differentiate the tangent angle

For a smooth planar parametrized curve, the derivative of the tangent angle is

$$\frac{d\theta}{dt}=\frac{x'(t)y''(t)-y'(t)x''(t)}{x'(t)^2+y'(t)^2}.$$

Substituting the ellipse derivatives gives

$$x'(t)y''(t)-y'(t)x''(t)=ab,$$

$$x'(t)^2+y'(t)^2=a^2\sin^2 t+b^2\cos^2 t,$$

hence

$$\frac{d\theta}{dt}=\frac{ab}{a^2\sin^2 t+b^2\cos^2 t}.$$

The denominator is always positive for $a,b>0$, so the angular speed stays positive and no sign ambiguity remains in the arc-length integral.

Step 4: Obtain the closed integrand for the center speed

Combining the previous two steps yields

$$v(t)=\frac{ab\sqrt{a^2\cos^2 t+b^2\sin^2 t}}{a^2\sin^2 t+b^2\cos^2 t}.$$

Therefore the total center-path length over one revolution is

$$C(a,b)=\int_0^{2\pi}\frac{ab\sqrt{a^2\cos^2 t+b^2\sin^2 t}}{a^2\sin^2 t+b^2\cos^2 t}\,dt.$$

A useful sanity check is the circular case $a=b=R$. Then the integrand becomes the constant $R$, so

$$C(R,R)=\int_0^{2\pi}R\,dt=2\pi R,$$

which is exactly the expected path length for the center of a rolling circle.

Step 5: Use symmetry to integrate only one quarter-turn

The integrand depends only on $\sin^2 t$ and $\cos^2 t$. Consequently

$$v(t+\pi)=v(t),\qquad v(\pi-t)=v(t).$$

These symmetries imply

$$C(a,b)=4\int_0^{\pi/2}\frac{ab\sqrt{a^2\cos^2 t+b^2\sin^2 t}}{a^2\sin^2 t+b^2\cos^2 t}\,dt.$$

This is the exact form evaluated by the implementation, and it is numerically convenient because the interval is short and the integrand is smooth on $[0,\pi/2]$.

Step 6: Worked Example for $C(2,4)$

For $a=2$ and $b=4$, the quarter-turn integrand becomes

$$v(t)=\frac{4\sqrt{\cos^2 t+4\sin^2 t}}{\sin^2 t+4\cos^2 t}.$$

At the three standard Simpson nodes we get

$$v(0)=1,\qquad v\left(\frac{\pi}{4}\right)=\frac{4\sqrt{10}}{5},\qquad v\left(\frac{\pi}{2}\right)=8.$$

A single Simpson panel on $[0,\pi/2]$ gives the rough estimate

$$S=\frac{\pi/2}{6}\left(1+4\cdot\frac{4\sqrt{10}}{5}+8\right)=\frac{\pi}{12}\left(9+\frac{16\sqrt{10}}{5}\right)\approx 5.0056.$$

Multiplying by $4$ would give about $20.0224$, which is not yet accurate enough. After adaptive refinement the numerical integral converges to

$$C(2,4)\approx 21.38816906,$$

showing why an adaptive method is used instead of a single coarse panel.

How the Code Works

The C++, Python, and Java implementations all follow the same algorithm. First they evaluate the closed-form speed formula above. Next they apply Simpson's rule on the quarter interval $[0,\pi/2]$ using the endpoint values and the midpoint value.

They then recurse adaptively. For a current interval, the implementation evaluates the speed at the two quarter points, compares the sum of the left and right Simpson estimates with the Simpson estimate on the whole interval, and interprets the difference as an error signal.

If the estimated local error already satisfies

$$|\Delta|\le 15\varepsilon,$$

or if the recursion limit has been reached, the code accepts the interval and applies the standard Richardson correction

$$S_{\text{refined}}=S_{\text{left}}+S_{\text{right}}+\frac{\Delta}{15}.$$

Otherwise it splits the interval in half and repeats the same process on both subintervals. The requested integral is obtained by multiplying the converged quarter-interval value by $4$. Finally the implementation evaluates $C(1,4)$ and $C(3,4)$, adds them, and prints the result to eight decimal places. The numerical parameters are $\varepsilon=10^{-13}$ and a maximum recursion depth of $30$.

Complexity Analysis

Adaptive Simpson integration has no fixed iteration count in advance, so it is best analyzed in terms of the number of accepted subintervals. If the recursion finishes with $N$ accepted leaves, the total time is $O(N)$ function evaluations up to constant factors, because each subdivision adds only a constant number of new samples. The memory usage is $O(D)$, where $D$ is the recursion depth; in the implementation $D\le 30$.

For this problem the integrand is smooth on $[0,\pi/2]$, so $N$ stays modest in practice. The adaptive strategy is efficient because it refines only where the shape of the integrand demands it, instead of spending the same resolution everywhere.

Footnotes and References

Problem page: https://projecteuler.net/problem=525
Ellipse: Wikipedia — Ellipse
Arc length: Wikipedia — Arc length
Adaptive Simpson's method: Wikipedia — Adaptive Simpson's method

Problem 525 source code

C++

#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>

namespace {

long double integrand_center_speed(const long double a, const long double b, const long double t) {
    const long double s = std::sin(t);
    const long double c = std::cos(t);
    const long double num = a * b * std::sqrt(a * a * c * c + b * b * s * s);
    const long double den = a * a * s * s + b * b * c * c;
    return num / den;
}

long double simpson(const long double fa, const long double fm, const long double fb, const long double a,
                    const long double b) {
    return (b - a) * (fa + 4 * fm + fb) / 6;
}

long double adaptive_simpson(long double (*f)(long double, long double, long double), const long double a_param,
                             const long double b_param, const long double a, const long double b,
                             const long double eps, const long double whole, const long double fa,
                             const long double fm, const long double fb, int depth) {
    const long double m = (a + b) / 2;
    const long double l = (a + m) / 2;
    const long double r = (m + b) / 2;

    const long double fl = f(a_param, b_param, l);
    const long double fr = f(a_param, b_param, r);

    const long double left = simpson(fa, fl, fm, a, m);
    const long double right = simpson(fm, fr, fb, m, b);
    const long double delta = left + right - whole;

    if (depth <= 0 || std::fabsl(delta) <= 15 * eps) {
        // Richardson extrapolation.
        return left + right + delta / 15;
    }
    return adaptive_simpson(f, a_param, b_param, a, m, eps / 2, left, fa, fl, fm, depth - 1) +
           adaptive_simpson(f, a_param, b_param, m, b, eps / 2, right, fm, fr, fb, depth - 1);
}

long double integrate_center_curve_length(const long double a, const long double b) {
    // From rolling kinematics: if r(t)=(a cos t, b sin t) is the contact point relative to the center,
    // and φ(t) is the tangent angle, then the center velocity is |C'(t)| = |φ'(t)| * |r(t)|.
    //
    // With r'(t)=(-a sin t, b cos t), we get:
    //   φ'(t) = (x y' - y x') / |r'(t)|^2 = ab / (a^2 sin^2 t + b^2 cos^2 t),
    //   |r(t)| = sqrt(a^2 cos^2 t + b^2 sin^2 t).
    //
    // Hence the arc-length integrand is:
    //   ab * sqrt(a^2 cos^2 t + b^2 sin^2 t) / (a^2 sin^2 t + b^2 cos^2 t).
    //
    // The integrand is π-periodic and symmetric, so one full turn is:
    //   C(a,b) = 4 * ∫[0, π/2] integrand(t) dt.
    const long double pi = acosl(-1.0L);
    const long double A = 0.0L;
    const long double B = pi / 2;
    const long double fa = integrand_center_speed(a, b, A);
    const long double fb = integrand_center_speed(a, b, B);
    const long double m = (A + B) / 2;
    const long double fm = integrand_center_speed(a, b, m);
    const long double whole = simpson(fa, fm, fb, A, B);
    const long double eps = 1e-13L;
    const long double quarter = adaptive_simpson(&integrand_center_speed, a, b, A, B, eps, whole, fa, fm, fb, 30);
    return 4 * quarter;
}

bool run_checkpoints() {
    const long double v = integrate_center_curve_length(2.0L, 4.0L);
    const long double expected = 21.38816906L;
    if (std::fabsl(v - expected) > 5e-9L) {
        std::cerr << std::setprecision(15) << "Checkpoint failed: C(2,4) got " << v << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main() {
    if (!run_checkpoints()) {
        return 1;
    }

    const long double c14 = integrate_center_curve_length(1.0L, 4.0L);
    const long double c34 = integrate_center_curve_length(3.0L, 4.0L);
    const long double ans = c14 + c34;
    std::cout << std::fixed << std::setprecision(8) << ans << '\n';
    return 0;
}

Python

import math

def integrand_center_speed(a, b, t):
    s = math.sin(t)
    c = math.cos(t)
    num = a * b * math.sqrt(a * a * c * c + b * b * s * s)
    den = a * a * s * s + b * b * c * c
    return num / den

def simpson(fa, fm, fb, a, b):
    return (b - a) * (fa + 4 * fm + fb) / 6

def adaptive_simpson(a_param, b_param, a, b, eps, whole, fa, fm, fb, depth):
    m = (a + b) / 2
    l = (a + m) / 2
    r = (m + b) / 2
    
    fl = integrand_center_speed(a_param, b_param, l)
    fr = integrand_center_speed(a_param, b_param, r)
    
    left = simpson(fa, fl, fm, a, m)
    right = simpson(fm, fr, fb, m, b)
    delta = left + right - whole
    
    if depth <= 0 or abs(delta) <= 15 * eps:
        return left + right + delta / 15
        
    return adaptive_simpson(a_param, b_param, a, m, eps / 2, left, fa, fl, fm, depth - 1) + \
           adaptive_simpson(a_param, b_param, m, b, eps / 2, right, fm, fr, fb, depth - 1)

def integrate_center_curve_length(a, b):
    A = 0.0
    B = math.pi / 2
    fa = integrand_center_speed(a, b, A)
    fb = integrand_center_speed(a, b, B)
    m = (A + B) / 2
    fm = integrand_center_speed(a, b, m)
    whole = simpson(fa, fm, fb, A, B)
    eps = 1e-13
    quarter = adaptive_simpson(a, b, A, B, eps, whole, fa, fm, fb, 30)
    return 4 * quarter

def solve():
    c14 = integrate_center_curve_length(1.0, 4.0)
    c34 = integrate_center_curve_length(3.0, 4.0)
    ans = c14 + c34
    return f"{ans:.8f}"

if __name__ == '__main__':
    print(solve())

Java

public class Euler525 {

    static double integrandCenterSpeed(double a, double b, double t) {
        double s = Math.sin(t);
        double c = Math.cos(t);
        double num = a * b * Math.sqrt(a * a * c * c + b * b * s * s);
        double den = a * a * s * s + b * b * c * c;
        return num / den;
    }

    static double simpson(double fa, double fm, double fb, double a, double b) {
        return (b - a) * (fa + 4 * fm + fb) / 6;
    }

    static double adaptiveSimpson(double aParam, double bParam, double a, double b, double eps, double whole,
            double fa, double fm, double fb, int depth) {
        double m = (a + b) / 2;
        double l = (a + m) / 2;
        double r = (m + b) / 2;

        double fl = integrandCenterSpeed(aParam, bParam, l);
        double fr = integrandCenterSpeed(aParam, bParam, r);

        double left = simpson(fa, fl, fm, a, m);
        double right = simpson(fm, fr, fb, m, b);
        double delta = left + right - whole;

        if (depth <= 0 || Math.abs(delta) <= 15 * eps) {
            return left + right + delta / 15;
        }

        return adaptiveSimpson(aParam, bParam, a, m, eps / 2, left, fa, fl, fm, depth - 1) +
                adaptiveSimpson(aParam, bParam, m, b, eps / 2, right, fm, fr, fb, depth - 1);
    }

    static double integrateCenterCurveLength(double a, double b) {
        double A = 0.0;
        double B = Math.PI / 2;
        double fa = integrandCenterSpeed(a, b, A);
        double fb = integrandCenterSpeed(a, b, B);
        double m = (A + B) / 2;
        double fm = integrandCenterSpeed(a, b, m);
        double whole = simpson(fa, fm, fb, A, B);
        double eps = 1e-13;
        double quarter = adaptiveSimpson(a, b, A, B, eps, whole, fa, fm, fb, 30);
        return 4 * quarter;
    }

    public static void main(String[] args) {
        double c14 = integrateCenterCurveLength(1.0, 4.0);
        double c34 = integrateCenterCurveLength(3.0, 4.0);
        double ans = c14 + c34;
        System.out.printf(java.util.Locale.US, "%.8f\n", ans);
    }
}