0
$\begingroup$

Let's consider the usual Poisson problem for FEM $$-u''(x)=1 \quad x \in [0,1]$$ with homogeneous Dirichlet boundary conditions. The solution is $u(x)=-\frac{1}{2}x(x-1)$

I know that the FEM solution (I'm assuming linear elements) is exact at nodes, as noted also here

I implemented my code with constant step size $h$. Looking at internal nodes, the stiffness matrix has entry $$(A)_{ij}=\int \phi_i' \phi_j' dx$$ For $i =j$ (diagonal elements) I obtain $A_{ii}=2/h$, while the extradiagonal are $A_{i,i-1}=A_{i -1,i}=-1/h$ .

The load vector has entries $$(f)_i=(1,\phi_i)=h$$ In the end, I obtain the same matrix as with finite differences and hence I would expect second order. Why is exact solution at grid points recovered? (Modulo machine precision)?

$\endgroup$
0
$\begingroup$

The finite difference scheme also gives exact solution at the nodes to the problem $-u''=1$ because $$ \frac{u_{j-1} - 2 u_j + u_{j+1}}{h^2} $$ is exact for a quadratic function.

For the more general case $$ -u''(x) = f(x) $$ if you compute the integral $(f,\phi_i)$ exactly in $P_1$ FEM, it gives exact solution at the nodes. The FD scheme will not be exact in the general case.

$\endgroup$
12
  • $\begingroup$ I've just checked my implementation and I was comparing solutions at grid points shifted by one, this was the source of the error. I've seen that indeed now I obtain the exact solution at nodes with P1 FEM whenever I'm integrating exactly, while with trapezoidal rule, for instance, I obtain a second order approximation, as expected. If I use P2 FEM, should I expect again the exact value at grid points? @cfdlab $\endgroup$ – bobinthebox Mar 1 at 11:01
  • $\begingroup$ Another last question: should I expect the exact solution also for the problem $-\alpha(x) u''(x)=f(x)$ with linear elements (assuming that I integrate exactly) both the $(f,\phi_i)$ and $(\alpha(x) \phi_i', \phi_j')$ ? @cfdlab $\endgroup$ – bobinthebox Mar 1 at 16:44
  • $\begingroup$ For this particular problem you consider, if you use P2 FEM, you actually would find the exact solution $u_h=u=0.5x(1-x)$. In general, if the exact solution and the source function are sufficiently smooth, and if you evaluate the bilinear and linear forms exactly, you would expect the numerical solution and the exact solution to match at the nodes. (See interpolatory properties of continuous Galerkin FEMs) $\endgroup$ – Abdullah Ali Sivas Mar 1 at 16:45
  • $\begingroup$ @AbdullahAliSivas Regarding my last question: is it possible that I'm seeing first order of convergence with linear elements for the problem $-\alpha(x)u''(x)=f(x) $ for $\alpha,f$ smooth functions like $\sin$ and $\cos(x)$? $\endgroup$ – bobinthebox Mar 1 at 16:49
  • $\begingroup$ @bobinthebox, you need to be careful with $\alpha(x)$, it has to be positive everywhere in the domain for the problem to be well-defined. Other than that, the integration technique you use should be sufficiently high-degree and usually higher than twice the degree of the FEM (if $\sin$ and $\cos$ are involved with P1, I would probably use Simpson's method, for example) to guarantee optimal convergence rates. $\endgroup$ – Abdullah Ali Sivas Mar 1 at 16:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.