# Does the convergence of finite element have limit?

I have been trying to solve a nonlinear PDE related to structural mechanics (nonlinear Timoshenko beam to be precise). I am doing both h-refinement and p-refinement to reach the solution. The nonlinear equations were solved using Newton-Raphson method and I am comparing my solution (transverse displacement at a specified load) with the benchmark solution (the benchmark solution was generated using a highly p-refined model and this benchmark solution was interpolated to the solution nodes using spline interpolation). The relative error was computed using both $L_2$-norm and $L_{\infty}$-norm.

$\text{Relative error}= \frac{||U_\text{FEM}-U_\text{benchmark}||}{||U_\text{benchmark}||}$.

The codes were written in MATLAB.

Having said all these, my problem is, the convergence of the solution (with benchmark) stops after a limit and start oscillating around that limit (See the figure attached for reference- LFEM is FEM with linear shape function, QFEM is FEM with quadratic shape function, PFEMs are different p-refined finite elements). I have verified the solution with commercial FEM package and it matches perfectly with it. I have combed the whole code for error with no avail.

What could be the sources of error? I have observed that the solutions are accurate till some significant digits, and that doesn't change even after the refinements. The behavior is as expected till the relative error reaches $10^{-6}$, and the problem starts afterward.

• Are you using single precision? $10^{−6}$ seems to be close to single precision machine tolerance. – Vikram Oct 30 '16 at 12:56
• I think its not. Matlab's default option is double. I haven't changed it explicitly anywhere to single precision. And i have checked the precision of each variable involved in the computation and all of them are in double precision. Even I had that doubt. Anyway, I will check it once more! Thanks :) – Darth_Vader_210382 Oct 30 '16 at 14:42
• This is an expected behaviour: the "condition number" increases with increasing number of "unknowns", so there should always be a limit where the round-off error becomes dominant. – Stefano M Oct 30 '16 at 23:53
• In double precision it seems rather extreme to have a condition number so bad that only 2 digits of correctness remain. – Nick Alger Oct 31 '16 at 4:21
• I guess, the bad condition number is the real culprit. The condition number of the stiffness matrix is $1.8635\times10^{10}$. I think I may have to use an iterative solver with some pre-conditioners. Any suggestions? – Darth_Vader_210382 Oct 31 '16 at 7:55

The overall error of a finite element solution often has many components. For example, these may include (i) the discretization error due to a finite mesh and/or finite polynomial degree, (ii) the error incurred in solving the linear system that defines the nodal values of your finite element solution, (iii) if your problem is nonlinear, the error incurred by terminating a nonlinear iteration (such as a Newton scheme) after only finitely many iterations, (iv) the error incurred by numerical round-off; and possibly others.

In order for your error truly to go to zero, you need to make sure that every single one of these error contributions goes to zero. In other words, as you refine the mesh, you also need to solve the linear system more accurately, run more nonlinear iterations, and possibly switch from single to double precision. If you don't control each one of these error sources, eventually one of them will come to dominate all of the others. For example, if you only do 3 Newton iterations, then it doesn't matter how accurately you solve each of them: the overall error will simply be dominated by the fact that you stop after 3 Newton iterations.

Finally, there is also the possibility that your "reference" solution is not completely accurate. For example, if you read it from a file into which you stored the nodal values of the reference solution with 6 digits of accuracy, you can not expect that the error between your numerical solution and this "reference" solution will be less than $10^{-6}$ simply because the last digit of the reference solution is no longer accurate.

• Thanks for this detailed reply. I also suspected the nonlinear iteration as a source of error. Therefore, I used a linear model. The error remains almost the same. My benchmark solution was constructed using a highly refined pFEM (polynomial order 100) and all computations were done in double precision. It seems like the error emanates either from the numerical round-off(most probably) or from the interpolation of the benchmark solution. The condition number of the stiffness matrix is in the order of $10^{10}$. Could this be the source of error? What should I do to get rid of it? – Darth_Vader_210382 Nov 1 '16 at 5:15
• Right now I'm using the backslash '\' operator in Matlab for solving the equations. Would it make sense to use some iterative solvers like conjugate gradient with pre-conditioners to solve my issue (owing to the bad condition number). If yes, what would be the right choice of pre-conditioner and the iterative solver? The stiffness matrix involved here is NOT diagonally dominant. – Darth_Vader_210382 Nov 1 '16 at 6:57
• I don't know. It depends on the equation you solve, the formulation you use for it, and other factors. You need to find out what people with this kind of equation typically do. – Wolfgang Bangerth Nov 1 '16 at 14:07

I don't really have a definitive answer to your question. However, after some experimentation, I believe your conclusion regarding the condition number of the stiffness matrix is correct.

I also believe you can improve the accuracy with an iterative solver. A very good preconditioner is to use the factors from the direct solution. In Octave/MATLAB, you can use the following code:

lt = chol(K);
u = pcg(K, rhs, 1e-16, 10, lt', lt);


This is basically a form of iterative refinement.

This isn't the final word on this question but it points in a direction for more experiments.