Численное решение эллиптического уравнения МКР и МКЭ

При численном решении одного эллиптического уравнения обнаружена следующая особенность. Уравнение решалось методом конечных разностей и методом конечных элементов (пакет FreeFem++), и решения, полученные двумя методами, сравнивались друг с другом. На некоторых сетках проявлялась особенность, что при измельчении сетки разность двух решений увеличивалась. В связи с этим возникли сомнения в правильности вычислений методом конечных разностей.
Поэтому МКР и МКЭ были протестированы на тестовом примере, для которого точное решение известно.
Рассмотрим уравнение

в области

с граничными условиями третьего рода

(условия на остальных гранях куба выглядят аналогично).

Точное решение:

В следующей таблице приводится максимум модуля разности приближённого и точного решений для метода конечных разностей (N - размер сетки).

N	погрешность
30	0.000514
50	0.000185
70	0.0000945
100	0.0000463
150	0.0000206
200	0.0000111

При дискретизации методом конечных разностей использовалась аппроксимация граничных условий второго порядка. Можно заметить, что при уменьшении шага в 2 раза погрешность уменьшается примерно в 4 раза, что свидетельствует о том, что разностная схема действительно сходится со вторым порядком.

В следующей таблице приводится погрешность для метода конечных элементов для сплайнов 1-й и 2-й степени.

N	погрешность P1	погрешность P2
30	0.11701	0.01105
50	0.16234	0.00982
70	0.14249	0.01122
100	0.176	-

Вычислить решение на сетке размером 100 со сплайнами P2 не удалось, FreeFem++ сообщил об ошибке.

Проведённые вычисления показывают, что FreeFem++ оказался неэффективным по сравнению с методом конечных разностей (по крайней мере, для трёхмерной задачи).

Поэтому интересно было бы сравнить FreeFem++ с другим конечноэлементным пакетом, чтобы выяснить, характерна ли такая погрешность вообще для метода конечных элементов или это особенность FreeFem++.

С исходным кодом в пакете FreeFem++ можно ознакомиться здесь.
Отметим, что для запуска FreeFem++ при большом N может понадобиться 64-битная версия программы.

Дополнение от 31.10.15

Проведем аналогичный численный эксперимент для двумерной задачи.
В таблице указана погрешность для FreeFem++ (сплайны P1 и P2) и для метода конечных разностей. Можно видеть, что погрешность метода конечных разностей уменьшается примерно в 4 раза при уменьшении шага сетки в 2 раза. Для МКЭ со сплайнами P1 коэффициент несколько меньше 4, а для сплайнов P2 уменьшение в 8 раз.

N	МКЭ P1	МКЭ P2	МКР
50	0.00256	5.0903e-6	0.000216
100	0.000776	6.4124e-7	5.4239e-5
200	0.000229	8.0454e-8	1.3561e-5
400	6.6212e-5	1.0066e-8	3.3903e-6

Исходный код можно посмотреть здесь.