Параллельные методы решения дифференциальных уравнений в частных производных
Дифференциальные уравнения в частных производных представляют собой широко применяемый математический аппарат при разработке моделей в самых разных областях науки и техники. К сожалению, явное решение этих уравнений в аналитическом виде оказывается возможным только в частных простых случаях, и, как результат, возможность анализа математических моделей, построенных на основе дифференциальных уравнений, обеспечивается при помощи приближенных численных методов решения.
Объем выполняемых при этом вычислений обычно является значительным, и использование высокопроизводительных вычислительных систем традиционно для данной области вычислительной математики.
Проблематика численного решения дифференциальных уравнений в частных производных является областью интенсивных исследований.
Рассмотрим в качестве учебного примера проблему численного решения задачи Дирихле для уравнения Пуассона, которая определяется как задача нахождения функции u=u(x,y), удовлетворяющей в области определения D уравнению
![\left\{
\begin{aligned}
& \frac{\partial^2u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}= f(x,y), & (x,y) \in D, \\
& u(x,y) = g(x,y), & (x,y) \in D^0,
\end{aligned}
\right.](/sites/default/files/tex_cache/b566fd6e4ba5df99880b05347dc19492.png)
Для простоты изложения материала в качестве области задания D функции u(x,y) далее будет использоваться единичный квадрат
![D=\{ (x,y)\in D:0\le x,y\le 1\}](/sites/default/files/tex_cache/1b7490c1de8e5440571c9e7baf4b0307.png)
11.1. Последовательные методы решения задачи Дирихле
Одним из наиболее распространенных подходов к численному решению дифференциальных уравнений является метод конечных разностей ( метод сеток ) (см., например, [ [ 6 ] , [ 13 ] , [ 60 ] ]). Следуя этому подходу, область решения D можно представить в виде дискретного (как правило, равномерного) набора ( сетки ) точек ( узлов ). Так, например, прямоугольная сетка в области D может быть задана в виде (рис. 11.1)
![\left\{
\begin{gathered}
D_h=\{(x_i,y_j):x_i=ih, y_i=jh, 0\le i, j\le N+1, \\
h=1/(N+1)
\end{gathered}
\right.](/sites/default/files/tex_cache/863918e1354f357775157785253da462.png)
Обозначим оцениваемую при подобном дискретном представлении аппроксимацию функции u(x,y) в точках (xi, yj) через uij. Тогда, используя пятиточечный шаблон (см. рис. 11.1) для вычисления значений производных, мы можем представить уравнение Пуассона в конечно-разностной форме
![\frac{u_{i-1,j}+u_{i+1,j}+u_{i,j-1}+u_{i,j+1}-4u_{ij}}{h^2} = f_{ij}](/sites/default/files/tex_cache/041b0d22ab73cddefcdc7e87b51c6517.png)
Данное уравнение может быть разрешено относительно uij:
uij=0,25(ui-1,j+ui+1,j+ui,j-1-h2fij).
Разностное уравнение, записанное в подобной форме, позволяет определять значение uij по известным значениям функции u(x,y) в соседних узлах используемого шаблона. Данный результат служит основой для построения различных итерационных схем решения задачи Дирихле, в которых в начале вычислений формируется некоторое приближение для значений uij, а затем эти значения последовательно уточняются в соответствии с приведенным соотношением. Так, например, метод Гаусса – Зейделя для проведения итераций уточнения использует правило
![u_{ij}^k=0,25(u_{i-1,j}^k + u_{i+1,j}^{k-1}+u_{i,j-1}^{k}+u_{i,j+1}^{k-1}-h^2 f_{ij})](/sites/default/files/tex_cache/e49d3addb34583cb9dd331b19cff3a03.png)
![Прямоугольная сетка в области D (темные точки представляют внутренние узлы сетки, нумерация узлов в строках слева направо, а в столбцах — сверху вниз)](/EDI/22_05_17_2/1495405295-1783/tutorial/295/objects/11/files/11_01.gif)
Рис. 11.1. Прямоугольная сетка в области D (темные точки представляют внутренние узлы сетки, нумерация узлов в строках слева направо, а в столбцах — сверху вниз)
Рассмотренный алгоритм (метод Гаусса – Зейделя) на псевдокоде, приближенном к алгоритмическому языку С++, может быть представлен в виде:
Алгоритм 11.1. Последовательный алгоритм Гаусса – Зейделя
// Алгоритм 11.1 do { dmax = 0; // максимальное изменение значений u for ( i=1; i<N+1; i++ ) for ( j=1; j<N+1; j++ ) { temp = u[i][j]; u[i][j] = 0.25*(u[i-1][j]+u[i+1][j]+ u[i][j-1]+u[i][j+1]–h*h*f[i][j]); dm = fabs(temp-u[i][j]); if ( dmax < dm ) dmax = dm; } } while ( dmax > eps );
(напомним, что значения uij при индексах i,j=0,N+1 являются граничными, задаются при постановке задачи и не изменяются в ходе вычислений).
Для примера на рис. 11.2 приведен вид функции u(x,y), полученной для задачи Дирихле при следующих граничных условиях:
![\left\{
\begin{aligned}
&\quad f(x,y)=0, \quad (x,y) \in D, \\
&\quad 100-200x, \quad y=0, \\
&\quad 100-200x, \quad x=0, \\
&-100+200x, \quad y=1, \\
&-100+200x, \quad x=1. \\
\end{aligned}
\right.](/sites/default/files/tex_cache/046921e8d25700ae6a1e893f1b45d866.png)
Общее количество итераций метода Гаусса – Зейделя составило 210 при точности решения eps=0,1 и N=100 (в качестве начального приближения величин uij использовались значения, сгенерированные датчиком случайных чисел из диапазона [-100, 100]).