Использование метода сеток для решения гиперболических уравнений
Решение гиперболических уравнений также можно осуществить с помощью разностных схем. Разностные схемы решения одномерного гиперболического уравнения рассмотрим на примере следующего уравнения:
Построим
сетку
,
в которой будем искать решение уравнения.
Производную
заменим соотношением:
,
остальные производные такие же, как и в предыдущем подразделе.
Получим следующую явную разностную схему решения уравнения:
Она
устойчива при
и по аналогии с разностной схемой
(см.подраздел 1.4.1.) может быть легко
запрограммирована в MATLAB.
В качестве примера рассмотрим следующую задачу.
Пример. Решить начально-граничную задачу:
Ниже представлена функция giperbol решения данного уравнения, а на рис.1.4 – график полученного решения. Параметры функции аналогичны рассмотренным ранее, т.е. в подпрограмме решения параболического уравнения.
function [u,x,t]=giperbol(N,K/L,T,a)
//Функция решения гиперболического уравнения с помощью явной разностной схемы. Входные данные: N - количество участков разбиения интервала (0,L); K - количество участков разбиения интервала (0,T); a - коэффициент теплопроводности; u - матрица решений в узлах сетки; массив x; массив t,
// Вычисляем шаг по х и по t
h=L/N; delta=T/K;
//Формируем массив x, первый и второй столбцы матрицы решений u из начального условия
for i=1:N+1
x(i)=(i-1)*h;
u(i,1)=fi(x(i));
u(i,2)=u(i,1)+delta*psi(x(i));
end
//Формируем массив t, первую и последнюю строки матрицы решений u из граничных условий
for j=1:K+1
t(j)=(j-1)*delta;
end
//Формируем первую и последнюю строки матрицы решений u из граничных условий
for j=2:K+1
u(1,j)=0;
u(N+1,j)=fi(L);
end
gam=a^2*delta^2/h^2;
//Формируем матрицу решений u
for j=2:K
for i=2:N
u(i,j+1)=-u(i,j-1)+gam*u(i-1,j)+(2-2*gam)*... u(i,j)+gam*u(i+1,j)+delta^2*f(x(i),t(j));
end
end
end
Рис.1.4. График решения гиперболического уравнения
Использование метода сеток для решения эллиптических уравнений
Разностные схемы решения одномерного эллиптического уравнения рассмотрим на следующем примере.
Пример.
Рассмотрим
разностную схему для эллиптического
уравнения в прямоугольной области
с
граничными условиями Дирихле на границе
Г:
Проведем
в области Ω прямые, параллельные осям
и
,
где
,
Для построения разностного уравнения заменим частные производные и граничные условия следующими соотношениями:
Преобразуем данную эллиптическую краевую задачу к следующей системе разностных уравнений:
Эту систему можно решать итерационными методами (например, методом Зейделя). В случае медленной сходимости итерационных процессов при решении сеточных уравнений, получаемых при аппроксимации гиперболических и эллиптических задач, можно попробовать заменить метод Зейделя градиентными методами (или методами релаксации). Ниже представлено решение эллиптического уравнения сеточным методом, а на рис.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');