Доклады Российской академии наук. Математика, информатика, процессы управления, 2021, T. 497, № 1, стр. 31-34
МЕТОД ПОИСКА РЕДУЦИРОВАННОГО БАЗИСА ДЛЯ НЕСТАЦИОНАРНЫХ ЗАДАЧ
И. В. Тимохин 1, 2, *, С. А. Матвеев 1, 2, академик РАН Е. Е. Тыртышников 1, 2, А. П. Смирнов 1
1 Московский государственный университет
имени М.В. Ломоносова
Москва, Россия
2 Институт вычислительной математики
им. Г.И. Марчука Российской академии наук
Москва, Россия
* E-mail: m@ivan.timokhin.name
Поступила в редакцию 16.02.2021
После доработки 16.02.2021
Принята к публикации 24.02.2021
Аннотация
Методы редукции модели позволяют заметно сократить вычислительные затраты при решении больших систем дифференциальных уравнений при помощи перехода к расчетам для специального пространства малой размерности. Эти методы требуют априорной информации о базисе такого маломерного пространства, которую возможно получить лишь при численном решении исходной системы высокой размерности. Основное наблюдение данной работы состоит в том, что на широком классе экспериментально рассмотренных задач агрегационной кинетики базис малой размерности существует, следовательно, редукция возможна. В данной работе мы предлагаем новый и эффективный алгоритм построения искомого базиса редуцированной модели без проведения полного расчета. Предложенный алгоритм позволяет существенно выиграть от использования методов редукции модели даже при решении единичной системы без существенной априорной информации о ней.
1. ПОСТАНОВКА ЗАДАЧИ
Мы будем рассматривать метод редукции модели применительно к общего вида системе обыкновенных дифференциальных уравнений
Идея метода в такой ситуации состоит в том, чтобы найти подпространство $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, найдем базис, аппроксимирующий решение на этом окне. Сравним полученный базис с уже накопленным (изначально пустым); если накопленный базис приближает новый достаточно хорошо, процесс можно останавливать, искомый базис найден; если достаточно плохо, дополним накопленный базис новым и продолжим расчеты.
В качестве меры качества приближения нового базиса накопленным будем использовать норму проекции
где Vi – базис текущего окна, а Ui – накопленный (все базисы предполагаются ортонормированными). За оценку качества приближения в алгоритме отвечают два параметра: $\varepsilon {\kern 1pt} ' > \varepsilon > 0$. Если это значение меньше $\varepsilon $, процесс останавливается; если больше $\varepsilon {\kern 1pt} '$ – базис расширяется; в противном случае процесс продолжается без расширения базиса; это важно для предотвращения переобучения на начальном отрезке решения. Подбор подходящего значения $\varepsilon {\kern 1pt} '$, таким образом, связан с необходимостью балансировки скорости построения базиса и опасности переобучения; в наших экспериментах мы использовали $\varepsilon {\kern 1pt} ' = {{10}^{3}}\varepsilon $.Для дополнения базиса ${{U}_{i}}$ новым базисом Vi прибегнем к сингулярному разложению. А именно, составим из векторов обоих базисов блочную матрицу $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$ до 109) желательно иметь возможность сократить вычислительные расходы дальше; нашей целью будет получить сжатое описание решения за время $o(N)$.
3.2. Применение редукции
Покажем теперь, как, построив искомое подпространство, воспользоваться им для ускорения вычислений. Для этого соберем коэффициенты при квадратичных членах в (2) в тензор $S \in {{\mathbb{R}}^{{N \times N \times N}}}$ с элементами
(3)
$\frac{{d{{n}_{k}}}}{{dt}} = {{J}_{k}} + \sum\limits_{i,j = 1}^N \,{{S}_{{ijk}}}{{n}_{i}}{{n}_{j}}.$Теперь, если нам известен искомый базис U, мы можем записать
и преобразовать уравнение (3) в следующий вид:(4)
${{\hat {S}}_{{\alpha \beta \gamma }}} = \sum\limits_{i,j,k} \,{{S}_{{ijk}}}{{U}_{{i\alpha }}}{{U}_{{j\beta }}}{{U}_{{k\gamma }}},$Заменив в (4) приближенное равенство на точное, получим систему ОДУ на координаты разложения nk по базису U.
Сложность вычисления правой части для (4) “в лоб” $O({{R}^{3}})$, что для достаточно малых R может оказаться лучше, чем $O(NlogN)$.
3.3. Численные результаты
Численные эксперименты проводились с системой (2) со следующими значениями параметров:
Для этой системы в [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 = 32 768, с различными значениями параметра ядра a и параметра точности ε (здесь $\varepsilon {\kern 1pt} ' = {{10}^{3}}\varepsilon $). Приведены время полного расчета на всем интервале ${{t}_{\Pi }}$, время расчета редуцированной системы на том же интервале ${{t}_{{\text{P}}}}$, размер полученного базиса R, момент времени T, на котором он был построен, и относительная погрешность редуцированного решения в конце интервала
a | ε | tΠ, c | tP, 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 × 103 | 75 | 49 | 56 | 2.7 × 10–3 |
0.55 | 10–11 | 4.8 × 103 | 1.9 × 103 | 121 | 256 | 9.0 × 10–8 |
0.55 | 10–13 | 4.8 × 103 | 4.5 × 103 | 330 | 256 | 5.0 × 10–10 |
0.65 | 10–9 | 4.8 × 103 | 83 | 51 | 62 | 3.8 × 10–2 |
0.65 | 10–11 | 4.8 × 103 | 3.5 × 103 | 144 | 224 | 2.8 × 10–9 |
0.65 | 10–13 | 4.8 × 103 | 4.5 × 103 | 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.
Дополнительные материалы отсутствуют.
Инструменты
Доклады Российской академии наук. Математика, информатика, процессы управления