converting $(a(x)u(x))_{xx}=f(x)$ to weak form in FEM

Question

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.

Answer 1

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.