How to handle non-linear terms on FEM for PDEs

Question

When using FEM for solving a PDE you first have to do a discretization. when we have a linear PDE is quite straightforward. You find the week form of the PDE and then make the discretization. But, what happens when you have nonlinear terms, for instance, picture the following equation $$ \frac{\partial u}{\partial t}=\nabla^2u+u-u^3 $$ How would you handle the cubic term on the discretization? You can assume some Dirichlet and Neumann boundary conditions on a domain $\Omega$. The approach I took was to assume the nonlinear term as a forcing term $$ \frac{\partial u}{\partial t}-\nabla^2u-u=F(u) $$ But because F depends explicitly on u, I don't know if that is legal to do.

Answer 1

If you are just interested in the numerical approximation, nothing really changes compared to the linear case. You multiply with a test function,

$$ v\partial_t u - v \Big( \Delta u + u - u^3 \Big) = 0,$$ integrate by parts $$ \int_\Omega v \partial_t u + \nabla v \cdot \nabla u -vu + vu^3 \mathrm d \boldsymbol{x} - \int_{\partial \Omega} (v \nabla u )\cdot \boldsymbol{n} \mathrm d \boldsymbol{s}= 0$$ choose a set of basis functions $ \phi_i(x)\: i = 1 \dots N $, substitute this into the weak form \begin{align} \int_\Omega & \phi_i(x) \partial_t \sum_{j=1}^N \widehat{u_j}(t) \phi_j(x) + \nabla \phi_i(x) \cdot \nabla \sum_{j=1}^N \widehat{u_j}(t) \phi_j(x) -\phi_i(x)\sum_{j=1}^N \widehat{u_j}(t) \phi_j(x) \\ + &\phi_i(x) \Bigg(\sum_{j=1}^N \widehat{u_j}(t) \phi_j(x) \Bigg)^3 \mathrm d \boldsymbol{x} - \int_{\partial \Omega} \phi_i(x) \nabla \sum_{j=1}^N \widehat{u_j}(t) \phi_j(x) \cdot \boldsymbol{n} \mathrm d \boldsymbol{s}= 0 \: \forall \: i = 1 \dots N \end{align} In principle, the FEM part is now already done, since you have now a nonlinear ODE in the coefficents $\widehat{u_j}(t)$ of the form $$\frac{\mathrm d }{\mathrm d t} \widehat{\boldsymbol{u}}(t) + \boldsymbol{F}\Big(\widehat{\boldsymbol{u}}(t)\Big) = 0.$$ Instead of a linear solver, however, you will now need to solve a nonlinear system and thus you would have to perform e.g. Newton iterations.

Of course, you can greatly simplify the integrals of the weak form by reducing the domain of integration for each $\phi_i \phi_j$ combination to their shared support. In particular, the support of any functions is not increased through multiplication, i.e., the occurence of the $\Big(\sum_{j=1}^N \widehat{u_i}(t) \phi_j(x) \Big)^3$ term does not make your life much worse. Consider for example piece-wise linears in 1D. Then, the only functions sharing support with $\phi_i$ are $\phi_i, \phi_{i \pm 1}$. Thus, you can reduce the "troublemaker" to \begin{align}\phi_i \Bigg(\sum_{j=i-1}^{i+1} \widehat{u_i}(t) \phi_j(x) \Bigg)^3\end{align} which can then be spelled out - it becomes an ugly, long expression, but essentially you end up with the usual scenario of time-dependent coefficients (now combined in a non-linear way) multiplied with some number depending on the integral of the base functions.