Should we always expect FEM error plots to be straight lines?

Question

The error estimates in FEM are usually of the form

$$||u^h-u||\leq Ch.$$

Taking logarithm on both sides, we obtain

$$\log ||u^h-u||\leq \log C + \log h.$$

This estimate implies that the error lies below the straight line given by $\log y=\log C + \log h$ in the log-log scale. This estimate does not say that the plot of error with respect to $h$ should be a straight line, only that it should lie below a straight line.

Then, why is it that most error plots which are published in journals show a very sharp straight line for the error? I would be grateful for some insight into this question, as I am a novice to scientific computing.

In particular, some computations that I performed on FEniCS showed error plots which are not straight lines, however the graph does lie below a straight line. Does a particular linear solver affect such behaviour? I am using mumps linear solver.

Edit: In the picture below, theoretical results show that the error (plotted on the y axis) should decay as $1/R$ where $R$ is plotted on the x axis. The second figure is in the log-log scale.

If you implement a method that converges with polynomial order N in the mesh size (and implement it correctly), that is what you should expect. — Paul, Jul 13 '19 at 01:34
Judging from the plots, it does not seem like you are refining the mesh but creating a new one with different size every time, is that correct? — nicoguaro, Jul 13 '19 at 20:39
@nicoguaro Yes, this is true. I am creating a new mesh every time. — Thangachelli Debopritama, Jul 14 '19 at 01:13
Have you tried refining the mesh instead, splitting by 4 every element every time, assuming it's a 2D problem. — nicoguaro, Jul 14 '19 at 01:25
@nicoguaro I will certainly try this, but is there a reason behind why this may cause such a behaviour? — Thangachelli Debopritama, Jul 14 '19 at 03:40
-- A direct linear solver (like MUMPS) has no impact for the level of errors you are considering (i.e. ~ $10^{-3}$) -- When the problem is coercive (i.e. $\sqrt{ a(v,v) }$ is equivalent to a norm), using a series of nested meshes will guarantee that the error decreases with $h$... but not that you get a straight line. — user7440, Jul 16 '19 at 04:34

score 3 · Accepted Answer · answered Jul 16 '19 at 09:33

While I cannot comment your particular case as you didn't provide enough details, I can give some guidance on the matter since I've spent several years of my life trying to make these convergence graphs look as nice as possible. I've only worked with elliptic and saddle-point problems and the guidance probably doesn't apply to parabolic and hyperbolic problems at all.

One thing should be clear already from the comments: there is no reason whatsoever for the graphs to always look like straight lines if you have an upper bound for the error.

Another thing to realize is that such error estimates are usually satisfied by some generic set of problems where the boundary and loading data is specified as an element of a Sobolev space. Thus, there are in reality lots of problems with different loadings, boundary conditions, etc. that satisfy the same estimate—typically with a different constant. When numerically solved on a sequence of meshes, only some of the problems stay near enough the asymptotic limit—as defined by the estimate—for the graphs to look linear.

A third thing to notice is that many times you have multiple asymptotic limits in place, e.g., from different boundary conditions, singularities, and varying regularity of the data. Obviously the estimate with the smallest power of $h$ will dominate in the limit $h \rightarrow 0$ but it is sometimes hard to know that when the other estimates still have an effect. This range where the non-dominating estimates still have an effect is sometimes referred to as the preasymptotic range.

Artificially making the preasymptotic range as small as possible and finding problems that mainly stay near the asymptotic limit is more of an art than a science. Here are some tips:

Make sure that the initial mesh is refined enough to reasonably represent all the nonlinearities in the data. For example, if the loading or boundary condition has $\sin(x)$, make sure that $h << \pi$.
Make sure that there are no singularities other than the ones that you're interested in studying.
Make sure that you don't make any accidental variational crimes, e.g., integrate all functions and data exactly if possible. If not possible, use at least $L^2$ projection to define the boundary and loading data. Be careful with low regularity data, i.e. discontinuous data. Be extra careful also when integrating non-polynomial data.
Make sure that the sequence of meshes is nested. This will often improve the monotonicity of the graph: the next mesh should always be better than the previous one.
Make sure you use high enough integration rule when evaluating the error. I've made this mistake more than once and it took a while to realize where the problem was.
Don't use integration points that reside on element boundaries if you have discontinous data over the element boundary. I've also made this mistake once and it took a while to debug.

In general, I think the error estimates are more useful when checking that your newly implemented finite element code is working correctly and not so useful for applied problems where the data and the boundary conditions are not necessarily nice enough for you to get linear convergence graphs.

If you provide more details on the problem at hand, I think I could provide more specific tips. However, I think these tips will get you started on finding the problem.

Should we always expect FEM error plots to be straight lines?

1 Answers1