Участник:Rpu6/Алгоритм построения приближённого решения нерегулярно вырождающегося эллиптического дифференциального уравнения

Алгоритм построения приближённого решения нерегулярно вырождающегося эллиптического дифференциального уравнения
Последовательный алгоритм
Последовательная сложность	[math]\displaystyle{ O(KN\mathrm{max}(K, N)) }[/math]
Объём входных данных	[math]\displaystyle{ KN + 2N }[/math]
Объём выходных данных	[math]\displaystyle{ KN }[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы	[math]\displaystyle{ O(\mathrm{max}(K, N^2)) }[/math]
Ширина ярусно-параллельной формы	[math]\displaystyle{ O(KN) }[/math]

Автор описания: Д. П. Емельянов.

Свойства и структура алгоритма

Общее описание алгоритма

Алгоритм позволяет построить приближения к точному решению эллиптического дифференциального уравнения с нерегулярным вырождением и аналитическими коэффициентами в прямоугольнике. Точные формулы решения уравнения и их строгое обоснование получены с использованием метода спектрального выделения особенностей. Решение строится в виде ряда. Погрешность алгоритма заключается в замене рядов на конечные суммы.

Алгоритм решает следующую задачу:

[math]\displaystyle{ \begin{cases} y^2u''_{yy}+u''_{xx}+c(y)u'_y-a(y)u=f(x,y),(x,y)\in(0,1)\times(0,b), \\ u(0,y)=u(1,y)=u(x,b)=0,|u(x,0)|\lt +\infty. \end{cases} }[/math]

Функции [math]\displaystyle{ a(y), c(y) }[/math] полагаем аналитическими в круге [math]\displaystyle{ |y|\lt R }[/math], где [math]\displaystyle{ R\gt b }[/math]. Кроме того, [math]\displaystyle{ c(y), a(y) \ge 0 }[/math] в [math]\displaystyle{ [0,b] }[/math], [math]\displaystyle{ c(0)=0 }[/math].

Математическое описание алгоритма

1. Приблизим функцию [math]\displaystyle{ f(x,y) }[/math] конечными суммами ряда Фурье по [math]\displaystyle{ x }[/math]:

[math]\displaystyle{ f(x,y)=\sum_{k=1}^{K}f_k(y)\sin \pi kx. }[/math]

[math]\displaystyle{ f_k(y)=2\int\limits_0^1f(x,y)\sin \pi kx\,dx. }[/math]

Для интеграла следует использовать квадратурную формулу.

2. Приблизим функции [math]\displaystyle{ f_k(y),a(y),c(y) }[/math] полиномами:

[math]\displaystyle{ a(y)=\sum_{n=0}^{N}a_ny^n, c(y)=\sum_{n=1}^{N}c_ny^n, f_k(y)=\sum_{n=0}^{N}f_{kn}y^n. }[/math]

Разложение можно проводить по аналогии с разложением в ряд Фурье, но с использованием полиномов Лежандра.

3. Введём в рассмотрение числа:

[math]\displaystyle{ \lambda_k=\pi^2k^2+a_0, k = \overline{1,K}, }[/math]

[math]\displaystyle{ r_k=\frac{1-c_1 + \sqrt{(c_1-1)^2+4\lambda_k}}{2}, k = \overline{1,K}. }[/math]

4. Построим функции [math]\displaystyle{ \eta_k(y),\varphi_k(y) }[/math], они будут использованы при построении приближения к решению:

[math]\displaystyle{ \eta_k(y)=\sum_{n=0}^{N}\eta_{kn}y^n, \stackrel{\circ}{\varphi}_k(y)=\sum_{n=0}^{N}\varphi_{kn}y^n, }[/math]

[math]\displaystyle{ \begin{matrix} \eta_{kn}=\frac{f_{kn}-\sum_{l=2}^{n+1}c_l(n+1-l)\eta_{k,n+1-l}+ \sum_{l=1}^{n}a_l\eta_{k,n-l}}{n(n-1)+nc_1-\pi^2k^2-a_0},n = \overline{1,N},\\ \varphi_{kn}=\frac{ \sum_{l=1}^na_l\varphi_{k,n-l}- \sum_{l=2}^{n+1}(n-l+1+r_k)c_l\varphi_{k,n-l+1} }{n(n-1)+2nr_k+nc_1} ,n= \overline{1,N},\\ \eta_{k,0}=\frac{f_{k,0}}{-\pi^2k^2-a_0}, \varphi_{k,0}=1, \varphi_k(y)=-\frac{\eta_k(b)}{\stackrel{\circ}{\varphi}_k(b)}\stackrel{\circ}{\varphi}_k(y). \end{matrix} }[/math]

5. Построим приближение к решению:

[math]\displaystyle{ u(x,y)=\sum_{k=1}^{K}\left(\left(\frac{y}{b}\right)^{r_k}\varphi_k(y)+\eta_k(y)\right)\sin \pi kx. }[/math]

За счёт выбора [math]\displaystyle{ K }[/math] и [math]\displaystyle{ N }[/math] можно добиться той или иной точности.

Вычислительное ядро алгоритма

Алгоритм решения задачи можно разделить на три крупные части примерно одинаковой вычислительной сложности:
1. Разложение заданной функции [math]\displaystyle{ f(x,y) }[/math] в ряд Фурье, разложение её коэффициентов Фурье в ряды Тейлора, разложение заданных функций [math]\displaystyle{ a(y), c(y) }[/math] в ряды Тейлора. На этом этапе известные сеточные функции скалярно умножаются на тригонометрические функции или полиномы Лежандра. В первом случае коэффициент является скалярным произведением, во втором - линейной комбинацией соответствующих коэффициентов полиномов Лежандра, помноженных на отвечающие данному полиному скалярные произведения.
2. Построение коэффициентов ряда Тейлора для функций [math]\displaystyle{ \eta_k(y), \varphi_k(y) }[/math] с последующим построением самих функций. Используются формулы этапа 4 пункта 1.2.
3. По построенным функциям [math]\displaystyle{ \eta_k(y), \varphi_k(y) }[/math] строится решение задачи [math]\displaystyle{ u(x,y) }[/math] (формула этапа 5 пункта 1.2).

Макроструктура алгоритма

Макроструктура этого алгоритма фактически совпадает с его вычислительным ядром (см. пункт 1.3).

Схема реализации последовательного алгоритма

Алгоритм принимает на вход сеточные функции [math]\displaystyle{ f(x_i, y_j), a(y_j), c(y_j); x_i = \frac{i}{K}, y_j = \frac{j}{N}; i = \overline{0,K-1}, j=\overline{0,N-1} }[/math].

[math]\displaystyle{ f_k(y_j) = \frac{2}{K}\sum_{l=1}^K f(x_l, y_j) sin (\pi kx_l), k=\overline{1,K},\\ f_{kn} = \sum_{m=0}^{N-1} (\mathrm{coeff}_n L_m) \frac{2}{N} \sum_{l=0}^{N-1} f_k(y_l) L_m(y_l), n=\overline{0,N-1},\\ a_{n} = \sum_{m=0}^{N-1} (\mathrm{coeff}_n L_m) \frac{2}{N} \sum_{l=0}^{N-1} a(y_l) L_m(y_l),\\ c_{n} = \sum_{m=0}^{N-1} (\mathrm{coeff}_n L_m) \frac{2}{N} \sum_{l=0}^{N-1} c(y_l) L_m(y_l). }[/math]

Тут [math]\displaystyle{ \mathrm{coeff}_n }[/math] - коэффициент при n-й степени y в полиноме-аргументе (подразумевается, что алгоритм их знает), [math]\displaystyle{ L_m }[/math] - нормированный полином Лежандра степени m.
Формулы для вычислений на других шагах алгоритма изложены в пункте 1.2. Достаточно заметить, что при все функции вычисляются в узлах заданной сетки, суммирование можно вести, запоминая [math]\displaystyle{ y^{n-1} }[/math] с предыдущего шага, вычисляется лишь [math]\displaystyle{ y^{n-1}y }[/math].

Последовательная сложность алгоритма

Будем оценивать количество операций умножения вещественных чисел. Произведём оценку отдельно для каждого шага алгоритма.
1. [math]\displaystyle{ (K + 1)KN }[/math] умножение для вычисления всех [math]\displaystyle{ f_k(y_j) }[/math] (N возникает в силу необходимости вычислять функцию в каждой точке сетки по y), [math]\displaystyle{ ((N + 1)N + 1)K }[/math] операций для вычисления [math]\displaystyle{ f_{kn} }[/math] (считаем, что значения гармоник и полиномов Лежандра известны заранее, учитываем равенство некоторых скалярных произведений при вычислении), [math]\displaystyle{ 2((N + 1)N + 1) }[/math] операций для вычисления [math]\displaystyle{ a_n, c_n }[/math]. 2. Для вычисления n-го коэффициента требуется [math]\displaystyle{ 5 + 3(n-1) }[/math] операций, для вычисления всех коэффициентов - [math]\displaystyle{ \sum_{n=1}^{N-1} (5 + 3(n-1))=5(N-1) + 3\frac{N-2}{2}(N-1)=(N-1)(5 + 1,5 (N-2)) }[/math], суммирование [math]\displaystyle{ \stackrel{\circ}{\varphi}_k(b) }[/math] требует [math]\displaystyle{ 2N }[/math] операций, коррекция [math]\displaystyle{ \varphi_k(y_j) }[/math] требует [math]\displaystyle{ N + 2 }[/math] операции, построение самих функций - [math]\displaystyle{ N(2N)=2N^2 }[/math] операций. Это всё нужно проделать [math]\displaystyle{ K }[/math] раз.
3. Вычисление значения функции в узле сетки стоит [math]\displaystyle{ 3K }[/math] операций умножения. Это нужно проделать в [math]\displaystyle{ KN }[/math] узлах.
Выделим главный член асимптотики по всем шагам: [math]\displaystyle{ K^2N + KN^2 + 3,5KN^2 + 3K^2N=O(KN\mathrm{max}(K,N)). }[/math]
Если положить [math]\displaystyle{ K=N }[/math], то сложность алгоритма имеет порядок [math]\displaystyle{ O(N^3) }[/math].

Информационный граф

Приведём информационные графы для каждой части алгоритма отдельно (иначе граф чрезмерно громоздкий).
1. Построение коэффициентов [math]\displaystyle{ f_{kn} }[/math] по функции [math]\displaystyle{ f(x,y) }[/math] (тут и далее K = N = 3).

2. Построение коэффициентов [math]\displaystyle{ a_n, c_n }[/math] осуществляется по общей схеме.

3. Построение коэффициентов [math]\displaystyle{ \eta_{kn}, \varphi_{kn} }[/math] при фиксированном k.

4. Построение самих функций [math]\displaystyle{ \eta_k(y), \varphi_k(y) }[/math].

5. Построение решения в некоторой точке [math]\displaystyle{ (x,y) }[/math].

Макроструктуру всего алгоритма можно описать следующим графом:

Ресурс параллелизма алгоритма

Будем оценивать параллельную сложность алгоритма в терминах блоков из пункта 1.7. Предполагаем для простоты, что любые данные, помеченные промежуточными на общем графе, доступны с любой параллельной ветви (для обеспечения этого в программе потребуется синхронизация на промежуточных этапах с обменом информацией). Длинные операции суммирования считаем выполняемыми последовательно (без параллелизма). На практике это имеет смысл в силу малой сложности суммы. Распараллеливание только замедлит программу за счёт обмена данными.
Итак, в первом блоке операций каждый коэффициент Фурье в каждой точке может быть вычислен независимо от других. Числа [math]\displaystyle{ f_{kn} }[/math] также вычисляются независимо при различных K и N. В этом случае параллельная сложность - [math]\displaystyle{ O(\mathrm{max}(K,N)) }[/math]. Допустив суммирование сдваиванием, получим [math]\displaystyle{ O(\log_2 \mathrm{max}(K,N)) }[/math]. Параллелизм по K и N - массовый, в случае суммирования сдваиванием имеем дополнительно скошенный параллелизм.
Вычисление [math]\displaystyle{ a_n, c_n }[/math] не зависит от [math]\displaystyle{ f_{kn} }[/math], его можно выполнить параллельно (конечный параллелизм). Сами же вычисления за счёт массового параллелизма по N упрощаются до параллельной сложности [math]\displaystyle{ O(N) }[/math].
Третий блок, как видно из общего графа, допускает массовый параллелизм по K. Однако, внутренние вычисления, как видно из третьего графа, весьма запутанные и не допускают координатного параллелизма (достаточно заметить в формулах, что каждый следующий n+1-й коэффициент требует знания всех предыдущих для вычисления). В этом случае параллельная сложность - [math]\displaystyle{ O(N^2) }[/math]. С использованием скошенного параллелизма можно потенциально сократить параллельную сложность до [math]\displaystyle{ O(N) }[/math], если все коэффициенты будут считаться параллельно. Но в этом случае после вычисления каждого следующего коэффициента, необходимо передать его в параллельные ветви, вычисляющие другие (со старшими номерами). Это на практике потребует существенных временных расходов на передачу данных и синхронизацию.
Вычисление же [math]\displaystyle{ \eta_k(y), \stackrel{\circ}{\varphi}_k(y) }[/math] независимо в каждой точке y и при каждом индексе k. Массовый параллелизм. Параллельная сложность - [math]\displaystyle{ O(N) }[/math]. Коррекция [math]\displaystyle{ \stackrel{\circ}{\varphi}_k(y) }[/math] при этом имеет параллельную сложность [math]\displaystyle{ O(1) }[/math].
И, наконец, вычисление решение можно вести независимо в каждой точке (массовый параллелизм). Параллельная сложность - [math]\displaystyle{ O(K) }[/math], в случае суммирования сдваиванием - [math]\displaystyle{ O(\log_2 K) }[/math].
Итоговая параллельная сложность алгоритма: [math]\displaystyle{ O(\mathrm{max}(K, N^2)) }[/math]. Использование суммирования сдваиванием и указанного подхода в третьем блоке позволяют уменьшить параллельную сложность до [math]\displaystyle{ O(\mathrm{max}(\log_2 K, N)) }[/math].

Входные и выходные данные алгоритма

На вход алгоритм получает заданные на сетке функции [math]\displaystyle{ f(x,y), a(y), c(y) }[/math] ([math]\displaystyle{ (K+2)N }[/math] вещественных чисел). На выходе у алгоритма заданное на сетке решение [math]\displaystyle{ u(x,y) }[/math] ([math]\displaystyle{ KN }[/math] вещественных чисел). Используется сетка, объявленная в пункте 1.5.

Литература

<references \>