Проблема с МКЭ на треугольниках

[Kemi]

Надеюсь, пишу в правильный раздел. Очень прошу помочь.
Задача простейшая, модельная. Решить МКЭ эллиптическое уравнение

$$\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + 2y(1-y) + 2x(1-x)=0$$

на квадратной области [0..1,0..1], разбитой на 4 треугольника (в виде конверта для писем). В узлах, находящихся по углам области заданы первые граничные условия, функция в этих точках равна нулю. То есть интересует только точка одна, в центре области. Линейные базисные функции. Входные данные: три файла.
coords.txt - содержит координаты узлов.
nvtr.txt - каждая строчка соответствует треугольнику, 3 числа - номера узлов глобальной нумерации.
cond_first.txt - граничные значения, глобальный номер узла и соответствуеющее ему значение.
Решаю в системе Maple, чтобы на первых этапах избежать проблем, относящихся исключительно к программированию.

Файл coords.txt

5
0 0
1 0
0 1
0.5 0.5
1 1

Файл nvtr.txt

4
0 1 3
1 4 3
4 2 3
2 0 3

Файл cond_first.txt

4
0 0
1 0
4 0
2 0

Код для системы Maple.

restart:
with(LinearAlgebra):

interface(displayprecision=4):

fd := fopen("coords.txt", READ, TEXT):
coord_count := parse(readline(fd));
coords := readdata(fd, 2);
fclose(fd):

fd := fopen("nvtr.txt", READ, TEXT):
fe_count := parse(readline(fd, integer));
nvtr := readdata(fd, 3, integer);
fclose(fd):

fd := fopen("cond_first.txt", READ, TEXT):
cond_count := parse(readline(fd));
cond_first := readdata(fd, [integer, float]);
fclose(fd):

# переход к индексации массивов с единицы
for i from 1 to fe_count do
  nvtr[i][1] := nvtr[i][1] + 1;
  nvtr[i][2] := nvtr[i][2] + 1;
  nvtr[i][3] := nvtr[i][3] + 1;
end do:

for i from 1 to cond_count do
  cond_first[i][1] := cond_first[i][1] + 1;
end do:

uf := x*(1-x)*y*(1-y);
Ff := -(diff(uf,x,x) + diff(uf,y,y));

u := (x,y) -> x*(1-x)*y*(1-y);
F := (x,y) -> 2*y*(1-y)+2*x*(1-x);

x := Vector(3):
y := Vector(3):

# сборка локальных матриц и векторов правой части
# элементы сразу разносятся по глобальной нумерации
# (то есть локальные матрицы 3x3 записываются как соответствующие им 5x5)
for i from 1 to fe_count do
  for j from 1 to 3 do
    x[j] := coords[nvtr[i][j]][1]:
    y[j] := coords[nvtr[i][j]][2]:
  end do:
  GLocal[i] := Matrix(coord_count):
  FLocal[i] := Vector(coord_count):
  dm := Matrix([[1,1,1],[x[1],x[2],x[3]],[y[1],y[2],y[3]]]):
  mes := abs(Determinant(dm)/2);
  alpha := dm^(-1):
  for k from 1 to 3 do
    for z from 1 to 3 do
      GLocal[i][nvtr[i][k],nvtr[i][z]] := (alpha[k,2]*alpha[z,2] + alpha[k,3]*alpha[z,3])*mes:
    end do:
  end do:
  for k from 1 to 3 do
    FLocal[i][nvtr[i][k]] := F(x[k],y[k]) * 2 * mes / 6:
  end do:
end do:

for i from 1 to fe_count do
  print(GLocal[i], FLocal[i]);
end do:

G := Matrix(5):
F := Vector(5):

# сборка глобальной матрицы и вектора правой части
for i from 1 to fe_count do
  G := G + GLocal[i];
  F := F + FLocal[i];
end do:

print(G, F);

for i from 1 to cond_count do
  for k from 1 to coord_count do
    G[cond_first[i][1],k] := 0;
  end do:
  G[cond_first[i][1],cond_first[i][1]] := 1;
  F[cond_first[i][1]] := cond_first[i][2];
end do:

print(G, F);

LinearSolve(G, F);

Получаемое мной решение в центральной точке (0.5,0.5) 0.0833, но по модельной функции x(1-x)y(1-y) должно быть 0.0625.

Скорее всего, я упустил какой-то момент в сборке локальной матрицы или вектора. Уже ознакомился с двумя книжками. Реализую вроде все так, но где-то затык...

Спасибо большое всем, кто откликнется.

[TOTAL]

Kemi писал(а):

Спасибо большое всем, кто откликнется.

На что откликаться? Вы решали, никаких проблем не возникло.

[Kemi]

Цитата:

На что откликаться? Вы решали, никаких проблем не возникло.

Но ведь решение неверно (еще больше оно отличается от нужного, если пытаться увеличить количество треугольников в сетке). Проблема-то в этом как раз. Не могу понять причину.

[TOTAL]

Разве численное решение должно совпадать с точным?
А вот то, что с увеличением числа элементов погрешность возрастает, свидетельствует об ошибке. Ошибка в составлении системы уравнений, либо в её решении. Для тестирования возьмите ещё линейную функцию (и подправьте уравнение так, что эта функция является решением), для которого линеёные элементы должны дать точный ответ.

[Zai]

Решите задачу для одного треугольника с двумя закрепленными узлами.

[Kemi]

TOTAL, а вы говорите, не на что откликаться. Проблема была как раз в том, что искомая функция не линейная. Точнее сказать, проблемы как таковой и не было. Как-то я упустил это по рассеянности. На линейной $x+y$ считается точно. Даже на большом количестве элементов.

Всем спасибо за помощь и желание.