Solving nonlinear equation resulting from finite element method

Question

Using the finite element method (for a uniform mesh in the spatial domain) I have the system with initial conditions $u_j(0)=\cos(x_j)$ for $j=1,\dots, N$ $$\frac{d\vec{u}}{dt}=A\vec{u}+B\vec{c},$$ where $A$ and $B$ are coefficient matrices and $\vec{c}$ takes the form $$\vec{c}=\begin{bmatrix}\displaystyle\int_0^1\left(\vec{u}\cdot\beta(x)\right)^2\beta_j'(x)dx\end{bmatrix}_{1\leq j\leq N}$$ and $u_j=u_j(t)\approx u(x_j,t)$ is the semi-approximation of the solution at time $t$ on the mesh node $x_j$ for $j=1,\dots, N$ and $\beta_j$ is a linear basis function. My book says that this system can be solved with a "standard ODE solver", but provides no other insight besides that. My question is how does one go about solving such a system? Is there a particular MATLAB command that is used to solve this system? After searching online I haven't been able to find something conrete. Any discussion or insight would be greatly appreciated.

Answer 1

After doing the FEM discretization in space, you obtain a semi-discretized problem: Your temporal derivative $\frac{\mathrm d u }{\mathrm dt}$ remained, but you managed to approximate the spatial derivatives. The important thing now to realize is that your coefficients of the FEM functions are functions of time, i.e., $\vec{u} = \vec{u}(t)$. To get a standard non-linear ODE, you want to get rid of the integral. You could do this in the standard procedure by transforming it to the reference element and then use a quadrature rule to transform the integral essentially into a sum. You could also simplify the integral a bit a priori by noting that $$\big(\vec{u}\cdot\beta(x)\big)^2 = \Bigg(\sum_k^N u_k \beta_k(x) \Bigg) \cdot\Bigg(\sum_i^N u_i \beta_i(x) \Bigg) = \sum_{i,k}^N u_ku_i\beta_k(x)\beta_i(x)$$ and thus $$c_j = \int_0^1 \sum_{i,k}^N \beta_j'(x) u_ku_i\beta_k(x)\beta_i(x) \mathrm d x = \sum_{i,j}^N u_ju_i \int_0^1 \beta_j'(x) \beta_k(x)\beta_i(x) \mathrm d x$$ This potentially allows you avoiding the quadrature since you know your basis functions $\vec{\beta}(x)$, i.e., you can precompute the integral.

In both cases, however, you end up with a problem like $$\frac{\mathrm d \vec{u}}{\mathrm d t} = \vec{F}\big(\vec{u}(t)\big), $$ a nonlinear ODE for the coefficients $\vec{u}(t)$. The standard command for this in Matlab would be ode45, although there are some other solvers for stiff systems, ...