Численные методы решения дифференциальных уравнений в частных производных
1.4.1. Использование метода сеток для решения параболических уравнений в частных производных
Одним из наиболее распространенных численных методов решения уравнений в частных производных является метод сеток. В методе сеток область Ω, в которой необходимо найти решение уравнения, разобьем прямыми, параллельными осям 
и 
, на прямоугольные области (рис.1.1), где
,
.
Точки, которые лежат на границе области Ω, называются внешними, остальные точки – внутренними. Совокупность всех точек называется сеткой, величины h и Δ – шагами сетки по хи t соответственно.
Рис.1.1. Сетка 
для области Ω
Идея метода сеток состоит в том, что вместо любой непрерывной функции 
будем рассматривать дискретную функцию 
, которая определена в узлах сетки , вместо производных функции будем рассматривать их простейшие разностные аппроксимации в узлах сетки. Таким образом, вместо системы дифференциальных уравнений в частных производных получим систему алгебраических уравнений. Чем меньше величины h иΔ, тем точнее получаемые алгебраические уравнения моделируют исходное дифференциальное уравнение в частных производных. В этом и последующих подразделах этой главы будет рассмотрен метод сеток для каждого из трех типов уравнений и его реализация в MATLAB. Знакомство с численными методами решения дифференциальных уравнений в частных производных начнем с разностных схем решения параболических уравнений.
Разностные схемы решения параболических уравнений будем рассматривать на примере следующего одномерного уравнения:
Построим сетку (см. рис.1.1). Для получения сеточного уравнения заменим производную 
приближенной разностной формулой:
.
В этой и последующих формулах 
– значение функции u в точке 
, 
– в точке 
, 
– в точке 
, 
– в точке 
, 
– в точке 
.
Для замены 
можно воспользоваться одной из приближенных разностных формул
;
.
Кроме того, заменим начальные и граничные условия их разностной аппроксимацией:
Заменив частные производные в задаче соответствующими разностными соотношениями, получим следующую вычислительную схему для расчета значений функции и в узлах сетки :
;
Это явная двухслойная разностная схема (рис.1.2).
Рис.1.2. Шаблон явной двухслойной разностной схемы
Учитывая, что на нулевом слое (при i = 0)все значения 
(как, впрочем, и 
) известны, можно сначала явно рассчитать значения 
, затем 
и так до 
.Для устойчивости разностной схемы значения шагов по t и х должны удовлетворять следующему условию:
.
Пример. Решить параболическое уравнение, описывающее распределение температуры в стержне длиной L,начальная температура стержня задается произвольной функцией 
. Температуры концов стержня равны 
.
Здесь а 2 – коэффициент температуропроводности; 
– коэффициент теплопроводности материала стержня; с – удельная теплоемкость; 
– плотность массы.
Подпрограмма решения задачи с помощью явной разностной схемы в MATLAB:
//Правая часть дифференциального уравнения
//Условие на левой границе
//Условие на правой границе
//Функция решения параболического уравнения методом сеток с помощью явной разностной схемы. N — количество участков разбиения интервала по х (0,1_); К – количество участков разбиения интервала по t (0,Т); а — коэффициент температуропроводности,
// Функция возвращает матрицу решений u вектора х, t
// Вычисляем шаг по х
// Вычисляем шаг по t
// Формируем массив х и первый столбец матрицы решений u
// из начального условия
//Формируем массив t, первую и последнюю строку матрицы решений u из граничных условий
gam = a ^ 2*delta/h ^ 2;
// Формируем матрицу решений u
Входными данными подпрограммы parabol решения задачи являются: N – количество интервалов, на которые разбивается отрезок (О, L); К – количество интервалов, на которые разбивается отрезок (О, Т); L – длина стержня, Т – интервал времени, а – параметр дифференциального уравнения. Функция возвращает три параметра: решение – сеточная функция u, определенная на сетке 
, массивы x и t.
График решения данной задачи изображен на рис.1.3.
Рис.1.3. График решения параболического уравнения при 
1.4.2. Использование метода сеток для решения гиперболических уравнений
Решение гиперболических уравнений также можно осуществить с помощью разностных схем. Разностные схемы решения одномерного гиперболического уравнения рассмотрим на примере следующего уравнения:
Построим сетку , в которой будем искать решение уравнения. Производную 
заменим соотношением:
,
остальные производные такие же, как и в предыдущем подразделе.
Получим следующую явную разностную схему решения уравнения:
Она устойчива при 
и по аналогии с разностной схемой (см. подразд. 1.4.1) может быть легко запрограммирована в MATLAB.
В качестве примера рассмотрим следующую задачу.
Пример. Решить начально-граничную задачу:
Ниже представлена функция giperbol решения данного уравнения, а на рис.1.4 – график полученного решения. Параметры функциианалогичны рассмотренным ранее, т.е. в подпрограмме решения параболического уравнения.
//Функция решения гиперболического уравнения с помощью явной разностной схемы. Входные данные: N — количество участков разбиения интервала (0,L); K — количество участков разбиения интервала (0,T); a — коэффициент теплопроводности; u — матрица решений в узлах сетки; массив x; массив t,
// Вычисляем шаг по х и по t
//Формируем массив x, первый и второй столбцы матрицы решений u из начального условия
//Формируем массив t, первую и последнюю строки матрицы решений u из граничных условий
//Формируем первую и последнюю строки матрицы решений u из граничных условий
//Формируем матрицу решений u
Рис.1.4. График решения гиперболического уравнения
1.4.3. Использование метода сеток для решения эллиптических уравнений
Разностные схемы решения одномерного эллиптического уравнения рассмотрим на следующем примере.
Пример. Рассмотрим разностную схему для эллиптического уравнения в прямоугольной области 
с граничными условиями Дирихле на границе Г:
Проведем в области Ω прямые, параллельные осям 
и 
, где 
; 
Для построения разностного уравнения заменим частные производные и граничные условия следующими соотношениями:
Преобразуем данную эллиптическую краевую задачу к следующей системе разностных уравнений:
Эту систему можно решать итерационными методами (например, методом Зейделя). В случае медленной сходимости итерационных процессов при решении сеточных уравнений, получаемых при аппроксимации гиперболических и эллиптических задач, можно попробовать заменить метод Зейделя градиентными методами (или методами релаксации). Ниже представлено решение эллиптического уравнения сеточным методом, а на рис.1.5 график полученного решения.
Рис.1.5. График решения эллиптического уравнения
// Функция ellip решения задачи
// Входные данные:R, a, b — значения, определяющие область решения задачи; Nx — количество участков разбиения интервала(R-b,R+b); Ny — количество участков разбиения интервала(-a,a); eps — точность решения уравнения методом Зейделя.
// psi — матрица решений в узлах сетки, массив x, массив y, k — количество итераций при решении разностного уравнения
// Вычисляем шаг по y и по x
// Формируем массив x, первый и последний столбцы матрицы решений psi из граничного условия
// Формируем массив y, первую и последнюю строки матрицы решений psi из граничного условия
// Вычисляем коэффициенты разностного уравнения
//Решение разностного уравнения методом Зейделя