- •5. Методы решения уравнений
- •5.1. Метод разделения переменных
- •Решить уравнение Пуассона
- •Учитывая, что
- •Результаты решения краевой задачи
- •Собственные значения дискретного оператора Лапласа
- •Преобразованные по Фурье правые части
- •5.2. Методы решения многомерных уравнений
- •Результаты решения краевой задачи
Решить уравнение Пуассона
в
квадрате
при нулевых граничных условиях:
.
Решение задачи производим в следующей последовательности:
Разбиваем стороны квадрата на интервалов величиной
.
Записываем решение задачи и правую часть дифференциального уравнения в виде разложения по собственным функциям дискретного оператора Лапласа:
;
;
.
Подставляя указанные разложения в дифференциальное уравнение, будем иметь
,
где
,
а
и
-
дискретные двумерные операторы Лапласа.
Учитывая, что
,
получим
.
С учётом ортогональности собственных функций уравнения для всех членов разложения записываются в виде
,
,
;
.
Заменяя в уравнениях оператор Лапласа конечно-разностным выражением, получим для каждого систему трёхчленных алгебраических уравнений
,
;
;
.
Собственные значения дискретного оператора Лапласа известны
.
Преобразовывая
каждое из уравнений системы для
,
приведём его к виду
,
Полученную
систему совместно с нулевыми граничными
условиями решаем методом прогонки с
минимальной затратой математических
операций. Подобные решения осуществляем
для всех
членов
разложения. Восстанавливаем искомую
функцию, используя приведённое выше
выражение
Результаты
решения при
,
;
представлены
ниже.
Программа решения краевой задачи
методом разделения переменных
n1=17; n2=17; hx=1./(n1-1); hy=1./(n2-1); f(1:n1,1:n2)=0.; f(8:10,8:10)=100;
for k=2:n2-1
r1=(k-1)*pi/(2.*(n2-1)); r2=sin(r1); lam(k-1)=4.*r2*r2/(hy*hy);
end
fi(1:n1,1:n2)=0.;
for i=8:10
for k=2:n2
fi(i,k-1)=0;
for j=2:n2
r1=(k-1)*pi*(j-1)/(n2-1); r2=sin(r1); r3=f(i,j);
fi(i,k-1)=fi(i,k-1)+r2*r3;
end
end
end
for k=2:n2-1
for i=1:n1
a(i)=1.; b(i)=1.; c(i)=2.+hx*hx*lam(k-1); ff(i)=hx*hx*fi(i,k-1);
end
alf(2)=0.; bet(2)=0.;
for i=2:n1
r=c(i)-alf(i)*a(i); alf(i+1)=b(i)/r; bet(i+1)=(ff(i)+bet(i)*a(i))/r;
end
y(n1,k-1)=0;
for i=n1-1:-1:1
y(i,k-1)=alf(i+1)*y(i+1,k-1)+bet(i+1);
end
end
for i=2:n1-1
for j=2:n2-1
r4=0;
for k=2:n2-1
r1=(k-1)*pi*(j-1)/(n2-1); r2=sin(r1); r3=y(i,k-1);
r4=r4+r2*r3; u(i,j-1)=2.*r4/(n2-1);
end
end
end
disp(u);
surf(u);
