Задача Дирихле для двумерного уравнения Лапласа [Delphi]

[fronnya]

Нужно решить граничную задачу для уравнения Лапласа, написав в Delphi соответствующую программу.
Постановка задачи:
$\frac{\partial^2u}{\partial x^2}+\frac{\partial^2u}{\partial y^2}=0$;

$u(0,y)=20y; u(x,0)=20x^2$;

$u(1,y)=20\cos(2\pi y); u(x;1)=20\cos(2\pi x)$

Предлагается решить задачу методом сеток. Сводится к решению следующей системы:

$u_{i,j}^{(k)}=\frac{1}{4}(u_{i+1,j}^{(k-1)}+u_{i,j+1}^{(k-1)}+u_{i-1,j}^{(k-1)}+u_{i,j-1}^{(k-1)}), i,j=\overline{1,n-1}$, $k$- номер итерации

$u_{0,j}=20y_j; u_{i,0}=20x_i^2; u_{50,j}=20\cos{2\pi y_j}; u_{i,50}=20\cos{2\pi x_i}; i,j=\overline{0,50}$

Область- квадрат, сетка равномерная. Условие выхода из цикла $\max|u_{ij}^{(k)}-u_{ij}^{(k-1)}|<\varepsilon$ Вот вариант моей программы в Delphi:

Используется синтаксис Pascal

Var h,e,max:real; a:TextFile;
    N,i,j,k:integer;
    u,u0,mas:array [0..50,0..50] of real;
    x,y:array [0..50] of real;
    s:string;
    t:byte;
begin
e:=StrToFloat(Edit2.Text);
h:=1/150;
memo1.Lines.Clear;
///////////////сетка////////////////////
x[0]:=0; y[0]:=0;
for i:=0 to 49 do
  begin
    x[i+1]:=x[i]+i*h;
    y[i+1]:=y[i]+i*h;
  end;
///////////////////////////////////////
//////////граничные условия/////////////
  for i:=0 to 50 do begin u0[0,i]:=20*y[i]; mas[0,i]:=0; end;
  for i:=0 to 50 do begin u0[i,0]:=20*sqr(x[i]); mas[i,0]:=0; end;
  for i:=0 to 50 do begin u0[50,i]:=20*Cos(2*pi*y[i]); mas[50,i]:=0; end;
  for i:=0 to 50 do begin u0[i,50]:=20*Cos(2*pi*x[i]); mas[i,50]:=0; end;
///////////////////////////////////////
/////////////Начальные приближения//////////
for i:=1 to 49 do
  begin
    for j:=1 to 49 do u0[i,j]:=0;
  end;
////////////////////////////////////////////
////////////Итерации///////////////////////
 repeat
 //---------------------Одна итерация---------------//
    for i:=1 to 49 do
      begin
        for j:=1 to 49 do
          begin
            u[i,j]:=1/4*(u0[i+1,j]+u0[i-1,j]+u0[i,j+1]+u0[i,j-1]);
            mas[i,j]:=ABS(u[i,j]-u0[i,j]); //Массив из которого будем выбирать максимум разности u-u0
            u0[i,j]:=u[i,j];
          end;
      end;
  //-----------------------------------------------//
  //--------------Выбор максимального элемента----------//
 max:=mas[0,0];
  for i:=1 to 50 do
    begin
      for j:=1 to 50 do
        begin
          if  mas[i,j]>=max then max:=mas[i,j];
        end;
    end;
  /////////////////////////////////////////////////////
  until max<e;

-- 20.03.2016, 13:26 --

Результат удалось вывести только в виде таблицы $(x_i,y_i,u_{ij})$, чиселки там сомнительные получились. Среди них есть числа порядка $o(10^{100})$ Может, есть ошибка в решении системы?

[Progger]

Первое, что заметил - Вы вычисляете u[i,j], а затем сразу присваиваете его в u0. Таким образом при вычислении следующего элемента используется только что вычисленное значение предыдущего элемента, а не значение, которое было после предыдущего шага. Нужно сделать присвоение u0[i,j]:=u[i,j] в отдельном цикле. Проверьте ещё, что Вы правильно заполнили массивы x и y, возможно в них не то, что было задумано. Для нахождения максимального элемента не обязательно заводить отдельный массив, можно его определять в основном цикле.

[fronnya]

Progger, спасибо, уже все исправил, работает!

-- 22.03.2016, 22:45 --

Кому интересно, вот решение, я набор точек "руками" экспортировал в mathematica: