Consider the second order operator
$Lu=\partial_i(a_{ij}\partial_j)u+b_i\partial_iu+cu$.
Can we find a functional $I[u]$ such that $Lu$ is the variation of $I[u]$ with respect to $u$? I have successfully dealt with the first and third term, but found difficulties with the second term. Is it even possible to do so?
Actually, Math604 is closer to the right answer. To see the correct condition, you need to pay attention to the placement of your indices. Your operator should be written in the form $$ Lu = \partial_i(a^{ij}\partial_ju) + b^k\partial_ku + c u = 0, $$ where one sums over repeated indices in opposition (i.e., `one up, one down'). Also, you didn't say this, but, usually, when one writes the operator this way, one assumes that the coefficients $a^{ij}$, $b^k$ and $c$ are functions of the independent variables (say, $x^i$), and that the matrix $A = (a^{ij})$ is both symmetric and invertible. If you didn't intend for us to make these assumptions, you should say so, and we can work out the more general case.
Under these conditions, there are unique functions $f_j$ satisfying the equations $$ a^{ij}f_j = b^i. $$
Then the condition you want is simply that the $1$-form $\phi = f_i\,\mathrm{d} x^i$ should be exact, i.e., that there should exist a function $f$ on your $x$-domain such that $\partial_i f = f_i$ for all $i$. Then, the Lagrangian you want is $$ I[u] = \frac12\int_D e^{f(x)}\bigl(-a^{ij}(x)\,(\partial_iu)(\partial_ju) + c(x)\, u^2\bigr)\,\mathrm{d} x $$ If the $1$-form $\phi$ is not exact, the equation $Lu=0$ is not the Euler-Lagrange equation of any first-order functional.
I'll assume the coefficients are independent of $u\,.$ There is no standard variational problem associated with this PDE for nonzero $b_i\,.$ The reason is that for a linear PDE to admit a variational formulation, at least on $\mathbb{R}^n\,,$ it needs to be symmetric. But $\langle v,b_i\partial_i u\rangle=\langle b_iv,\partial_i u\rangle=-\langle (\partial_i b_i)v+b_i\partial_i v,u\rangle.$ So for it to be symmetric you need $-(\partial_i b_i)v-b_i\partial_i v=b_i\partial_i v\iff -v\,\partial_ib_i=2b_i\partial_i v$ for all $v\,.$ Taking $v=1$ in some neighbourhood, this gives that $\partial _i b_i=0$ there, hence everywhere since the location was arbitrary, and this implies that $b_i\partial_iv=0$ for all $v\,,$ which implies that $b_i=0\,.$
Notice this also explains why PDE's in divergence form as you've stated have a variational formulation: $\langle v,\partial_i(a_{ij}\partial_j)u\rangle=-\langle \partial_i v,a_{ij}\partial_ju\rangle=\langle \partial_j(a_{ij}\partial_i) v,u\rangle=\langle \partial_i(a_{ij}\partial_j) v,u\rangle$ if $a_{ij}$ is symmetric, and so the operator is symmetric as well.
Edit: Since this was mentioned in the other answer, maybe I should make clear what I mean by standard variational formulation: A variational formulation for which there is a functional such that $I'[u]$ is the PDE operator (which was the question), with respect to the $L^2$ inner product. Since no inner product was specified I just assumed it was $L^2$. However lots of PDEs admit variational formulations that aren't of this form that are typically said to not have variational formulations, such as the heat equation.
