Использование метода сеток для решения эллиптических уравнений

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

Пример. Рассмотрим разностную схему для эллиптического уравнения в прямоугольной области с граничными условиями Дирихле на границе Г:

Проведем в области Ω прямые, параллельные осям и , где ,

Для построения разностного уравнения заменим частные производные и граничные условия следующими соотношениями:

Преобразуем данную эллиптическую краевую задачу к следующей системе разностных уравнений:

Эту систему можно решать итерационными методами (например, методом Зейделя). В случае медленной сходимости итерационных процессов при решении сеточных уравнений, получаемых при аппроксимации гиперболических и эллиптических задач, можно попробовать заменить метод Зейделя градиентными методами (или методами релаксации). Ниже представлено решение эллиптического уравнения сеточным методом, а на рис.1.5 график полученного решения.

Рис.1.5. График решения эллиптического уравнения

function [psi,x,y,k]=ellip(R,a,b,Nx,Ny,eps)

// Функция ellip решения задачи

// Входные данные:R, a, b - значения, определяющие область решения задачи; Nx - количество участков разбиения интервала(R-b,R+b); Ny - количество участков разбиения интервала(-a,a); eps - точность решения уравнения методом Зейделя.

// Выходные данные:

// psi - матрица решений в узлах сетки, массив x, массив y, k - количество итераций при решении разностного уравнения

// Вычисляем шаг по y

hy=2*a/Ny;

// Вычисляем шаг по x

hx=2*b/Nx;

// Формируем массив x, первый и последний столбцы матрицы решений psi из граничного условия

for i=1:Nx+1

x(i)=R-b+(i-1)*hx;

psi(i,1)=0;

psi(i,Ny+1)=0;

end;

// Формируем массив y, первую и последнюю строки матрицы решений psi из граничного условия

for j=1:Ny+1

y(j)=-a+(j-1)*hy;

psi(1,j)=0;

psi(Nx+1,1)=0;

end;

// Вычисляем коэффициенты разностного уравнения

A=2/hy^2+2/hx^2;

D=1/hy^2;

for i=2:Nx+1

B(i)=1/hx^2+5/(2*hx*x(i));

C(i)=1/hx^2-5/(2*hx*x(i));

end

//Решение разностного уравнения методом Зейделя с точностью eps

p=1;k=0;

while p>eps

for i=2:Nx

for j=2:Ny

V=1/A*(B(i)*psi(i-1,j)+C(i)*psi(i+1,j)+D*(psi(i,j-1)...

+psi(i,j+1))+2);

R(i,j)=abs(V-psi(i,j));

psi(i,j)=V;

end

end

p=R(2,2);

for i=2:Nx

for j=2:Ny

if R(i,j)>p

p=R(i,j);

end

end

end

k=k+1;

end

endfunction

//Вызов функции решения задачи.

[PSI,X,Y,K]=ellip(18,3,6,32,16,0.01);

//Построение графика функции

surf(X,Y,PSI');

xlabel('X');

ylabel('Y');