7
$\begingroup$

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

| cite | improve this question | |
$\endgroup$
  • 2
    $\begingroup$ 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. $\endgroup$ – Geoff Oxberry Aug 7 '13 at 22:25
  • $\begingroup$ 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). $\endgroup$ – Martin Vymazal Aug 8 '13 at 10:57
  • $\begingroup$ 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?? $\endgroup$ – sebas Mar 6 '14 at 21:59
  • $\begingroup$ 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. $\endgroup$ – Martin Vymazal Mar 7 '14 at 9:16
2
$\begingroup$

It may help to tell us more about your choice of interpolation points. [1] has a 1D discussion on conditioning issues for poor choices of points; the close-to-zero terms aren't usually an issue unless the matrix itself is poorly conditioned or close to singular. Warburton has a paper on optimal choices for interpolation points in 2D and 3D, which I think the Trilinos library bases its higher order Lagrange shape functions off of.

If this isn't the issue, you could always check by switching to a modal representation for a simple test case. Your "slightly off" nodal polynomials can themselves be considered basis functions, and most Galerkin schemes are agnostic to choices of basis function. If you tried both modal and nodal approaches you could see if being inexactly nodal causes instabilities.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ I used equidistant point set. I read about Fekete/warpblend point sets and I agree with you that the choice of interpolation point set is important. What disturbs me is the fact that the entries V[i][j] are not computed accurately enough. I tested with lower-order polynomials so far (up to order 4), so the ill-conditioning of the matrix shouldn't be of extereme importance here. $\endgroup$ – Martin Vymazal Aug 8 '13 at 11:03
  • 2
    $\begingroup$ It seems you might be asking for too much in the accuracy if an error on the order of machine precision is not accurate enough. If that error is not acceptable, then maybe quad precision would be more appropriate for your purposes. You will see errors of this nature though no matter how you decide to evaluate the Vandermonde, even if you did it analytically and hardcoded the values exactly in your code (roundoff). $\endgroup$ – Reid.Atcheson Aug 9 '13 at 19:44
2
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you for this remark. I was thinking about it as well, but I was not sure whether it is safe to clamp off the small values ... $\endgroup$ – Martin Vymazal Aug 8 '13 at 11:04
  • 1
    $\begingroup$ You can do it, but the gain from making this truncation is suspicious. Since floating point errors of this magnitude are inherent in any code which utilizes floating point, there is little justification for truncating entries to force the behavior of the Lagrange basis as expected in exact arithmetic. The only reason I can see really for doing this is if you are assembling a large sparse matrix, and controlling the number of nonzeros is important for memory. $\endgroup$ – Reid.Atcheson Aug 15 '13 at 19:04
  • 1
    $\begingroup$ 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. $\endgroup$ – Wolfgang Bangerth Aug 16 '13 at 2:13

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.