[proofplan]
We compute the first variation of the quadratic functional along an arbitrary admissible direction $v \in H^1_0(U)$. The expansion of $I[u+tv]$ is a quadratic polynomial in the real parameter $t$, whose linear coefficient is exactly the difference between the bilinear Dirichlet term and the forcing term. Therefore vanishing of every directional derivative is equivalent to the displayed weak Euler-Lagrange identity for every test direction.
[/proofplan]
[step:Verify that the directional energy is a finite quadratic polynomial]
Fix $u \in H^1_0(U)$ and $v \in H^1_0(U)$. Define the directional energy map $\Phi_v: \mathbb{R} \to \mathbb{R}$ by $\Phi_v(t) := I[u + tv]$ for $t \in \mathbb{R}$.
Since $u,v \in H^1_0(U)$, we have $u,v \in L^2(U)$ and $\nabla u,\nabla v \in L^2(U;\mathbb{R}^n)$. Since $a \in L^\infty(U)$, the functions $a|\nabla u|^2$ and $a|\nabla v|^2$ belong to $L^1(U)$. For the mixed term, the pointwise estimate
\begin{align*}
|a(x)\nabla u(x)\cdot\nabla v(x)| \leq \|a\|_{L^\infty(U)}|\nabla u(x)||\nabla v(x)|
\end{align*}
holds for $\mathcal{L}^n$-a.e. $x \in U$, and the [Cauchy-Schwarz Inequality](/theorems/432) in $L^2(U;\mathbb{R}^n)$ gives $|\nabla u||\nabla v| \in L^1(U)$. Hence $a\,\nabla u\cdot\nabla v \in L^1(U)$. Since $f,v,u \in L^2(U)$, the [Cauchy-Schwarz Inequality](/theorems/432) in $L^2(U)$ also gives $fu,fv \in L^1(U)$. Hence every integral below is finite.
Using linearity of the weak gradient,
\begin{align*}
\nabla(u+tv)(x) = \nabla u(x) + t\nabla v(x)
\end{align*}
for $\mathcal{L}^n$-a.e. $x \in U$. Expanding the Euclidean square gives
\begin{align*}
|\nabla u(x) + t\nabla v(x)|^2 = |\nabla u(x)|^2 + 2t\nabla u(x)\cdot\nabla v(x) + t^2|\nabla v(x)|^2
\end{align*}
for $\mathcal{L}^n$-a.e. $x \in U$. Substituting this identity into the definition of $I$ and using linearity of the [Lebesgue integral](/page/Lebesgue%20Integral) gives
\begin{align*}
\Phi_v(t) = I[u] + t\left(\int_U a(x)\nabla u(x)\cdot\nabla v(x)\,d\mathcal{L}^n(x) - \int_U f(x)v(x)\,d\mathcal{L}^n(x)\right) + \frac{t^2}{2}\int_U a(x)|\nabla v(x)|^2\,d\mathcal{L}^n(x)
\end{align*}
for every $t \in \mathbb{R}$.
[guided]
Fix $u \in H^1_0(U)$ and choose an arbitrary admissible direction $v \in H^1_0(U)$. The path through $u$ in direction $v$ is the map $\Phi_v: \mathbb{R} \to \mathbb{R}$ defined by $\Phi_v(t) := I[u + tv]$ for $t \in \mathbb{R}$.
Before expanding, we check that the expressions are legitimate. Since $u,v \in H^1_0(U)$, both functions belong to $L^2(U)$ and their weak gradients $\nabla u,\nabla v$ belong to $L^2(U;\mathbb{R}^n)$. The coefficient $a$ is essentially bounded, so multiplying an $L^1$ function by $a$ preserves integrability. In particular, $a|\nabla u|^2$ and $a|\nabla v|^2$ are integrable over $U$. For the mixed term, we use the pointwise Euclidean estimate
\begin{align*}
|a(x)\nabla u(x)\cdot\nabla v(x)| \leq \|a\|_{L^\infty(U)}|\nabla u(x)||\nabla v(x)|
\end{align*}
for $\mathcal{L}^n$-a.e. $x \in U$. The [Cauchy-Schwarz Inequality](/theorems/432) in the [Hilbert space](/page/Hilbert%20Space) $L^2(U;\mathbb{R}^n)$ gives $|\nabla u||\nabla v| \in L^1(U)$, so $a\,\nabla u\cdot\nabla v \in L^1(U)$. Also, $f \in L^2(U)$ and $u,v \in L^2(U)$, so the [Cauchy-Schwarz Inequality](/theorems/432) in $L^2(U)$ gives $fu,fv \in L^1(U)$. Thus all terms in the expansion of $I[u+tv]$ are finite.
The weak gradient is linear, so for $\mathcal{L}^n$-a.e. $x \in U$,
\begin{align*}
\nabla(u+tv)(x) = \nabla u(x) + t\nabla v(x)
\end{align*}
Expanding the Euclidean norm square pointwise gives
\begin{align*}
|\nabla u(x) + t\nabla v(x)|^2 = |\nabla u(x)|^2 + 2t\nabla u(x)\cdot\nabla v(x) + t^2|\nabla v(x)|^2
\end{align*}
for $\mathcal{L}^n$-a.e. $x \in U$. Substituting this into the definition of $I$ separates the constant, linear, and quadratic powers of $t$:
\begin{align*}
\Phi_v(t) = I[u] + t\left(\int_U a(x)\nabla u(x)\cdot\nabla v(x)\,d\mathcal{L}^n(x) - \int_U f(x)v(x)\,d\mathcal{L}^n(x)\right) + \frac{t^2}{2}\int_U a(x)|\nabla v(x)|^2\,d\mathcal{L}^n(x)
\end{align*}
This is the central computation: the first variation is the coefficient of $t$ in this quadratic polynomial.
[/guided]
[/step]
[step:Identify the derivative at the origin]
Since $\Phi_v$ is a real quadratic polynomial in $t$, it is differentiable at $t=0$, and the formula above gives
\begin{align*}
\Phi_v'(0) = \int_U a(x)\nabla u(x)\cdot\nabla v(x)\,d\mathcal{L}^n(x) - \int_U f(x)v(x)\,d\mathcal{L}^n(x)
\end{align*}
for every $v \in H^1_0(U)$.
[/step]
[step:Equate weak criticality with the variational identity]
By definition, $u$ is a weak critical point of $I$ if and only if $\Phi_v'(0)=0$ for every $v \in H^1_0(U)$. Using the derivative formula just obtained, this condition is equivalent to
\begin{align*}
\int_U a(x)\nabla u(x)\cdot\nabla v(x)\,d\mathcal{L}^n(x) - \int_U f(x)v(x)\,d\mathcal{L}^n(x) = 0
\end{align*}
for every $v \in H^1_0(U)$. Rearranging the equality gives
\begin{align*}
\int_U a(x)\nabla u(x)\cdot\nabla v(x)\,d\mathcal{L}^n(x) = \int_U f(x)v(x)\,d\mathcal{L}^n(x)
\end{align*}
for every $v \in H^1_0(U)$. This is precisely the asserted weak Euler-Lagrange equation, and the equivalence follows in both directions.
[/step]
