The PDE2D High Degree Elements

Everyone knows that high degree elements give better accuracy than low degree elements, for "sufficiently large" problems, when a finite element method is used to solve a PDE. Many people do not realize, however, what a huge effect an increase in degree may have; nor is it commonly known that high degree elements are generally much better even for moderate size problems, and often even for problems with rough solutions.

For 2D problems, PDE2D uses up to 4th degree isoparametric triangular elements, which means you can expect up to O(h5) accuracy in your solution even for problems with curved boundaries or curved interior interfaces. Compare this with the order, O(h2), of the error when centered finite differences or linear finite elements are used.

To illustrate the importance of the use of higher order elements, I solved the eigenvalue problem Uxx+Uyy = λ U in the unit circle, and have recorded below the errors in the first eigenvalue which resulted when different triangulations and element degrees were used. The Harwell sparse matrix routine MA27 (ISOLVE=4), one of the fastest direct solvers in existence, was used to solve the linear systems, and the CPU time reported is on a Cray J-90.

Although it is not fair to compare the different degree elements on the same triangulation, notice that 1000 quartic elements gave an error which is nearly 1 million times smaller than that produced by 10000 linear elements, yet required less than half the CPU time!

2D Problem (true eigenvalue = -5.78318596295)
triangles unknowns degree error CPU time (sec)
100 41 linear 0.24 0.74
100 180 quadratic 0.001 2 0.98
100 418 cubic 0.000 010 1.54
100 755 quartic 0.000 001 7 2.39
1000 469 linear 0.020 4.20
1000 1937 quadratic 0.000 020 6.62
1000 4405 cubic 0.000 000 011 12.32
1000 7873 quartic 0.000 000 002 3 21.17
10000 4873 linear 0.001 7 57.48

For 3D problems, PDE2D uses tri-cubic Hermite basis functions, which results in O(h4) accuracy. Here are some results on a 3D eigenvalue problem (Uxx+Uyy+Uzz=λ U, with U=0 on the boundary of a cylinder of radius 1 and height 1), illustrating the high accuracy you can usually expect for 3D problems. Cylindrical coordinates (P1=r, P2=θ, P3=z) were used.

3D Problem (true eigenvalue = -15.652790364039)
NP1GRID,NP2GRID,NP3GRID unknowns first eigenvalue error CPU time (sec)
5 1000 -15.6548069 0.002 016 5 10
10 4000 -15.6528712 0.000 080 8 30
20 16000 -15.6527944 0.000 004 1 192

This time the parallel band solver was used, and times are on an IBM P690, using 8 processors.

The 1D (collocation) solver also uses cubic Hermite basis functions, so it also produces solutions with O(h4) accuracy. Since the first eigenfunction of the Laplacian on the unit circle is a function of r only, we can recompute the first eigenvalue by solving the 1D problem: Urr + Ur/r = λ U, with U(1)=0 and no boundary condition at r=0. The results are shown in the table below:

1D Problem
NXGRID unknowns error CPU time (sec)
25 50 0.000 000 267 0.12
50 100 0.000 000 015 4 0.17
100 200 0.000 000 000 50 0.27