ЧИСЛЕННОЕ РЕШЕНИЕ ДИФФУЗИОННО-ВОЛНОВЫХ
УРАВНЕНИЙ ДРОБНОГО ПОРЯДКА НА КЛАСТЕРНЫХ СИСТЕМАХ
Лукащук С.Ю., Костригин И.В.
Уфимский государственный авиационный технический университет

Введение

В последние 10 лет наблюдается повышенный интерес к исследованию диффузионных процессов, не подчиняющихся закону Фика и не описываемых классическим уравнением диффузии [1,2]. Такие процессы наблюдаются на практике, в частности, в высокоэнергетической плазме, при переносе во фрактальных средах и в аморфных полупроводниках. Оказалось, что для математического описания таких процессов удобно использовать аппарат дробного интегродифференцирования. Однако, наличие в уравнении дробных производных приводит к существенным трудностям при построении как аналитических, так и численных решений. В связи с этим является актуальной задача разработки эффективных алгоритмов численного решения краевых задач для дифференциальных уравнений в частных производных дробного порядка.

Численное решение уравнений дробного порядка сопряжено с большими вычислительными трудностями, связанными, в первую очередь, с нарастанием объемов вычислений по мере увеличения значения временной переменной. Это обусловлено самой природой таких уравнений, которые учитывают так называемые эффекты памяти. В результате время расчета даже одномерной задачи на относительно грубой сетке может составлять несколько часов. Поэтому представляется оправданным привлечение для исследования уравнений с дробными производными многопроцессорной вычислительной техники, что влечет за собой необходимость соответствующей адаптации существующих методов и алгоритмов.

В данной работе рассматриваются некоторые подходы к численному решению одного класса дифференциальных уравнений дробного порядка, называемых диффузионно-волновыми.

1 Алгоритмы численного решения краевых задач
для диффузионно-волнового уравнения

1.1 Постановка краевой задачи

В работе рассматривается численное решение одномерного уравнения аномальной диффузии

$$\frac{\partial^2}{\partial t^2} y(x,t) = D^{1-\gamma}_{t} \frac{\... ...c{\partial}{\partial x}y(x,t) \right), \ \gamma = \alpha - 1, \ 1<\alpha\leq2.$$ (1.1)

в области $\Omega: 0 \leq x \leq L$, где

$$D^{\alpha}f(t)=\frac{d^m}{dt^m} \left[ \frac{1}{\Gamma(m-\alpha)}\int\limits_0^t\frac{f(\tau)}{(t-\tau)^{\alpha+1-m}}d\tau \right],$$

производная Римана - Лиувилля, $m=1,2,...,$ и $0\leq m-1<\alpha\leq m$. Уравнение (1.1) дополняется граничными условиями первого рода (условия Дирихле)

$$\left\{ \begin{aligned} &y(0,t)=g_{0}(t)\\ &y(L,t)=g_{L}(t), \end{aligned} \right.$$

и начальными значениями

$$\left\{ \begin{aligned} &y(x,0)=p_{0}(x)\\ &\frac{\partial}{\partial t}y(x,0)=p_{1}(x). \\ \end{aligned} \right.$$

1.2 Описание и анализ используемых численных схем

Разобьем область $\left[0,L\right]$ на $N+1$ отрезков с узлами $x_{i}=i\Delta x,\ (i=0,...,N)$, где $\Delta x = L/N$. Временной интервал $\left[0,T\right]$ разобьем на части с шагом $\Delta t$ и временными узлами $t_{m}=m\Delta t,\ m=0,1,...,$. Для аппроксимации вторых производных по пространственной и временной переменным используем центральную разность, имеющую второй порядок точности:

$$\frac{\partial^2}{\partial \eta^2}f(\eta_i)=\frac{f(\eta_{i-1})-2f(\eta_i)+f(\eta_{i+1})}{\Delta \eta^2}+O(\Delta \eta^2), \ \ \eta = t, \ x.$$ (1.4)

Для аппроксимации дробной производной Римана-Лиувилля воспользуемся формулой Грюнвалда [9] и L1-аппроксимацией [7]. В результате в первом случае получаем метод конечных разностей (МКР) для дробных производных, а во-втором случае приходим к L1-схеме.

1.2.1 Метод конечных разностей для дробных производных

Известно [3], что для любой функции $f(t)$, которая допускает разложение в степенной ряд, дробная производная порядка $1 - \gamma$ в любой точке сходимости этого ряда может быть записана в виде производной Грюнвалда -Летникова:

$$D^{1-\gamma}_{t} f(t)= \lim\limits_{h\rightarrow0}\frac{1}{h^{1-\... ...a}} \sum\limits_{k=0}^{\left[\frac{t}{h}\right]}\omega_{k}^{(1-\gamma)}f(t-kh).$$ (1.5)

Если $f(t)$ непрерывна, а $f^{'}(t)$ интегрируема на отрезке $[0,t]$, тогда для любого порядка $0 < 1 - \gamma < 1$ производные Римана - Лиувилля и Грюнвалда - Летникова существуют и совпадают для любого момента времени из $[0,t]$ [2].

Определение Грюнвалда - Летникова позволяет численно находить производную Римана - Лиувилля:

$$D^{1-\gamma}_{t} f(t)\approx \frac{1}{h^{1-\gamma}} \sum\limits_{k=0}^{\left[\frac{t}{h}\right]}\omega_{k}^{(1-\gamma)}f(t-kh).$$ (1.6)

Порядок аппроксимации зависит от выбора коэффициентов $\omega_{k}^{(1-\gamma)}$. Аппроксимация будет иметь первый порядок точности, если $\omega_{k}^{(\alpha)}$ - $k$-ый коэффициент в разложении $(1-z)^\alpha$, т.е.

$$(1-z)^\alpha = \sum_{k=0}^{\infty}\omega_{k}^{(\alpha)}z^k.$$ (1.7)

Таким образом,

$$\omega_{k}^{(\alpha)}=(-1)^k\binom{\alpha}{k}=\frac{1}{\Gamma(-\alpha)}\frac{\Gamma(k-\alpha)}{\Gamma(k+1)}.$$

Для $\omega_{k}^{(\alpha)}$ справедливы рекуррентные соотношения:

$$\omega_{0}^{(\alpha)}=1, \ \ \omega_{k}^{(\alpha)}=\left(1-\frac{1+\alpha}{k}\right)\omega_{k-1}^{(\alpha)}.$$

Важно отметить, что аппроксимация (1.6) справедлива, либо при
$t/h \gg 1$, либо когда $f(t)$ достаточна гладкая в окрестности $t=0$ [8].

Для простоты рассмотрим случай $k_{\alpha}\equiv$   const$=1$. Уравнение (1.1), согласно (1.4) и (1.6), заменяется явным конечно-разностным уравнением $(0<j<N)$:

$$y_j^{k+1}-2y_j^k+y_j^{k-1}=\frac{\Delta t^{1+\gamma}}{\Delta x^2} \sum_{l=0}^{k}\omega_l^{1-\gamma}(y_{j+1}^{k-l}-2y_j^{k-l}+y_{j-1}^{k-l}).$$ (1.8)

Следуя [9], отметим, что при указанных коэффициентах $\omega_{k}^{(1-\gamma)}$ схема будет иметь первый порядок точности.

Утверждение 1. Для устойчивости схемы (1.8) достаточно выполнения неравенства

$$\frac{\Delta t^{1+\gamma}}{\Delta x^2} \leq \frac{1}{2^{1-\gamma}}.$$ (1.9)

Д о к а з а т е л ь с т в о. Для анализа устойчивости разностной схемы воспользуемся методом фон Неймана. Предположим, что решение имеет вид $y_j^k=\delta_ke^{ijq\Delta x}$, где $q$ - вещественное пространственное волновое число. Подставим это представление в выражение (1.8) и полученное выражение разделим на $e^{ijq\Delta x}$. C учетом того факта, что

$$e^{iq\Delta x}-2+e^{-iq\Delta x} = \sin^2\left(\frac{q\Delta x}{2}\right),$$

получим

$$\delta_{k+1}-2\delta_{k}+\delta_{k-1} ={} - 4\frac{\Delta t^{1+\g... ...2\left(\frac{q\Delta x}{2}\right) \sum_{l=0}^k\omega_l^{1-\gamma}\delta_{k-l}.$$ (1.10)

Схема безусловно устойчива, если

$$\left\vert \frac{\delta_{k+1}}{\delta_k} \right\vert<1, \ \forall k, q.$$

Положим

$$\delta_{k+1}=\xi \delta_k$$

и предположим, что $\xi=\xi(q)$ не зависит от времени. Тогда из (1.10) получаем уравнение на $\xi$:

$$\xi - 2 + \frac{1}{\xi} = {}- 4\frac{\Delta t^{1+\gamma}}{\Delta ... ...n^2\left(\frac{q\Delta x}{2}\right) \sum_{l=0}^k \omega_l^{1-\gamma} \xi^{-l}.$$ (1.11)

Если $\left\vert\xi\right\vert > 1$ для некоторого $q$, то схема неустойчива. Из (1.11) следует, что если

$$\xi - 2 + \frac{1}{\xi} > 0 \ \ (\Leftrightarrow \xi > 0),$$

то

$$\sum\limits_{l=0}^k\omega_l^{1-\gamma}\xi^{-l} < 0.$$

Заметим, что при $\xi\geq1$ имеем

$$\sum\limits_{l=0}^k\omega_l^{1-\gamma}\xi^{-l} \geq 0$$

и, следовательно, (1.11) выполнено только в области устойчивости. Действительно, в силу (1.7) имеем:

$$\left( 1-\frac{1}{\xi} \right)^\alpha = \sum_{k=0}^{\infty}\omega_{k}^{(\alpha)}\xi^{-k},$$

откуда непосредственно видно, что при $\xi\geq1$

$$\sum\limits_{l=0}^k\omega_l^{1-\gamma}\xi^{-l} \geq 0.$$

Теперь, если $\xi < 0$, то из (1.11) необходимо

$$\sum\limits_{l=0}^k\omega_l^{1-\gamma}\xi^{-l} \geq 0.$$

Используя подход, предложенный в [9], рассмотрим экстремальное значение $\xi=-1$. В результате получим условие устойчивости:

$$\frac{\Delta t^{1+\gamma}}{\Delta x^2} \leq \frac{1}{\sum\limits_{l=0}^k\omega_l^{1-\gamma}(-1)^{-l}}.$$ (1.12)

Обозначим

$$S_{\gamma,k} = \frac{1}{\sum\limits_{l=0}^k\omega_l^{1-\gamma}(-1)^{-l}}.$$

В [9] показано, что

$$\lim\limits_{k\rightarrow\infty}S_{\gamma,k} = S_{\gamma} = \frac{1}{2^{1-\gamma}},$$

на основе чего из (1.12) немедленно получаем (1.9). Утверждение 1 доказано.

1.2.2 L1 схема

Дробная производная аппроксимируется, используя L1 схему, которая справедлива для $0<\gamma \leq 1$. L1 аппроксимация для дробной производной порядка $1 - \gamma$ дается формулой [7]:

$$\begin{split} D^{1-\gamma}_{t} f(t_k) &\approx \frac{\Delta... ...(t_{l-1}))((k-l+1)^\gamma-(k-l)^\gamma) \right\}. \end{split}$$ (1.13)

Оценка точности произведена в статье [10], где показано, что данная аппроксимация имеет первый порядок точности. Опять для простоты рассмотрим случай $k_{\alpha}\equiv$   const$=1$. Запишем явную разностную схему для дробного уравнения аномальной диффузии (1.1) $0<i<N$:

$$\begin{split} y_{j}^{k+1} - 2y_j^k+y_j^{k-1} &= \frac{\Del... ...\nabla y_i^{l-1})\alpha_{k-l+1}(\gamma) \right\}, \end{split}$$ (1.14)

где

$$\begin{split} &\alpha_s(\gamma)=(s)^\gamma-(s-1)^\gamma,\\ & \nabla y_i^k = y_{i+1}^k - 2 y_i^k + y_{i-1}^k. \end{split}$$ (1.15)

Утверждение 2. Для устойчивости схемы (1.14) достаточно выполнения неравенства

$$\frac{\Delta t^{1+\gamma}}{\Delta x^2} \leq \frac{2\Gamma(1+\gamma)}{\gamma2^{\gamma-1}-2^{1+\gamma}+4}.$$ (1.16)

Д о к а з а т е л ь с т в о. Анализ устойчивости разностной схемы проведем аналогично предыдущему пункту. Предположим, как и ранее, что решение имеет вид $y_j^k=\delta_ke^{ijq\Delta x}$. Подставим это представление в выражение (1.14) и разделим результат на $e^{ijq\Delta x}$:

$$\begin{split} \delta_{k+1}-2\delta_{k}+\delta_{k-1} &= {}- 4... ...a_l-\delta_{l-1})\alpha_{k-l+1}(\gamma) \right\}, \end{split}$$ (1.17)

где

$$\rho=\frac{\Delta t^{\gamma+1}}{\Delta x^2 \Gamma(1+\gamma)}.$$

Положим

$$\delta_{k+1}=\xi \delta_k$$

и предположим, что $\xi=\xi(q)$ не зависит от времени. Тогда из (1.17) получаем уравнение на $\xi$:

$$\xi - 2 + \frac{1}{\xi} = {}-4\rho \sin^2\left(\frac{q\Delta x}{2... ...1-\gamma}}+\sum_{l=1}^k(\xi^{l-k}-\xi^{l-k-1})\alpha_{k-l+1}(\gamma) \right\}.$$ (1.18)

Сумма в правой части последнего уравнения может быть упрощена:

$$\begin{split} \sum_{l=1}^k(\xi^{l-k}-\xi^{l-k-1})\alpha_{k-l... ...alpha_l=(\xi-1)\sum_{l=1}^k\frac{\alpha_l}{\xi^l}. \end{split}$$ (1.19)

Выражение (1.18) умножим на $\frac{\xi}{\xi-1}$ и, с учетом (1.19), получим

$$\xi=1-4\rho\sin^2\left(\frac{q\Delta x}{2}\right) \left\{ \frac... ...-\gamma}\xi^{k-1}(\xi-1)}+\sum_{m=0}^{k-1}\frac{\alpha_{m+1}}{\xi^m} \right\}.$$

Рассматривая последнее выражение при экстремальном значении $\xi=-1$, получим оценку

$$\rho \leq \frac{1/2}{ \displaystyle\frac{\gamma}{k^{1-\gamma}(-1)^{k}2}+\sum\limits_{m=0}^{k-1}\displaystyle\frac{\alpha_{m+1}}{(-1)^m} }.$$ (1.20)

Знаменатель в (1.20) убывает для четных $k$ и возрастает для нечетных. Поэтому можно записать:

$$\frac{\Delta t^{1+\gamma}}{\Delta x^2} \leq \frac{\frac{1}{2}\Gamma(1+\gamma)}{\frac{\gamma}{2^{2-\gamma}}+2-2^\gamma}.$$

Последнее выражение после несложных преобразований в точности совпадает с (1.16). Утверждение 2 доказано.

З а м е ч а н и е. Заметим, что как в выражении (1.9), так и в (1.16), при $\gamma \rightarrow 1$ получим

$$\frac{\Delta t^2}{\Delta x^2} \leq 1,$$

что совпадает с условием устойчивости классической разностной схемы для волнового уравнения. В случае $\gamma\rightarrow 0$ получим

$$\frac{\Delta t^2}{\Delta x^2} \leq \frac{1}{2},$$

что также совпадает с условием устойчивости разностной схемы для уравнения теплопроводности.

1.3 Описание параллельного алгоритма

Для метода конечных разностей и для L1 схемы был разработан параллельный алгоритм, реализующий принцип пространственной декомпозиции. Отрезок $\left[0,L\right]$ делится на число отрезков, равное число процессоров, так что области соседних процессоров (номера которых отличаются на единицу) пересекаются по двум точкам, как показано на рисунке 1.

Рис. 1: Разделение отрезка на области.

На каждом временном шаге происходит обмен между двумя соседними процессорами значениями искомой функции в крайних точках отрезках (2 числа типа double) (см. рисунок 2).

Рис. 2: Схема пересылок.

Отметим особенности разработанного параллельного алгоритма:

• количество вычислительных операций пропорционально номеру временному слоя $k$ (см. таблицу 1);
• объем межпроцессорного обмена данными не зависит от номера временного слоя и одинаков в обеих схемах.

Таблица 1: Количество операций, необходимых для вычисления функции в одном узле пространственной сетки.

	Количество операций сложения	Количество операций умножения
МКР	2+3k	4+2k
L1	5+6k	6+3k

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

2 Результаты численных расчетов

Проведем исследование эффективности параллельной реализации схем МКР и L1 на ряде тестовых задач.

2.1 Задача с постоянным коэффициентом диффузии

Рассмотрим случай, когда коэффициент диффузии не зависит от координаты. Не ограничивая общности, можно положить его равным единице $k_{\alpha}=1.$ Начальные и граничные условия поставим следующим образом:

$$\left\{ \begin{aligned} &g_0(t)=g_L(t)=0,\\ &p_0(x)=\sin\left(\frac{\pi x}{L}\right),\\ &p_1(x)=0. \end{aligned} \right.$$

Отметим, что краевая задача (1.1), (2.1) аналогична задаче колебаний струны с закрепленными концами и допускает аналитическое решение вида:

$$y(x,t) = E_\alpha\left(-\frac{\pi^2}{L^2}t^\alpha \right)\sin\left(\frac{\pi x}{L}\right),$$ (2.2)

где

$$E_{\alpha}(z)=\sum_{n=0}^{\infty}\frac{z^n}{\Gamma(\alpha n + 1)}, \alpha>0, z \in \mathbf{C}$$

- функция Миттаг-Леффлера.

2.1.1 Результаты расчетов на кластерной системе УГАТУ

Численные расчеты производились на кластере УГАТУ, имеющем следующие параметры:

533 MHz Alpha 21164EV5/256 Mb/100 Mb Fast Ethernet FD

Рис. 3: Решение (1.1),(2.1) для различных моментов времени

Рис. 4: Временная эволюция центральной точки $x=L/2$ задачи (1.1),(2.1)

Результаты численного решения задачи при $\gamma = 0.5$ для различных моментов времени показаны на рисунке 3. На рисунке 4 показана временная эволюция центральной точки $x=L/2$. Хорошо видно затухание колебаний с течением времени. В этом проявляется существенное отличие диффузионно-волнового уравнения от обычного волнового. В первом случае затухание колебаний оказывается обусловлено структурой среды (напомним, что поскольку коэффициент диффузии постоянный, то среда изотропна). Для классического же волнового уравнения затухание может быть вызвано либо действием внешних сил (источниковых членов), либо неоднородностью среды.

Рис. 5: Изменение абсолютной погрешности численного решения задачи (1.1),(2.1)

На рисунке 5 показано изменение абсолютной погрешности численного решения задачи по схеме L1 для трех моментов времени. Из рисунка видно, что с увеличением времени погрешность убывает. Это объясняется свойствами как самого уравнения ("обладающего памятью"), так и свойствами разностной схемы.

Таблица 2: Время счета на одном процессоре Alpha 533MHz

	$\Delta x =0.01$	$\Delta x=0.005$
Метод конечных разностей	24 с	523 с
L1 схема	119 с	1676 с

В таблице 2 приведены значения времени решения задачи по схемам МКР и L1, при двух различных значениях шагов по пространственной переменной. Отметим, что уменьшение шага в два раза приводит к весьма существенному (более чем в 20 раз) увеличению времени счета. Это связано с резким увеличением количества суммируемых членов ряда, что, в свою очередь, следует из полученных в разделе 2 условий устойчивости.

Рис. 6: Зависимость ускорения от числа процессоров для метода конечных разностей (слева) и L1 схемы (справа)

Рис. 7: Зависимость эффективности от числа процессоров для метода конечных разностей (слева) и L1 схемы (справа)

На рисунках 6 и 7 показаны зависимости ускорения и эффективности параллельных алгоритмов от числа процессоров, полученные при тестовых расчетах на кластерной системе УГАТУ. Видно, что схема L1 показывает более высокие значения ускорения по сравнению с МКР. Более того, в L1 наблюдается почти линейная зависимость ускорения от числа процессоров, слабо зависящая от шага по пространству. Эффективность схемы L1 также оказывается выше. По результатам тестирования можно сделать вывод, что использование процедуры распараллеливания для численных схем решения задач с дробными производными оказывается весьма эффективным.

2.2 Задача с линейным коэффициентом диффузии

Рассмотрим теперь задачу с коэффициентом диффузии, линейно зависящем от пространственной переменной:

$$k_{\alpha}(x)=x.$$

Начальные и граничные условия оставим теми же, что и в предыдущей задаче (см. (2.1)).

На рисунке 8 приведен график численного решения задачи (1.1),(2.1) при $\gamma = 0.5$ в различные моменты времени. Чтобы еще раз подчеркнуть отличие диффузионно-волнового уравнения от аналогичного волнового

$$\frac{\partial^2}{\partial t^2}y(x, t)=\frac{\partial}{\partial x}\left( x \frac{\partial}{\partial x} y(x,t)\right),$$ (2.3)

приведем его решение с такими же начальными и граничными условиями (см. рисунок 9).

Из графиков видно, что решение диффузионно-волнового уравнения затухает со временем, что обусловлено наличием диффузионных эффектов в уравнении (1.1).

Рис. 8: Решение нелинейной диффузионно-волновой задачи (1.1),(2.1), $\gamma = 0.5$

Рис. 9: Решение нелинейной волновой задачи (2.3), (2.1) $\gamma = 0.5$

Список литературы

1

Metzler R., Klafter J. The random walk's guide to anomalous diffusion: a fractional dynamics approach // Phys. Rep. - 2000. - V. 339. - pp. 1-77.

2

Mainardi F., Paradisi P., Gorenflo R. Probability distributions generated by fractional diffusion equations. // In J. Kertesz and I. Kondor (Editors). Econophysics: an Emerging Science (Kluwer, Dordrecht). - 1999. - pp. 312-350.

3

Самко С.Г., Килбас А.А., Маричев О.И. Интегралы и производные дробного порядка и некоторые их приложения. - Минск: Наука и техника, 1987. - 688 c.

4

Золотарев В.М., Учайкин В.В., Саенко В.В. Супердиффузия и устойчивые законы//ЖЭТФ - 1999. - Т.115, в. 4. - C. 1411-1425.

5

Учайкин В.В. Субдиффузия и устойчивые законы//ЖЭТФ - 1999. - Т. 115, в. 6. - С. 2113-2132.

6

Hilfer R. Fractional diffusion based on Riemann-Liouville fractional derivatives//J. Phys. Chem. B. - 2000. - V. 104. - pp. 3914 - 3922.

7

Oldham K., Spanier J. The fractional calculus: theory and applications of differentiation and integration to arbitrary order // Mathematics in Science and Engineering. - Academic Press, New York and London, 1974. - 234 p.

8

Gorenflo R. Fractional calculus: some numerical methods //In: Fractals and Fractional Calculus in Continuum Mechanics, eds. A. Carpinteri and F. Mainardi. - Springer Verlag, Wien, 1997. - N. 378 - pp. 277-290.

9

Yuste S.B., Acedo L. On an explicit finite difference method for fractional diffusion equations // SIAM J. Numer. Anal. - 2005. - V. 42, N. 5. - pp. 1862-1874.

10

Langlands T.A.M., Henry B.I. The accuracy and stability of an implicit solution method for fractional diffusion equation // Journal of Computational Physics. - 2005. - V. 205. - pp. 719-736.

11

Mainardi F., Gorenflo R. On Mittag-Leffler-type functions in fractional evolution processes// Journal of Computational and Applied Mathematics. - 2000. - V. 118 - pp. 283-299.