Переход из квадратичной в прямоугольную область в задаче Дирихле для уравнения Лапласа

2. Численное решение двумерной задачи Дирихле для уравнения Лапласа
% Листинг функции Iteration.m 
% Уравнение Лапласа: d^2phi/dx^2+d^2phi/dy^2=0 в квадратной области [0,1]x[0,1]
% Краевые условия: phi(x,0)=10, phi(x,1)=10, phi(0,y)=phi(1,y)=5
function z=interationL(N,Omega,Number_of_Iteration,phi)
% Возвращает значение потенциала на каждом шаге итерационного процесса
h=1/N; % шаг сетки
% Вычисление координат узлов сетки
i=1:N+1;
x(i)=(i-1)*h;
j=1:N+1;
y(j)=(j-1)*h;
% Итерационный цикл
for k=1:Number_of_Iteration
% Проход по узлам сетки
for j=2:N 
    for i=2:N
% Релаксация
phi(i,j)=(1-Omega)*phi(i,j)+Omega/4*(phi(i+1,j)+phi(i-1,j)+phi(i,j+1)+phi(i,j-1));
    end;
end;
if k==1
    q=phi; % Сохранение решения на 1-м шаге итерации
else
    q=cat(2,q,phi); % сохранение каждого следующего решения
end;
end;
z=q;
 
% Команды вызова
N=15; % Число узлов сетки
i=1:N+1;
j=1:N+1;
mu(i,1)=10; % потенциал на левой границе
mu(i,N+1)=10; % потенциал на правой границе
mu(1,j)=5; % потенциал на нижней границе
mu(N+1,j)=5; % потенциал на верхней границе
% Задание начального приближения
kx=2:N;
ky=2:N;
mu(kx,ky)=12;
Omega=1; % Параметр релаксации
Niter=200; % Число итераций
z=interationL(N,Omega,Niter,mu); % Решение уравнения Лапласа
% Вычисление векторов и матриц для построения линий уровня
% Формируется массив z [16x3200], где [16x16] -1 итерация, [16x16] -2
% итерация и т. д. [16x16] - 200 итерация
x(i)=(i-1)/N;
y(j)=(j-1)/N;
[x1 y1]=meshgrid(x,y); % формируется 2 массива 16x16 x1 - одинаковые
% столбцы, y1 - одинаковые строки
K=100; % Номер итерации для вывода линий уровня
N1=(N+1)*K+1; % (1601)
N2=(N+1)*(K+1); % (1616)
A=z(1:N+1,N1:N2); % Выделение K-того решения из общей матрицы решений
% (матрица 16x16)
[C,h]=contour(x1,y1,A,17); % (K100_10_10.bmp)
 
% Параметрия
K=199;% Вывод решения на 199 итерации (K199_10_10.bmp)
N1=(N+1)*K+1;
N2=(N+1)*(K+1);
A=z(1:N+1,N1:N2);
[C,h]=contour(x1,y1,A,17);
 
K=1;% Вывод решения на 1 итерации (K1_10_10.bmp)
N1=(N+1)*K+1;
N2=(N+1)*(K+1);
A=z(1:N+1,N1:N2);
[C,h]=contour(x1,y1,A,17);
 
[C,h]=contour(x1,y1,mu,17); % Вывод начального приближения (K0_10_10.bmp)