Application of high-performance computing to solve the Cauchy problem with the fractional Riccati equation using an nonlocal implicit finite-difference scheme

Full Text

Введение

Численное решение уравнений при математическом моделировании динамических процессов может создавать высокую вычислительную нагрузку на узлы ЭВМ. В частности при решении задач возникающих в дробной динамике может возникать нелокальность и даже переменная нелокальность т.е. зависимостью текущего значения решения от конечного числа предыдущих на временном отрезке. Данное явление называют эффектом памяти [?, ?, ?]. Что будет серьезно замедлять вычисления, так как для его корректного обсчета в алгоритме могут возникать строго последовательные участки, которые нельзя распараллелить. Определение и математическое описание эффекта памяти, также эквиваленитное понятиям: эредитарности, запаздывания, наследственности, было дано В. Вольтерра в 1912 г. [?]. Память с точки зрения математики, может описываться с помощью дробной производной. На сегодняшний день существуют самые разные определения понятия оператора дробного дифференцирования: Лиувилля, Римана-Лиувилля, Вейля, Грюнвальда-Летникова и многие другие. Причём для некоторых из них существуют ещё большие обобщения на случай VO (variable order) порядка [?, ?], что в свою очередь приводит к переменной нелокальности. Обобщение известных и разработка новых математических моделей с учетом эредитарности позволяет заметно уточнить известные результаты [?, ?]. На сегодняшний день уже во многих областях науки и техники успешно применяются производные и интегралы нецелых порядков [?, ?].

Так например, в данном исследовании рассматривается прямая задача Коши для нелинейного дробного уравнения Риккати с непостоянными коэффициентами, которая разрешается численно с помощью методов конечно-разностных схем [?]. Память описывается с помощью оператора дробной производной типа Герасимова-Капуто переменного порядка по времени. Такие модели еще называют нелокальными по времени [?]. Одним из способов решения является нелокальная явная конечно-разностная схема (EFDS). И с помощью технологий OpenMP удалось уменьшить в 3-5 раз время вычисления по EFDS [?]. Другой способ заключается в использовании неявных конечно-разностно схем (IFDS) и методов их решения [?].

В данной работе представлен анализ эффективности параллельной программной реализации решения IFDS методом Ньютона, с использованием программно-аппаратной архитектуры OpenMP – открытого стандарта программного интерфейса для параллельных систем с общей памятью [?]. OpenMP является набором директив для компиляторов, библиотек функций и переменных окружения, и реализует параллельные вычисления с помощью идеи многопоточности – нескольких параллельно выполняемых задач CPU потоках. В качестве основного языка программирования использован язык высокого уровня С, из-за его универсальности, широкими возможностями по работе с памятью и официальной поддержке работы с OpenMP [?]. Анализ проводится на основе данных полученных в результате серии вычислительных экспериментов с использованием разработанного программного кода. Вычисления проводились на вычислительном сервере NVIDIA DGX STATION, расположенном в институте математики имени В.И. Романовского АН РУз, а также персональном ноутбуке.

План исследования статьи следующий: Сначала представляется математическое описание задачи. Дается определение дробной производной типа Герасимова-Капуто переменного порядка и приводится задача Коши для дробного уравнения Риккати; Далее описывается нелокальная неявная конечно-разностная схема (IFDS), на примере дробного уравнения Риккати. Приводится метод Ньютона для решения IFDS и критерий сходимости метода; Даются определения понятий необходимых для анализа. Приводятся в графиках и таблицах результаты анализа эффективности параллельной OpenMP реализации IFDS на основе данных серии вычислений, в сравнении с лучшей последовательной реализацией IFDS.

Постановка задачи

Вопросами о возможности обобщения операций интегрирования и взятия производной нецелого порядка задавались ещё Г. Лейбниц, Г. Лопиталь и Я. Бернули в конце 17 века. В частности, их интересовал вопрос о смысле и значении теоремы о дифференцировании произведения двух функций $d^{p}f\left(x\right)/d{x}^p$, при нецелом p значении. Тема активно развивается, с тех пор и по наши дни. К настоящему моменту, существует не один десяток определений дробного дифференцирования, имеющих тот или иной смысл, и находящих свое примирение. Например, Лиувилля, Римана-Лиувилля, Вейля, Грюнвальда-Летникова и многие другие. Причем для некоторых из них существуют обобщения на случай VO (переменного) порядка [?, ?]. Можно отметить труд Самко С., Килбаса А., Маричева О. от 1987 года <<Интегралы и производные дробного порядка и некоторые их приложения>> [?], который подводит итоги исследования данного направления математики за последние 300 лет, и позволяет ознакомится с историей возникновения дробно-дифференциального исчисления.

В этом разделе мы рассмотрим постановку задачи Коши, численные способы решения которой, были ранее подробно изучены автором в работах [?, ?, ?]. Отметим, что дробное исчисление тесно взаимосвязано с понятием наследственности или памяти, рассматриваемой системы, когда текущее ее состояние зависит от конечного числа предыдущих состояний.

Определение [1] Дробной производной типа Герасимова-Капуто [?] где $0<\alpha \left(t\right)<1$ – переменный порядок дробной производной, а функция $u\left(t\right)\in C\left[\mathrm{0,}T\right]$ будем называть дробную производную вида:

${\partial }_{0t}^{\alpha \left(t\right)}u\left(\sigma \right)=\frac{1}{\text{Γ}\left(1-\alpha \left(t\right)\right)}{\int }_0^t\frac{\dot{u}\left(\sigma \right)}{{(t-\sigma )}^{\alpha \left(t\right)}}d\sigma ,$ (1)

где производная $\dot{u}\left(t\right)=\frac{du}{dt}$, а $t\in \left[\mathrm{0,}T\right]$ – текущее время моделирования, $T>0$ – общее время моделирования, а $\text{Γ}\left(.\right)$ – гамма-функция Эйлера.

Данное обобщение (1) на случай $\alpha \left(t\right)$ известной дробной производной Герасимова-Капуто [?, ?] проведено авторами в работе [?]. Согласно [?] можно обобщить (1) на случай $n<\alpha \left(t\right)<n+1$.

Рассмотрим задачу Коши [?] для квадратично нелинейного уравнения Риккати с дробной производной переменного порядка [?]:

${\partial }_{0t}^{\alpha \left(t\right)}u\left(\sigma \right)-b\left(t\right)u\left(t\right)=-a\left(t\right)u{\left(t\right)}^2+c\left(t\right), u\left(0\right)=u_0,$ (2)

где ${\partial }_{0t}^{\alpha \left(t\right)}u\left(\sigma \right)$ – оператор VO вида (1); $T>0$ – время моделирования; $t\in \left[\mathrm{0,}T\right]$ – текущее время моделирования; $u\left(t\right)\in C^2\left[\mathrm{0,}T\right]$ – неизвестная функция решения; $a\left(t\right)$, $b\left(t\right)$, $c\left(t\right)$ $\in C\left[\mathrm{0,}T\right]$ – заданные функции, определяющие значения коэффициентов в каждый момент времени.

Замечание [1] В работе [?] рассмотрено уравнение (2) с более общей правой частью $-a\left(t\right)u{\left(t\right)}^2+c\left(t\right)=f\left(u\left(t\right),t\right)$. Такая задача уже не обязательно представляет уравнение Риккати, и может описывать широкий класс динамических процессов с переменной памятью в насыщенных средах [?]. Это результат использовался нами в работах [?, ?, ?] при математическом моделировании процессов объёмной активности радона (RVA), где каждый член (2) будет иметь конкретный физический смысл.

Схема численного решения

Решение задачи (2) будет будет даваться с помощью теории конечно-разностных схем [?, ?, ?]. Самое сложное в переходе от дробного уравнения в (2) к его дискретному представлению. В работе [?, p.7] нами предложен несколько вариантов такого перехода для VO Герасимова-Капуто (1) при $t\in \left[t_i,t_{i+1}\right]$, в том числе для случая, когда в (1) подразумевается, что интегрирование производится слева.

Определение [2] Нелокальная неявная конечно-разностная схема (IFDS) (от англ. Implicite Finite-Difference Sheme) примет вид:

$\begin{array}{l}{A}_i\sum _{}_{j=0}^{i-1}{w}_j^i\left(u_{i-j}-u_{i-j-1}\right)+a_i{u}_i^2-b_i{u}_i-c_i=\mathrm{0,}\\ A_i=\frac{h^{-{\alpha }_i}}{\text{Γ}\left(2-{\alpha }_i\right)}, w_j^i={(j+1)}^{1-{\alpha }_i}-j^{1-{\alpha }_i},\\ i=\mathrm{1,...,}N, u_0=C,\end{array}$ (3)

где $u_0=C$ – начальное условие задачи Коши (2); C – известная константа; $h=T/N$ – шаг дискретизации равномерной сетки; N – число узлов сетки; $i=\mathrm{2,...,}N-1$ – номер узла вычислительной сетки; $a_i$, $b_i$, $c_i$ – сеточные аналоги функций коэффициентов задачи Коши.

Достоверность вычислений с применением (3) обоснована использованием теории дробного исчисления, методов вычислительной математики, функционального анализа, а также строгостью математических рассуждений. Отметим, что в работе (см. [?, p.7-10]) сформулированы и доказаны ряд ключевые теоремы о сходимости и устойчивости численной схемы IFDS, но в более общем случае при $-a\left(t\right)u{\left(t\right)}^2+c\left(t\right)=f\left(u\left(t\right),t\right)$.

IFDS (3) представляет собой систему нелинейных алгебраических уравнений, которую будем решать модифицированным методом Ньютона (MNM) (от англ. Modification Newton Method). Метод начинается с представления схемы (3) как итерационной функции:

$F\left(u_i\right)=A_i\sum _{}_{j=0}^{i-1}{w}_j^i\left(u_{i-j}-u_{i-j-1}\right)+a_i{u}_i^2-b_i{u}_i-c_i,$ (4)

для которой составляется итерационный процесс:

$U^{m+1}=U^m-\frac{F\left(U^m\right)}{J\left(U^m\right)},$ (5)

где

• m – номер (или шаг) итерационного процесса;

 

• $U^m={\left(u_1^m\mathrm{,...,}{u}_N^m\right)}^T$ – вектор значений функции решения во всех узлах сетки $i=\mathrm{1..}N$ известного на итерационном m (текущем) шаге;

 

• $F\left(U^m\right)=\left(f_1\left(u_1^m\right)\mathrm{,...,}{f}_N\left(u_N^m\right)\right)$ – значения вычисляемые по (4) с учетом значений $U^m$ решения на m шаге;

 

• $U^0$ – начальное приближение на самом первом шаге $m=0$, вектор необходимый для старта (5) процесса. Некое предполагаемое решение (3) на первом шаге;

 

• $J\left(U^m\right)=\left(J_{ij}\right),i=\mathrm{1,...,}N,j=\mathrm{1,...,}N$ – якобиан представляет собой нижнетреугольную матрицу вида:

$J_{ij}=\left\{\begin{array}{ll}& \mathrm{0,}j\ge i+\mathrm{1,}\\ & A_i-b_i+2a_i{u}_i,j=i,\\ & A_i\left({\omega }_{i-j}^i-{\omega }_{i-j-1}^i\right),j\le i-1.\end{array}\right.$ (6)

на каждой m итерации цикла.

Замечание [2] Итерационный процесс (5) не сработает, если будет иметь место сингулярность в нуле (точка, где функция не определёна или нестабильна) для итерационной функции (4).

Что касается сходимости, то стоит отметить что:

$\left|J\left(U^0\right)\right|=\prod _{}_{i=1}^N\left(A_i{\omega }_1^i-b_i+2a_i{u}_i^0\right)\ne 0.$            

т.е. определитель у матрицы Якоби (6) должен быть отличен от нуля, потому что только тогда для якобиана $J\left(U^0\right)$ существует обратная матрица $J{\left(U^0\right)}^{-1}$, т.е. матрица Якоби (6) – невырожденная. И только в таком случае итерационный процесс (5) будет локально сходится порядком $O\left(h\right)$.

В качестве начального приближения $U^0$, в самом простом случае, можно взять вектор составленный из значений $u_0$ – начального условия задачи Коши (3). Однако лучше будет использовать последнее значение $u_{N-1}$ вычисленное, например, другим методом (схемой EFDS ) [?, ?], как на рисунке 1, при условии её сходимости [?]. Это немного ускорит сходимость MNM, а также решит вопрос с сингулярностью, т.к. функция в точке $u_{N-1}$ будет определена, и сама по себе является некоторым возмущением.

Определение [3] Итерационный процесс (5) продолжается до тех пор пока не будет достигнута заданная точность решения, характеризуемая $\epsilon$. Критерием выхода из цикла будет G наименьшая ошибка, полученного решения, такая что $G\le \epsilon$.

Анализ вычислительной эффективности

Отметим, что в рамках анализа вычислительной эффективноси использовались следующие вычислительные системы. Вычислительный сервер NVIDIA DGX STATION с характеристиками: CPU: Intel(R) Xeon(R) CPU E5-2698 v4 @ 2.20GHz(1.20GHz) (20 CPUs, 40 Threads); RAM: 263855 Mb. Вычислительный сервер находится в Институте математики имени В.И. Романовского АН Уз, (г. Ташкент, Узбекистан). Также в экспериментах использовался ноутбук HP Pavilion Gaming Laptop Z270X 15-dk0xxx с характеристиками: CPU: Intel(R) Core(TM) i5-9300H CPU @ 2.40GHz (4 CPUs, 8 Threads); RAM: 8192 Mb.

Анализ вычислительной эффективности будем проводить на основе данных по замерам времени вычисления, полученным при множестве запусков исследуемых на эффективность алгоритмов. Из леммы Брента [?, ?, ?] следует что любой алгоритм реализованный для выполнения в системах типа PRAM можно модифицировать таким образом, что его можно выполнить в системах PRAM с меньшим числом потоков. Здесь PRAM – параллельная машина с произвольным доступом к памяти, объединяющая  потоков (ядер или CPU), общую память и устройство управления. Это одна из самых распространенных моделей вычислительных систем [?]. Из данной леммы следует, что мы можем сравнивать время вычисления (5) на системах с: одним ядром CPU (последовательный алгоритм) и разным количеством ядер (параллельный алгоритм).

Для замеров времени вычисления рассмотрим тестовый пример, на основе модели (2) со следующими параметрами: $\alpha \left(t\right)=0.9-\frac{0.1t}{T}$;  $a\left(t\right)=\frac{{\text{cos}}^2\left(t\right)}{T}$; $b\left(t\right)=1-\frac{0.1t}{T}$; $c\left(t\right)=\frac{{\text{sin}}^2\left(t\right)}{T}$; $N=2000$; $T=20$. Результат вычисления по которому представлен на рисунке 1 ниже:

 

Рис. 1. Результат вычисления итерационной процедуры MNM (5) для IFDS различными алгоритмами

[Figure 1. Calculation result of iterative procedure MNM (5) for IFDS by different algorithms ]

 

Исследуем, насколько быстрее выполняются алгоритмы реализующие итерационный процесс (5), для чего обратимся к таким понятиям как среднее время работы, ускорение, эффективность и стоимость алгоритма.

В силу того, что мы будем получать немного различные $T_c$ при каждом новом численном эксперименте, а количество экспериментов конечно, то $T_c$ можно считать случайной величиной с некоторой функцией распределения. Тогда для вычисления $T_1\left(N\right)$ и $T_c\left(N\right)$ прибегнем к понятию математического ожидания (среднее значение) [?] для дискретной случайной величины:

$\mathbb{E}\left(T_c\left(N\right)\right)=\sum _{}_{i=1}^L\frac{T_c^i\left(N\right)}{L},$ (7)

где i – индекс эксперимента; L – объем выборки.

 

Рис. 2. Среднее E(·) время расчёта в секундах, при разном c – числе потоков CPU

[Figure 2. Average E(·) calculation time in seconds for different c – number of CPU threads ]

 

Далее в таблице представим данные по $\mathbb{E}\left(\cdot \right)$ среднее (7) из выборки в $L=10$ значений при разном количестве задействованных c потоков CPU, где:

• пусть $N=2000$ – число узлов расчётной сетки для (3), характеризует сложность алгоритма;

 

• $T_1\left(N\right)$ – время вычислений по задаче сложностью N последовательным (Sequential) алгоритмом;

 

• $T_c\left(N\right)$ – время вычислений по задаче сложностью N параллельным OpenMP алгоритмом, на машине $c>1$ потоками CPU;

 

Table 1: Затраченные время  (сек.) и память RAM (Mb) для расчёта итерационной процедуры MNM (5) для IFDS (3) алгоритмом на основе OpenMP [Time spent  (sec.) and RAM memory (Mb) for calculating the iterative procedure MNM (5) for IFDS (3) by an algorithm based on OpenMP.]

IFDS	Notebook		Super PC
i	T1(N)	T6(N)	T1(N)	T6(N)	T17(N)	T28(N)	T39(N)
1	96.74	30.06	98.2	21.32	7,86	11.06	8.3
2	88.58	30.28	97.24	21.26	8.0	10.86	8.3
3	94.07	29.09	98.02	21.37	7.86	10.78	7.46
4	97.04	30.67	97.71	21.2	7.9	11.02	8.37
5	91.13	30.25	97.38	20.86	8.19	10.93	8.35
6	95.19	30.74	99.20	20.94	7.96	10.86	8.49
7	93.50	30.74	96.96	21.0	7.99	10.91	8.5
8	93.73	30.91	97.1	21.04	7.87	10.78	8.39
9	92.35	29.89	97.5	21.84	8.05	10.88	8.3
10	93.00	30.71	98.54	20.81	7.93	10.83	8.4
E(T)	93.53	30.34	97.78	21.16	7.961	10.89	8.286
RAM	24.08

 

Теперь с использованием данных о времени потраченного на выполнение кода мы можем исследовать эффективность параллельного алгоритма. Под эффективностью понимается оптимальное соотношение в координатах: ускорение вычислений – объём занимаемой RAM памяти, по сравнению с последовательной версией алгоритма. Ускорение определяется формулой:

$A_c\left(N\right)=\frac{\mathbb{E}\left(T_1\left(N\right)\right)}{\mathbb{E}\left(T_c\left(N\right)\right)},$ (8)

На основе которой можно определить эффективность от распараллеливания:

$V_c\left(N\right)=\frac{\mathbb{E}\left(T_1\left(N\right)\right)}{c\cdot \mathbb{E}\left(T_c\left(N\right)\right)}=\frac{A_c\left(N\right)}{c}.$ (9)

Результаты вычисления по формулам (8) и (9) представлены на рис. 3 и рис. 4 соответственно:

 

Рис. 3. Ускорение Ac(N) при разном c – числе потоков CPU

[Figure 3. Acceleration Ac(N) at different c – number of CPU threads ]

 

Рис. 4. Эффективность Vc(N) при разном c – числе потоков CPU

[Figure 4. Efficiency Vc(N) at different c – number of CPU threads ]

 

Стоимостно–оптимальный показатель $C{O}_c\left(N\right)$ на рис. 5, можно определить как произведение количества используемых потоков и $T_c\left(N\right)$ времени решения параллельном алгоритмом, пропорциональной сложности наилучшего последовательного алгоритма. Данный показатель определяется формулой:

$C{O}_c\left(N\right)=\frac{\mathbb{E}\left(T_c\left(N\right)\right)\cdot c}{\mathbb{E}\left(T_1\left(N\right)\right)}.$                                               

 

Рис. 5. Стоимостно–оптимальный показатель COc(N) при разном c – числе потоков CPU

[Figure 5. Cost-optimal COc(N) for different c – number of CPU threads ]

 

Заключение

Из результатов анализа эффективности параллельного OpenMP алгоритма IFDS видно, что имеет место существенный прирост скорости вычисления, но при увеличении числа потоков для параллельной обработки, их эффективность заметно снижается. При одинаковых параметрах численного эксперимента обе вычислительные системы выдают искомый результат за приблизительно равное время вычисления. Однако наибольшая эффективность достигается при использовании 10-15 потоков. Из анализа данных видно, что OpenMP параллельная программная реализация IFDS показывает ускорение работы от 9 - 12 раз в зависимости от количества задействованных ядер CPU.

Благодарности

Авторы выражают особую благодарность д.ф.-м.н., проф., Хайотову Абдулло Рахмоновичу за оказанную помощь при организации стажировки в Институте математики им. В.И. Романовского Академии наук Республики Узбекистан.

Авторы выражают особую благодарность к.ф.-м.н., Болтаеву Азизу Кузиевичу за оказанную помощь при организации работ на вычислительном сервере NVIDIA DGX STATION в Институте математики им. В.И. Романовского Академии наук Республики Узбекистан.

Аббревиатуры

CPU	Central Processing Unit
OpenMP	Open Multi-Processing
GPU	Graphics Processing Unit
CUDA	Compute Unified Device Architecture
EFDS	Explicit Finite-difference Method
IFDS	Implicit Finite-difference Method
RVA	Radon Volumetric Activity
NM	Newton Method
RAM	Random Access Memory
PC	Personal Computer
PRAM	Parallel Random Access Machine

About the authors

Dmitrii A. Tverdyi

Institute for Cosmophysical Research and Radio Wave Propagation FEB RAS

Email: romanparovik@gmail.com
ORCID iD: 0000-0001-6983-5258

Ph. D. (Phys. & Math.), Researcher, Laboratory of Electromagnetic Radiation

Russian Federation, 684034, Paratunka, Mirnaya st., 7

Roman I. Parovik

Institute for Cosmophysical Research and Radio Wave Propagation FEB RAS

Author for correspondence.
Email: romanparovik@gmail.com
ORCID iD: 0000-0002-1576-1860

D. Sci. (Phys. & Math.), Associate Professor, Leading Researcher at the Laboratory for Modeling Physical Processes

Russian Federation, 684034, Paratunka, Mirnaya st., 7

References

Supplementary files