converting $(a(x)u(x))_{xx}=f(x)$ to weak form in FEM Suppose I have an ODE of the form
$u(x) - [a(x)u(x)]_{xx}=f(x)$
and I want to solve it with FEM, how do I convert it to the weak form? Using integration by part, the middle term becomes
$(v,-(au)'')=(v',(au)')=(v',a'u+au')=(v',a'u)+(v',au')$
where $v=v(x)$ is the test function. If we let $v=\phi_i(x)$, and $u(x)=\sum_j u_j \phi_j(x)$, then we end up with something like
$\sum_j \left[ (\phi_i,\phi_j) + (\phi_i',a'\phi_j) + (\phi_i',a\phi_j')\right] u_j = (\phi_i,f)$
But I'm not sure if there are any subtleties I have missed. Also now looks like the matrix $A_{ij}$ is no longer symmetric, which is different from other more standard problems I have come across like the usual poisson equation $-(au')'=f$. So I would like to know how to properly convert equations like this to weak form / solving with FEM.
 A: This problem is often referred to as putting a differential equation in 'divergence form' as $b(x)u(x)-(c(x)u'(x))' = g(x)$ since this is a self-adjoint (think symmetric) operator:
$$(v,-(cu')') = (v',cu') = (cv',u') = (-(cv')',u)$$
for all sufficiently smooth $u, v$. Therefore this form is preferred whenever the weak formulation is used, such as with finite elements or when using an energy method. Your given problem $u(x) - (a(x)u(x))'' = f(x)$ is not in divergence form and the resulting operator is not self-adjoint, so we cannot expect the corresponding matrix $A_{ij}$ to be symmetric. Sometimes an equation not in divergence form can be converted to one that is, but unfortunately the procedure does not work in this case unless $a$ is constant. To see why, we calculate:
$$g = bu-(cu')' = bu - c'u' - cu''$$
$$f = u - (au)'' = (1-a'')u - 2a'u' - au''$$
and so equating coefficients of the derivatives of $u$ we must have both that $a = c$ and $2a' = c'$, which is not possible unless $a' = c' = 0$ and so $a$ is constant. Of course if this were the case then your equation would already be in divergence form by pulling $a$ outside of the $x$ derivatives.
However, the finite element method still works in the non-symmetric case as long as a weak solution is guaranteed to exist. One can show that this is true for this particular problem using the Lax-Milgram theorem but I'll omit the proof unless specifically asked. But the procedure works exactly the same as usual, just the matrix $A_{ij}$ to be inverted is no longer symmetric.
