Доклады Российской академии наук. Математика, информатика, процессы управления, 2021, T. 497, № 1, стр. 31-34

МЕТОД ПОИСКА РЕДУЦИРОВАННОГО БАЗИСА ДЛЯ НЕСТАЦИОНАРНЫХ ЗАДАЧ

И. В. Тимохин^1, 2, *, С. А. Матвеев^1, 2, академик РАН Е. Е. Тыртышников^1, 2, А. П. Смирнов¹

¹ Московский государственный университет имени М.В. Ломоносова
Москва, Россия

² Институт вычислительной математики им. Г.И. Марчука Российской академии наук
Москва, Россия

^* E-mail: m@ivan.timokhin.name

Поступила в редакцию 16.02.2021
После доработки 16.02.2021
Принята к публикации 24.02.2021

DOI: 10.31857/S2686954321020065

Полный текст (PDF)

Аннотация

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

Ключевые слова: уравнение Смолуховского, редукция модели, метод снимков

1. ПОСТАНОВКА ЗАДАЧИ

Мы будем рассматривать метод редукции модели применительно к общего вида системе обыкновенных дифференциальных уравнений

(1)

$\frac{{df}}{{dt}} = F(f(t),t),\quad f(t) \in {{\mathbb{R}}^{N}}.$

Идея метода в такой ситуации состоит в том, чтобы найти подпространство $V \subset {{\mathbb{R}}^{N}}$ размерности существенно меньшей N, на котором хорошо приближается f(t), и построить по (1) систему на координаты такого приближения; при этом численная сложность решения такой системы может оказаться заметно меньше сложности решения полной системы. Построение такого пространства в общем случае требует некоторой априорной информации о решении, что затрудняет использование метода для единичных систем.

В данной раобте мы предложим прямолинейный алгоритм построения такого подпространства на основе метода POD [4] по ходу решения исходной системы, с автоматическим критерием завершения построения. При этом мы будем предполагать дополнительно, что f(t) (приближенно) пробегает искомое пространство за некоторый начальный временной интервал $t \in [0;{{T}_{H}}]$, меньший интересующего нас (например, решение циклично начиная с некоторого момента) $t \in [0;{{T}_{\Pi }}]$, ${{T}_{H}} < {{T}_{\Pi }}$. Работоспособность метода будет продемонстрирована на примере уравнения Смолуховского.

2. ОПИСАНИЕ АЛГОРИТМА

Разобьем весь временной интервал на окна фиксированной длины и будем решать (обычным методом в N-мерном пространстве) систему (1) на каждом окне по очереди. Решив систему на очередном окне, с помощью классического метода снимков и POD, найдем базис, аппроксимирующий решение на этом окне. Сравним полученный базис с уже накопленным (изначально пустым); если накопленный базис приближает новый достаточно хорошо, процесс можно останавливать, искомый базис найден; если достаточно плохо, дополним накопленный базис новым и продолжим расчеты.

В качестве меры качества приближения нового базиса накопленным будем использовать норму проекции

${\text{||}}(I - {{U}_{i}}U_{i}^{*}){{V}_{i}}{\text{|}}{{{\text{|}}}_{2}},$

где V_i – базис текущего окна, а U_i – накопленный (все базисы предполагаются ортонормированными). За оценку качества приближения в алгоритме отвечают два параметра: $\varepsilon {\kern 1pt} ' > \varepsilon > 0$. Если это значение меньше $\varepsilon $, процесс останавливается; если больше $\varepsilon {\kern 1pt} '$ – базис расширяется; в противном случае процесс продолжается без расширения базиса; это важно для предотвращения переобучения на начальном отрезке решения. Подбор подходящего значения $\varepsilon {\kern 1pt} '$, таким образом, связан с необходимостью балансировки скорости построения базиса и опасности переобучения; в наших экспериментах мы использовали $\varepsilon {\kern 1pt} ' = {{10}^{3}}\varepsilon $.

Для дополнения базиса ${{U}_{i}}$ новым базисом V_i прибегнем к сингулярному разложению. А именно, составим из векторов обоих базисов блочную матрицу $A = ({{U}_{i}}\,{\text{|}}\,{{V}_{i}})$ и возьмем в качестве нового базиса левые сингулярные векторы A, соответствующие сингулярным значениям, большим некоторого $\delta > 0$.

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

3. КОНКРЕТНЫЙ ПРИМЕР НЕСТАЦИОНАРНОЙ СИСТЕМЫ

3.1. Формулировка задачи

Мы рассматриваем модель пространственно-однородной системы частиц, эволюционирующей посредством слипания частиц с образованием более крупных. Базовой моделью агрегации частиц в таком случае является формально бесконечная система уравнений, предложенная Смолуховским [5]. В данной работе мы концентрируемся на модели агрегации частиц с источниками и стоком частиц массы большей $N$:

(2)

$\begin{gathered} \frac{{d{{n}_{k}}}}{{dt}} = {{J}_{k}} + \frac{1}{2}\sum\limits_{i + j = k} \,{{C}_{{ij}}}{{n}_{i}}{{n}_{j}} - {{n}_{k}}\sum\limits_{j = 1}^N \,{{C}_{{jk}}}{{n}_{j}}, \\ k = 1,2,...,N. \\ \end{gathered} $

Здесь ${{n}_{k}}$ – концентрация частиц массы $k$, ${{C}_{{ij}}}$ – ядро, характеризующее вероятность столкновения частиц размеров $i$ и j, ${{J}_{k}}$ – мощность источников частиц размера k.

Время вычисления правой части непосредственно по формуле (2) составляет $O({{N}^{2}})$, что задает для решения прикладных задач с (2) непозволительно высокую нижнюю границу сложности. В [3, 6] для отдельных ядер предложены алгоритмы сложности $O(NlogN)$, но для достаточно больших систем (а на практике существует необходимость решать системы с $N$ до 10⁹) желательно иметь возможность сократить вычислительные расходы дальше; нашей целью будет получить сжатое описание решения за время $o(N)$.

3.2. Применение редукции

Покажем теперь, как, построив искомое подпространство, воспользоваться им для ускорения вычислений. Для этого соберем коэффициенты при квадратичных членах в (2) в тензор $S \in {{\mathbb{R}}^{{N \times N \times N}}}$ с элементами

${{S}_{{ijk}}} = \frac{1}{2}\left( {{{\delta }_{{i + j,k}}} - {{\delta }_{{i,k}}} - {{\delta }_{{j,k}}}} \right){{C}_{{ij}}},$

где ${{\delta }_{{i,j}}}$ – символ Кронекера, и запишем уравнение (2) в эквивалентной форме

(3)

$\frac{{d{{n}_{k}}}}{{dt}} = {{J}_{k}} + \sum\limits_{i,j = 1}^N \,{{S}_{{ijk}}}{{n}_{i}}{{n}_{j}}.$

Теперь, если нам известен искомый базис U, мы можем записать

${{n}_{k}} \approx \sum\limits_{i = 1}^R \,{{U}_{{ki}}}{{x}_{i}}$

и преобразовать уравнение (3) в следующий вид:

$\frac{{d{{x}_{\alpha }}}}{{dt}} \approx {{\hat {J}}_{\alpha }} + \sum\limits_{\beta ,\gamma } \,{{\hat {S}}_{{\alpha \beta \gamma }}}{{x}_{\beta }}{{x}_{\gamma }},$

(4)

${{\hat {S}}_{{\alpha \beta \gamma }}} = \sum\limits_{i,j,k} \,{{S}_{{ijk}}}{{U}_{{i\alpha }}}{{U}_{{j\beta }}}{{U}_{{k\gamma }}},$

${{\hat {J}}_{\alpha }} = \sum\limits_i \,{{U}_{{i\alpha }}}{{J}_{i}}.$

Заменив в (4) приближенное равенство на точное, получим систему ОДУ на координаты разложения n_k по базису U.

Сложность вычисления правой части для (4) “в лоб” $O({{R}^{3}})$, что для достаточно малых R может оказаться лучше, чем $O(NlogN)$.

3.3. Численные результаты

Численные эксперименты проводились с системой (2) со следующими значениями параметров:

$\begin{gathered} N = 32\,768, \\ {{J}_{k}} = {{\delta }_{{k1}}}, \\ {{C}_{{ij}}} = {{i}^{a}}{{j}^{{ - a}}} + {{i}^{{ - a}}}{{j}^{a}}, \\ \end{gathered} $

с окнами размером $\tau = 2$ и 65 снимками на каждом окне на интервале [0, 256].

Для этой системы в [1, 2, 7] было найдено циклическое поведение решения, начиная с некоторого момента времени при $a > 1{\text{/}}2$.

В табл. 1 представлено время построения полного решения ${{t}_{\Pi }}$ и редуцированного ${{t}_{{\text{P}}}}$, а также размеры использовавшихся при этом базисов R и момент времени T, до которого строился базис. Как видно из табл. 1, размерность построенного подпространства, а значит и ${{t}_{{\text{P}}}} = O({{R}^{3}})$, существенно зависят от $\varepsilon $, и при высокой точности может ${{t}_{{\text{P}}}}$ даже превышать время построения полного решения. При близком к пограничному для образования циклов значения a для построения базиса с большой точностью потребовалось использовать весь интервал времени решения. Данное наблюдение можно считать весьма естественным, так как по мере приближения значения $a$ к 0.5 быстро растет и период колебаний в системе. Остальные случаи показывают возможность получения приемлемого уровня точности даже с очень небольшим базисом, построенным по достаточно короткому начальному отрезку. Точность решения при этом зависит от $\varepsilon $ в основном опосредовано, через T, а потому не демонстрирует никакой простой зависимости.

Таблица 1.

Сравнение полного решения $n(t)$ и редуцированного решения $\tilde {n}(t)$ системы (2) для N = 32768, с различными значениями параметра ядра a и параметра точности ε (здесь $\varepsilon {\kern 1pt} ' = {{10}^{3}}\varepsilon $). Приведены время полного расчета на всем интервале ${{t}_{\Pi }}$, время расчета редуцированной системы на том же интервале ${{t}_{{\text{P}}}}$, размер полученного базиса R, момент времени T, на котором он был построен, и относительная погрешность редуцированного решения в конце интервала

a	ε	t_Π, c	t_P, c	R	T	$\frac{{{\text{\|\|}}n(T) - \tilde {n}(T){\text{\|}}{{{\text{\|}}}_{2}}}}{{{\text{\|\|}}n(T){\text{\|}}{{{\text{\|}}}_{2}}}}$
0.55	10^–9	4.8 × 10³	75	49	56	2.7 × 10^–3
0.55	10^–11	4.8 × 10³	1.9 × 10³	121	256	9.0 × 10^–8
0.55	10^–13	4.8 × 10³	4.5 × 10³	330	256	5.0 × 10^–10
0.65	10^–9	4.8 × 10³	83	51	62	3.8 × 10^–2
0.65	10^–11	4.8 × 10³	3.5 × 10³	144	224	2.8 × 10^–9
0.65	10^–13	4.8 × 10³	4.5 × 10³	345	250	1.5 × 10^–9

Для решения всех систем дифференциальных уравнений использовался явный метод средней точки с шагом по времени 2^–12; правая часть в полной системе при этом вычислялась по методу, описанному в [3]. Все расчеты проводились на системе с процессором Intel Core i7-7700K, с использованием библиотеки Intel MKL.

ЗАКЛЮЧЕНИЕ

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

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

Zagaynov V.A., Denisenko K., Moskaev A., Lushnikov A.A. Periodical regimes in source-enhanced coagulating systems with sinks // J. Aerosol Sci. 2001. V. 32. P. S983–S984.
Brilliantov N.V., Otienom W., Matveev S.A., Smirnov A.P., Tyrtyshnikov E.E., Krapivsky P.L. Steady oscillations in aggregation-fragmentation processes // Phys. Rev. E. 2018. V. 98. № 1.
Matveev S.A., Smirnov A.P., Tyrtyshnikov E.E. A fast numerical method for the Cauchy problem for the Smoluchowski equation // J. Comput. Phys. 2015. V. 282. № FEB. P. 23–32.
Pinnau R. Model Reduction via Proper Orthogonal Decomposition / In W.H.A. Schilders, H.A. van der Vorst, J. Rommes, ed. Model Order Reduction: Theory, Research Aspects and Applications. Mathematics in Industry. B.-Heidelberg: Springer, 2008, V. 13.
Smoluchowski M.V. Drei vortrage uber diffusion, Brownsche bewegung und koagulation von kolloidteilchen // Zeitschrift fur Physik. 1916. V. 17. P. 557–585.
Timokhin I.V., Matveev S.A., Siddharth N., Tyrtyshnikov E.E., Smirnov A.P., Brilliantov N.V. Newton method for stationary and quasi-stationary problems for Smoluchowski-type equations // J. Comput. Phys. 2019. V. 382. P. 124–137.
Ball R.C., Connaughton C., Jones P.P., Rajesh R., Zaboronski O. Collective oscillations in irreversible coagulation driven by monomer inputs and large-cluster outputs // Phys. Rev. Lett. 2012. V. 109. № 16. P. 168304.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Доклады Российской академии наук. Математика, информатика, процессы управления

Архивы выпусков Информация о журнале Отправить рукопись в журнал