# Evaluation of Vandermonde matrix

I would like to construct a Finite Element basis by using a generalized Vandermonde matrix. The idea is to compute the values of a suitable modal basis ('prime basis') at a set of points in reference space and express the values of Finite Element basis functions as a linear combination of the members of the prime basis as explained in the references [1] or [2].

I have the following problem. The primal basis of choice on triangles/tetrahedra is Dubiner basis and uses Jacobi polynomials. When I evaluate this basis, some entries in the Vandermonde matrix are nearly zero ( ~ 10^-17), but it appears that they should be in fact exactly zero. Later on, when I use the Vandermonde matrix to compute the values of Lagrange shape functions, some of the resulting values are again slightly perturbed. The consequence of this is that the sum of all Lagrange shape functions at certain points is not exactly 1, and the sum of derivatives is not exactly 0.

I am wondering if anyone has ever encountered a problem like this and if there's a remedy. I searched in books and on google, but with no success. As I would like to use higher-order shape functions in a CFD code, I'm afraid that the errors in the shape function evaluation will introduce instabilities.

To give a bit of implementation details: I implemented the computation of Dubiner modes in C++ in two different ways:

a) 'direct' implementation using a recursive formula for Jacobi polynomials as explained in [1] and [2]. This is known to fail at the vertex (-1,1) of the reference triangle due to a singularity in the collapsed coordinate transformation.

b) a 'workaround' implementation described in [3]. In accordance with the observations made in this paper, the errors obtained by b) are slightly smaller than in case a).

[1] Jan S. Hesthaven, Tim Warburton: Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications

[2] George Em Karniadakis, Spencer J. Sherwin: Spectral/hp Element Methods for Computational Fluid Dynamics

[3] Robert C. Kirby: Singularity-free evaluation of collapsed-coordinate orthogonal polynomials

• Aren't Vandermonde matrices typically ill-conditioned? I know that in some cases, the ill-conditioning isn't a problem, but this issue night be worth keeping in mind. – Geoff Oxberry Aug 7 '13 at 22:25
• Yes, they can be ill-conditioned. You can avoid this by selecting a 'good' set of interpolation points (described in [1] or the paper mentioned by JLC in his answer). Anyway, my problem is that the values of the matrix are not computed precisely enough (before I actually try to invert it). – Martin Vymazal Aug 8 '13 at 10:57
• I have read on your comments to other answers that you are using equidistant points to evaluate the polynomials. It is a bad idea for degrees higher that 5 or 6. Have you seen this link with some matlab scripts?? – sebas Mar 6 '14 at 21:59
• Thanks for the link. I think it also depends on your shape functions. I didn't use the equidistant point set with Lagrange shape functions of degree higher than 4 for actual computation and I would advise you to use a better basis than standard Lagrange for higher p. For scaling of condition number of Vandermonde matrix with modal bases, look for example in the thesis of Gonçalo Pena, page 33. – Martin Vymazal Mar 7 '14 at 9:16

This happens all the time due to numerical round-off and can't be avoided. In deal.II, we typically just set small elements (relative to some tolerance, say $10^{-14}$) to zero and rescale the matrix rows or columns if we know that their sum should be one.
• That's one reason. The other is that if we know that analytically certain entries need to sum to one but don't, then the error in these entries is biased rather than random. We can eliminate this bias by rescaling things. It's true that at the level of $10^{-14}$ this really makes no difference (we will always have errors of this magnitude) but it also doesn't cost very much. – Wolfgang Bangerth Aug 16 '13 at 2:13