Rectangle of maximum area under the curve $y=x(x-1)^2$ A rectangle's bottom is $y=0$. While the top corners are on the curve $y=x(x-1)^2$ between $x=0$ and $x=1$. Find the maximum area of this rectangle.
My Progress
Defining $f(x) =x(x-1) ^2$
I first starting by assuming that there exist $a$ and $b$ such that $0<a, b<1$ and $f(a) =f(b)$.
The area of this rectangle is given by $|b-a|f(a) $.
We know that $a(a-1) ^2=b(b-1) ^2$. Neglecting the $a=b$ case we get that $$a^2+b^2+ab-2(a+b) +1=0$$
Now it remains to maximise $$|b-a|a(a-1) ^2$$ under the constraint $$a^2+b^2+ab-2(a+b) +1=0$$
But now I'm pretty much stuck here.
 A: There is really only one free variable, which is without loss of generality, the $x$ coordinate of the vertex of the rectangle that is closest to the origin.  If we call this value $a$, then the other vertex coordinates are uniquely defined:  in total, we have
$$(a,0), (a, a(1-a)^2), (b,a(1-a)^2), (b, 0)$$  where $$b = \frac{2-a - \sqrt{4a - 3a^2}}{2}$$ is the solution to $f(a) = f(b)$ for which $b > a$.  The resulting area of the rectangle is
$$R(a) = a(1-a)^2 \cdot (b-a) = \frac{1}{2}\left( a(1-a)^2 \left(2-3a - \sqrt{4a-3a^2}\right)\right).$$  A tedious differentiation with respect to $a$ yields
$$R'(a) = 1 - 7a + 12a^2 - 6a^3 - (3 - 10a + 6a^2)(1-a) \sqrt{\frac{a}{4-3a}}.$$  So any critical points must satisfy $R'(a) = 0$, or $$(1-7a+12a^2-6a^3)^2(4-3a) = (3-10a+6a^2)^2(1-a)^2 a.$$  This eventually leads to the condition
$$0 = (1-a)^3 (1 - 14a + 58a^2 - 84a^3 + 36a^4).$$  Since $a < 1$, this latter factor yields the minimal polynomial $f(z) = 1 - 14z + 58z^2 - 84z^3 + 36z^4$ for $a$, which actually admits a solution by square roots.  We find this by presuming that $f$ admits a factorization into quadratics of the form
$$f(z) = (1 + pz + 6z^2)(1 + qz + 6z^2)$$
for undetermined constants $p, q$.  Then expanding and equating like coefficients yields the conditions $$p+q = -14, \quad pq + 12 = 58.$$  This gives $(p,q) = -7 \pm \sqrt{3}$ in any order, hence the roots of $f$ are
$$z \in \left\{ \frac{-p \pm \sqrt{p^2 - 24}}{12}, \quad \frac{-q \pm \sqrt{q^2 - 24}}{12} \right\}.$$
Numeric evaluation of the four roots shows that only one satisfies $0 < a < b < 1$, namely $$a = \frac{7 + \sqrt{3} - \sqrt{14(2+\sqrt{3})}}{12} = \frac{7 + \sqrt{3} - \sqrt{7} - \sqrt{21}}{12} \approx 0.12531031679570389137.$$
The corresponding maximum area is then $$R(a) = \frac{7 \sqrt{7} - 10}{108 \sqrt{3}} \approx 0.045547906756179342935.$$
By no means do I claim that this is the most elegant or simplest solution; it is simply what I was able to compute on a first try at this question.
A: This problem is ideally suited for optimization by method of Lagrange multipliers due to its multidimensional nature with boundary curve--especially since the boundary curve is quite difficult to solve explicitly or parametrize. Thus my solution will use this method.
As with one-dimensional optimization, you have an objective function $f(x)$, and  this function is minimized/maximized when $f'(x)=0$. Similarly, when you have a parametrized constraint curve $\vec c(t)$ in $N$-D space, the objective function $f(\vec x)$ is maximized when $f'(\vec c(t))=0$, that is, when the derivative along the curve is zero. Now for an objective function in $N$-D space, the derivative in a direction is zero iff that direction is perpendicular to the gradient $\nabla f(\vec x)$. So we need to make $\nabla f(\vec x)$ perpendicular to the constraint curve.
Here's where the trick comes in. Instead of trying to explicitly parametrize the constraint, leave it implicit by writing it as a level curve of a function $g(\vec x)$. That is, your constraint is the equation $g(\vec x)=c$. Since this function is constant everywhere along the constraint, its derivative in the direction of the constraint is zero--meaning $\nabla g(\vec x)$ is perpendicular to the constraint curve. Since $\nabla f(\vec x)$ and $\nabla g(\vec x)$ are both perpendicular to the constraint, they must be parallel to each other: $\nabla f(\vec x)=\lambda\nabla g(\vec x)$.

Solving the problem at hand, we have the objective function $f(a,b)=(b-a)a(1-a)^2$ and the constraint $g(a,b) = a^2+b^2+ab-2(a+b)+1 = 0$.
Finding the gradients:
$$\nabla f(a,b)=\left\langle\frac{df}{da},\frac{df}{db}\right\rangle=\left\langle-4a^3+3a^2b+6a^2-4ab-2a+b,a^3-2a^2+a\right\rangle$$
$$\nabla g(a,b)=\left\langle\frac{dg}{da},\frac{dg}{db}\right\rangle=\left\langle2a+b-2,a+2b-2\right\rangle$$
And now we have our system of equations:
$$\left\{\begin{align}
a^2+b^2+ab-2a-2b+1 &= 0\\
-4a^3+3a^2b+6a^2-4ab-2a+b &= \lambda(2a+b-2) \\
a^3-2a^2+a &= \lambda(a+2b-2)
\end{align}\right.$$
This system is difficult to solve by hand (combining them gives a 4th order polynomial), but its solution can be found numerically to be $a=0.1253,b=0.6004$.
