Bad accuracy of eigenvalue, how solve it?

Dear freefem community,

Could you help me to understand my mistake or to advise me the best way in freefem++ today to make a sure eigenvalue solving ?

Indeed I would like like to make eigenvalue analysis for thermoelasticity and thermopiezoelectricity. I have made the weak forms and now I am testing with the simplest case, well-known…of elasticity.

You will find here the *.msh from GMSH and the *.edp file.
maillage.7zip (368.2 KB)
eigenvalue test.edp (1.6 KB)

The frequencies from analytical solution for the flexural modes are:
158 kHz / 993 kHz / 2,78 Mhz / 5,46 MHz
with P1 elements it is surestimated by a factor of 1,35!
with P2 elements it yields to: 165,8 kHz (~4,9% error)/ 1,03 MHz (~3,7%)/ 2,9 MHz (~4,3%) / …
with P3 elements it is similar to P2 case results.

Why such a variation between analytical results and the finite element model ?
But the form of eigenvectors are good (saw on Gmsh postprocessing), the main problem is the accuracy of frequencies! Moreover is the P1 element forbidden for these calculations ?

Every advice is welcome. Thanks in advance.
Bruno