1
$\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)?

Improve this question
$\endgroup$

1 Answer 1

Reset to default
1
$\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.

Improve this answer
$\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$
    – FEGirl
    Commented Mar 1, 2021 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$
    – FEGirl
    Commented Mar 1, 2021 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
    Commented Mar 1, 2021 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$
    – FEGirl
    Commented Mar 1, 2021 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
    Commented Mar 1, 2021 at 16:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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