Матрицы специального вида и быстрые алгоритмы при численном решении уравнений типа Смолуховского тема диссертации и автореферата по ВАК РФ 05.13.18, кандидат наук Тимохин Иван Владиславович

• Тимохин Иван Владиславович
• кандидат наук
• 2022

Оглавление диссертации кандидат наук Тимохин Иван Владиславович

1.1 Тема исследования

1.2 Существующие подходы

1.3 Цели и задачи

1.4 Научная новизна

1.5 Теоретическая и практическая значимость

1.6 Методология и методы исследования

1.7 Положения, выносимые на защиту

1.8 Апробация результатов

1.9 Личный вклад автора

1.10 Основное содержание работы по главам

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

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

2.1 Оператор правой части

2.2 Быстрый алгоритм вычисления

3 Метод Ньютона

3.1 Решение для однородной системы

3.2 Метод Ньютона

3.2.1 Предобусловливатель

3.2.2 Выбор шага

3.3 Численные эксперименты

4 Тензоризация

4.1 Общие сведения о TTT-формате

4.1.1 Тензоризация

4.1.2 Операции с тензорами в TTT-формате

4.2 Тензор оператора Смолуховского

4.3 Метод Ньютона

4.4 Численные эксперименты

5 Редукция модели

5.1 Редукция модели

5.2 Построение базиса

5.3 Численные эксперименты

5.4 Усечение базиса

5.4.1 Оценка разности решений

5.4.2 Сравнение с решением бесконечной системы

5.4.3 Численные примеры

6 Редукция в уравнении с переносом

6.1 Постановка

6.2 Численный метод

6.2.1 Перенос

6.2.2 Общая схема

6.2.3 Редукция

6.2.4 Perfect Matching Layer

6.3 Численные эксперименты

7 Заключение

Глава

Введение диссертации (часть автореферата) на тему «Матрицы специального вида и быстрые алгоритмы при численном решении уравнений типа Смолуховского»

Введение

1.1. Тема исследования

Основной проблемой, решению частных случаев которой посвящены последующие главы, является вопрос численного моделирования процесса агрегации частиц. Системы соединяющихся каким-либо образом частиц встречаются в природе довольно часто; примерами могут служить образование небесных тел из протопланетных дисков [1, 2] и сворачивание крови, образование частиц смога в атмосфере [3] и полимеров в растворах [4], створаживание молока и образование мочевины [5]. Кроме того, используемый для моделирования таких процессов математический формализм оказывается полезен и в других, казалось бы, имеющих достаточно отдалённое отношение к слипающимся частицам задачах. К ним относятся, например, рост трещин в структуре материалов, или формирование социальных связей [6].

Существуют различные подходы к моделированию таких систем; так, в простых случаях, когда крупные кластеры обмениваются лишь частицами единичного размера, но не взаимодействуют между собой напрямую, для описания таких систем используют уравнения Беккера-Дёринга, поведение которых достаточно хорошо изучено [7].

Мы же будем рассматривать более общую модель, основанную на пред-

ложенной в 1916 году Марианом Смолуховским [8, 9] формально бесконечной системе обыкновенных дифференциальных уравнений, описывающей эволюцию концентраций частиц различных масс в пространственно-однородной системе, при абсолютно неупругих столкновениях отдельных частиц в результате их хаотического движения. За прошедший век были предложены различные модификации базовой модели — например, в работе [10] была предложен вариант с непрерывным спектром масс — и различные расширенные версии — например, совмещение с распадом частиц [2]. Отдельно отметим обобщение модели Смолуховского на пространственно-неоднородные системы и совмещение её с уравнением переноса [11, 12] — вполне естественное в контексте некоторых из обозначенных выше применений, но весьма затратное с точки зрения построения численного решения.

Первая же сложность, возникающая при попытке численного решения системы уравнений Смолуховского это формально бесконечное число уравнений в ней — классическая модель не подразумевает никакой верхней границы на массу частиц в системе и, в действительности, естественным образом порождает частицы сколь угодно большой массы. (Даже не вдаваясь в подробности устройства модели, это легко увидеть «на пальцах»: поскольку частицы всё время слипаются между собой, максимальная масса частиц в системе постоянно растёт.) Численное решение такой системы «в лоб», конечно, невозможно, а потому на практике приходится ограничивать максимальный размер частиц, «обрезая» систему до конечной.

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

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

Вне зависимости от происхождения конечной системы, таким образом, возникает естественное желание численно решать рассматриваемую систему уравнений как можно большего размера, с как можно большей допустимой массой частиц. Это оказывается достаточно сложным — в исходной формулировке, простое вычисление правой части системы ОДУ типа Смо-луховского с N уравнениями требует 2) операций, что очевидным образом накладывает достаточно жёсткие ограничения на размер решаемых систем.

1.2. Существующие подходы

Для борьбы с настолько стремительным ростом сложности в контексте необходимости решать системы с десятками тысяч, а, в отдельных задачах, и триллионами уравнений, предлагались весьма разнообразные подходы. Так, например, в работах [13, 14, 15] предлагается «огрублять» исходную модель, объединяя большие массы частиц в «окна» растущего с массой размера, и приближая концентрации частиц внутри каждого такого окна какой-то достаточно малопараметрической зависимостью от массы, полученной исходя из некоторых априорных соображений о структуре решения. Такой подход позволяет существенно сократить число параметров, используемых для описания системы — с N до 0(log N) — однако требует доста-

точно точного представления о структуре искомого решения; в отсутствие такой априорной информации точность подхода заметно падает [16]. В эту же категорию «полуаналитических» методов отнесём методы, основанные на производящей функции моментов решения [17].

В некотором роде противоположны рассмотренным методы типа Монте-Карло [18, 19, 20], которые можно рассматривать как варианты непосредственного моделирования системы сталкивающихся частиц; такие методы успешно применялись и к пространственно-неоднородным задачам [21].

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

Опорным результатом, на котором будет строиться большая часть предлагаемых методов, является полученный в работе [22] быстрый алгоритм вычисления правой части в системе уравнений Смолуховского, достигающий на широком классе практически значимых параметров сложности O(N log N) — существенное улучшение по сравнению с O(N2) «наивного» алгоритма, достигнутое без необходимости строить дополнительные предположения об устройстве самого решения. Более подробно этот метод будет рассмотрен в разделе 2.2, вместе с идейно незначительным, но важным для наших целей обобщением на билинейный вариант изначально квадратичного оператора.

1.3. Цели и задачи

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

1.4. Научная новизна

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

1.5. Теоретическая и практическая значимость

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

С теоретической точки зрения работа содержит ряд наблюдений, связанных со структурой уравнений Смолуховского, как например исследование структуры матрицы Якоби оператора правой части и оценки на ранги её составляющих (раздел 3.2.1) или оценки на TT-ранги оператора (раздел 4.2), а также эмпирические по своей природе но всё же представляющие общий теоретический интерес наблюдения, такие как существование маломерного линейного пространства, содержащего в себе решение на всём временном интервале (глава 5) или на части пространственной области (глава 6), а также возможность переноса первого между задачами разной размерности (раздел 5.4), вкупе формирующих потенциальную основу для дальнейших теоретических исследований рассматриваемой модели.

1.6. Методология и методы исследования

Приводимые численные методы строились главным образом на основе известных в вычислительной линейной алгебре подходов и разложений. Для эмпирической проверки их работоспособности все методы были реализованы автором в виде программ на языке ^—Ь и протестированы на пробных задачах.

1.7. Положения, выносимые на защиту

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

— обобщение алгоритма быстрого вычисления квадратичного оператора [22] на билинейный, и связанный с этим алгоритм быстрого вычисления произведения матрицы Якоби оператора агрегации на вектор,

— построенный для матрицы Якоби предобусловливатель, и

— программную реализацию построенного метода, с помощью которой эффективность метода была проверена на численных экспериментах;

• Реализация метода Ньютона для той же задачи на основе тензорных разложений, включая

— явный вид TT-разложения «скелета» оператора агрегации, и

— программную реализацию метода, с помощью которой была проверена работоспособность метода в отдельных случаях;

• Проверка применимости редукции модели к решению системы ОДУ типа Смолуховского, включая

— алгоритм построения базиса редукции с автоматическим обнаружением циклов в решении,

— программную реализацию указанного алгоритма в сочетании с редуцированной формой системы, и

— эмпирическое наблюдение о связи базисов редукции в системах разного размера;

• Численный метод моделирования агрегации с переносом на основе редукции по пространству, включая

— метод пересчёта решения на одном шаге по времени с контролируемой погрешностью, и

— программная реализация.

1.8. Апробация результатов

Основные результаты работы представлены в следующих публикациях:

1. Newton method for stationary and quasi-stationary problems for Smolu-chowski-type equations / I. V. Timokhin, S. A. Matveev, N. Siddharth et al. // Journal of Computational Physics. — 2019. — Vol. 382. — P. 124-137.

2. I. Timokhin. Tensorisation in the solution of smoluchowski type equations // Large-Scale Scientific Computing. — Vol. 11958 of Theoretical Computer Science and General Issues. — Springer International Publishing, 2020. — P. 181-188.

3. МЕТОД ПОИСКА РЕДУЦИРОВАННОГО БАЗИСА ДЛЯ НЕСТАЦИОНАРНЫХ ЗАДАЧ / И. В. Тимохин, С. А. Матвеев, Е. Е. Тыр-тышников, А. П. Смирнов // Доклады Российской академии наук. Математика, информатика, процессы управления. — 2021. — Т. 497, № 2_21. — С. 35-38.

4. Model reduction for Smoluchowski equations with particle transfer / I. Timokhin, S. Matveev, E. Tyrtyshnikov, A. Smirnov // Russian Journal of Numerical Analysis and Mathematical Modelling. — 2021. — Vol. 36, no. 3. — P. 177-181.

1.9. Личный вклад автора

В работе [23] автору принадлежат быстрый алгоритм вычисления произведения матрицы Якоби на вектор (не вошёл в текст диссертации, так как впоследствии оказался частным случаем быстрого алгоритма вычисления билинейного варианта оператора агрегации), конструкция предобусловли-вателя и лемма о единственном решении конечной однородной задачи, а также программная реализация описанного метода.

В работе [24] автором построено TTT-разложение для «скелета» оператора агрегации, и реализован полученный метод.

В работе [25] автором получен редуцированный вид уравнения Смолу-ховского и построен и реализован метод построения базиса редукции для частного класса ядер агрегации.

В работе [26] автором построен и реализован описанный алгоритм численного решения уравнения Смолуховского с переносом на основе редукции по пространству.

1.10. Основное содержание работы по главам

1. Введение.

2. Формальная постановка задачи. Обсуждение быстрого метода вычисления оператора Смолуховского и его обобщения на билинейную форму оператора.

3. Метод Ньютона-Крылова для поиска стационарного решения задачи. Обсуждение устройства матрицы Якоби оператора агрегации, построение предобусловливателя для неё. Демонстрация эффективности полученного метода.

4. Тот же метод Ньютона, но на основе тензорных разложений. Обсужде-

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

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

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

7. Заключение.

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

Автор выражает благодарность академику РАН Тыртышникову Евгению Евгеньевичу за научное руководство и поддержку в исследованиях, доценту факультета ВМК МГУ Матвееву Сергею Александровичу и доценту факультета ВМК МГУ Смирнову Александру Павловичу за неоценимую помощь и плодотворное сотрудничество, старшему научному сотруднику ИВМ РАН Замарашкину Николаю Леонидовичу и младшему научному сотруднику ИВМ РАН Желткову Дмитрию Александровичу за ценные замечания и обсуждение.

Глава 2

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

Рассматривается модель системы частиц, в изначальной формулировке пространственно-однородной, состоящей из частиц с кратными массами т, 2т, 3т,...; от нормировки в задаче ничего не зависит, так что в дальнейшем будем полагать т = 1, а массы — натуральными числами.

Поскольку система пространственно однородна, в каждый момент времени она полностью описывается концентрациями частиц различных масс, которые будем обозначать щ(£), к = 1, 2,...; из физических соображений, Пк(^ > 0. Требуется, таким образом, получить систему дифференциальных уравнений на концентрации, описывающую их эволюцию со временем. Для этого перечислим интересующие нас эффекты и взаимодействия между частицами.

Основным элементом нашей модели, который и рассматривался в основополагающей работе [8], будет слипание частиц; время от времени, частицы массой % и з сталкиваются и слипаются с образованием частицы массой % + з. Мы будем предполагать (полное обоснование выходит за рамки данной работы, и может быть найдено в той же [8]), что частота этого события пропорциональна произведению концентраций частиц масс % и з; коэффициент пропорциональности является первым параметром модели, и обозначается С{з. Из физических соображений ясно, что Cij = > 0.

Кроме того, будем предполагать, что запас частиц в системе пополняется со временем: а именно, что существует постоянный приток частиц массы к в объёме ^. В численных расчётах мы будем предполагать = 0 при к = 1 — то есть приток есть только для частиц единичной массы — однако большинство рассматриваемых алгоритмов никак не используют такой специальный вид при построении.

Наконец, мы также предусмотрим возможность деструктивных столкновений, при которых соударяющиеся частицы распадаются на частицы единичной массы. В общем случае такое предположение может показаться нереалистичным — в действительности спектр образующихся в результате распада частиц существенно более сложен [27], однако и такой относительно простой модели деструктивных столкновений оказывается во многих случаях достаточно для описания поведения системы в присутствие распада более общего вида [2, 28]. В общем случае коэффициент для частоты таких столкновений может быть не связан с коэффициентом для конструктивных столкновений Cij, но мы, следуя [2], будем полагать его равным АС^-, Л > 0.

Отметим, что в литературе существуют и другие модели фрагментации — например, самопроизвольный распад частиц [29, 30], существенно

*_» *_» и 1 1 и

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

Окончательно, правая часть уравнения для /(И представляет собой сумму слагаемых, отвечающих перечисленным эффектам; при этом для столкновений необходимо учесть как появление частиц в результате столк-

новении, так и расход на участие в них:

&пк_ 1 ^

— "к + 2 / у Сгк-гпгпк-г +

г=1

Д [к — 1]

+ [ 2= ] ^(г + 3 )Оц пгщ - (1 + Х)Пк^ Скгпг, (2.1)

1,0=1 г=1

где к — 1, 2,3 ....

Здесь, и в дальнейшем, мы используем обозначение [•] — скобку Айвер-сона [31] — для следующей функции

1, если Р,

[Р] Н (2.2)

0, иначе.

Конкретный частный случай её использования в (2.1) можно было бы записать и более привычно, через символ Кронекера 5к 1, однако

• скобка Айверсона всё равно понадобится нам в дальнейшем с несколько более сложными условиями, и

• в дальнейшем мы так же столкнёмся с равенствами между более сложными выражениями; на субъективный взгляд автора в таких ситуациях скобка Айверсона читается существенно легче, чем эквивалентная запись в терминах символа Кронекера,

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

Как уже отмечалось во введении, решение системы (2.1) напрямую затруднительно, по причине наличия бесконечного числа уравнений, поэтому на практике нас будет интересовать «обрезанная» система

Апк _ т , 1 V1 п

— "к + 2 / у Сгк-1П1Пк-1 +

1=1

Д [к — 1] N N

+--2- (г + 3)С*3пп - (1 + Л)п^ Скгпг, (2.3)

1,3=1 г=1

где к = 1, 2,...,Ы. Хотя ранее и говорилось, что такая система может в известной степени рассматриваться как модель физической системы с выводом достаточно массивных частиц, следует отметить, что сам механизм вывода достаточно грубый — все новообразованные частицы массы более N немедленно и бесследно исчезают из системы, в резком контрасте с никак не выводящимися частицами массы N — 1. По-видимому, систему (2.3) допустимо рассматривать как полноценную модель физической системы безо всяких оговорок насчёт «обрезанных» аппроксимаций лишь в случае, когда общее количество частиц массы более N в действительности исчезающе мало; пример такой системы мы увидим в дальнейшем.

2.1. Оператор правой части

Для целей анализа уравнения (2.3) нам будет удобно объединить зависящие от п слагаемые в правой части в оператор 5: ^ :

1 к—1 Д [к =1] ^ N

5(п)к = С^—¿ппк—--2-п*п — (1+А)пк^Скгпг. (2.4)

¿=1 ¿,7=1 ¿=1

Каждая компонента 5 (п) представляет собой квадратичную форму от вектора п, а потому мы может поставить ему в соответствие симметричный билинейный оператор, который мы так же обозначим 5: х ^ :

1 к—1 А [к = 1] N 5 (х,у)к = С*к-¿^Ук—г + -2- + 3 —

¿=1 ¿,.7=1

1+ \( N N \

--2^ ( Хк СкгУг + Ук ^ СкгжЛ . (2.5)

\ ¿=1 ¿=1 )

Здесь билинейность и симметрия понимаются в обычном смысле: 5 (ах +

ву, г) = а5(х, г) + в5(у, г) и 5(х, у) = 5(у, х).

2.2. Быстрый алгоритм вычисления

В [22] для специального класса ядер Cгj был предложен алгоритм вычисления оператора (2.4) сложности существенно меньшей, чем получаемые при вычислении «в лоб» O(N2). Он легко обобщается на билинейный оператор (2.5), что и будет проделано в данном разделе.

Предположим, что ядро Cij, рассматриваемое как матрица RNхN, имеет малый ранг Я. В этом случае для него можно, с учётом симметрии, записать сокращённое собственное разложение

я

Сгз ^ ^ и г адаиз а.1 (2.6)

а=1

где и — матрица ортонормированных собственных векторов, а да — собственные числа С. На практике его легко найти из любого скелетного разложения для С (например, если известно скелетное разложение С — ХУТ, и X — QЯ, то собственные числа да — собственные числа матрицы ЯУ Т Q размером Я х Я).

Подставив это выражение в (2.5), после элементарных алгебраических преобразований получаем

1 я к-1

£ (х,У)к — 2^2 ^^2 (хгиг а) (Ук-гик-га) + а=1 г=1 Д [к — 1] я N

+--2- ^2 Да ((гХгига) (Уз из а) + (Хг^а) (Ш и) а))

а=1 г,3=1

1 + Д....... .......... .

, хг,'

2

я / N N \

( (Хк ика)^2 (Угига) + (Ук ^а) ^ (Жг^а ) I , (2.7)

а=1 \ г=1 г=1 /

или, введя дополнительные обозначения

окончательно

Xk a xk Uk a

Xk о. kxk Uk а N

Xa ^ ^ Xk а k=1

N

xa ^ ^ xk а k=1

Д k- 1

yka = yk Uk а

1/ka = kyk Uka,

N

ya ^ ^ yk a k=1

N

ya ^ ^ yk a k=1

(2.8a) (2.8b)

(2.8c) (2.8d)

S(X,y)k = 2^ Ma^ ayk-г a +

a=1 г=1

+

Л [k = 1]

Д

^ ^ Ma (xaya + Xaya)

a=1

1 + Л 2

Д

У^ Ma (Xkaya + ykaXa) . (2.9)

a=1

Для вычисления Xa, Xa и аналогичных векторов для y по формулам (2.9) требуется O(NR) операций. После этого второе слагаемое теперь очевидно вычисляется за O(R) (поскольку оно не обращается в ноль только при k = 1); третье за O(NR) для всех k. Первое же слагаемое представляет собой сумму R свёрток вида, а потому вычисляется при помощи быстрого преобразования Фурье за O(NR log N) операций; такова же итоговая сложность вычисления всего оператора.

В рассматриваемых в следующих главах численных алгоритмах нередко приходится вычислять S(x,y) многократно для одного и того же x и различных у, или даже, как в главах 5 и 6, вычислять такие произведения для всех x и у из некоторого набора векторов. Чтобы понять, как такие вычисления можно организовать быстрее, и что на практике даёт переход к симметричному разложению (2.6), рассмотрим процесс вычислений по формуле (2.9) подробнее, предполагая, что свёртки вычисляются с помощью быстрого преобразования Фурье.

2

1. Для каждого из векторов x и y необходимо вычислить все вспомогательные векторы и матрицы из (2.8). Это требует O(NR) операций.

2. Для каждого столбца матриц xka и yka, должным образом расширенного, необходимо вычислить прямое преобразование Фурье — всего 2R преобразований Фурье с общей сложностью O(NR log N).

3. Необходимо перемножить полученные в предыдущем пункте образы Фурье столбцов; на это уйдёт ещё O(NR) операций.

4. От произведений образов Фурье необходимо теперь вычислить обратное преобразование Фурье; это ещё R преобразований с той же общей сложностью O(NR log N).

5. Наконец, если свёртки, фигурирующие в первом слагаемом, уже вычислены, остаток вычислений по формуле (2.9) «в лоб» требует ещё O(NR) операций.

Общая сложность при этом, как и обещалось, составляет O(NR log N) операций. Но теперь интересно посмотреть, как этот алгоритм изменится в некоторых практически важных частных случаях:

• Если x = y, первые два шага достаточно выполнить один раз вместо двух. Здесь как раз важна симметрия в (2.6): при использовании обычного скелетного разложения Cij = ^R=i UiaVja при x и y в (2.8) стояли бы разные матрицы, и даже при совпадающих x и y производные от них промежуточные векторы и матрицы оказались бы различными. В данном же случае такое совпадение позволяет сэкономить вычисления: если мерять сложность исключительно по наиболее затратной операции — преобразованию Фурье — то в таком случае достаточно вычислить лишь R прямых и R обратных преобразований, вместо 2R прямых и R обратных в общем случае — суммарный выигрыш в полтора раза.

• Если необходимо многократно вычислять 5(ж, у) с одним и тем же х и различными у, первые два шага для х достаточно выполнить один раз, после чего подставлять вычисленные векторы — включая образы Фурье — при каждом вычислении оператора. В этом случае достаточно будет выполнять первые два шага лишь для у, что, как и в предыдущем пункте, даёт ускорение примерно в полтора раза — хотя в этом случае и добавляются единовременные затраты на подготовку х.

• Наконец, в случае, если и х и у повторяются многократно — например, требуется для некоторого набора векторов вычислить все значения оператора на парах векторов (как в уравнениях (5.7) и (6.14)) — можно достичь ещё большего выигрыша. Предположим, что имеется п векторов х(г), и для каждой пары из них необходимо вычислить 5(х(г),х(^)). В этом случае можно сначала выполнить вычисления по первым двум шагам для каждого х(г) — пЯ преобразований Фурье — после чего для кажой пары векторов можно будет эти шаги пропустить вовсе, требуя, таким образом, лишь Яп(п + 1)/2 обратных преобразований Фурье для завершения вычислений. Для сравнения, при выполнении всего алгоритма для каждой пары векторов «с нуля» потребовалось бы суммарно 3Яп(п + 1)/2 прямых и обратных преобразований Фурье. В зависимости от п, таким образом, выигрыш от выделения из алгоритма подготовительной стадии составляет до трёх раз. Даже при условии выделения подготовительной части, алгоритм с обычным, несимметричным скелетным разложением требовал бы для каждого х(г) подготовки двух наборов промежуточных матриц и образов Фурье, что потребовало бы на подготовительном этапе 2пЯ преобразований вместо пЯ в предлагаемом варианте — выигрыш при больших п не столь заметный (поскольку всего преобразований 0(п2Я)), но тем не менее существенный.

Список литературы диссертационного исследования кандидат наук Тимохин Иван Владиславович, 2022 год

Литература

[1] Стояновская О.П., Снытников В.Н. Численное моделирование образования уединенных вихрей повышенной плотности в околозвездном диске // Вычислительные методы и программирование. — 2012. —Июнь. — Т. 13, № 33.-С. 377-383.

[2] Size distribution of particles in Saturn's rings from aggregation and fragmentation / Brilliantov N. V., Krapivsky P. L., Bodrova A., Spahn F., Hayakawa H., Stadnichuk V., and Schmidt J. // PNAS. — 2015.—Vol. 112, no. 31. —P. 9536-9541.

[3] Shrivastava R. C. A simple model of particle coalescence and breakup // J. Atom. Sci. —1982. —Vol. 39. —P. 1317.

[4] Blatz P. J., Tobolsky A. V. Note on the kinetics of systems manifesting simultaneous polymerization-depolymerization phenomena // The Journal of Physical Chemistry. —1945. —Vol. 49, no. 2. —P. 77-80.

[5] A numerical method for the simulation of an aggregation-driven population balance system / Hackbusch Wolfgang, John Volker, Khachatryan Aram, and Suciu Carina // International Journal for Numerical Methods in Fluids.—2012.—Aug.—Vol. 69. —P. 1646-1660.

[6] Ben-Naim E., Krapivsky P. L. Popularity-driven networking // Europhys. Lett. —2012.—Vol. 97, no. 4.—P. 48003.

[7] Slemrod Marshall. The Becker-Döring Equations // Modeling in Applied Sciences: A Kinetic Theory Approach / ed. by Bellomo Nicola, Pulvirenti Mario. — Boston, MA : Birkhäuser Boston, 2000. — P. 149-171. —ISBN: 978-1-4612-0513-5.

[8] Smoluchowski M. V. Drei vortrage uber diffusion, Brownsche bewegung und koagulation von kolloidteilchen // Zeitschrift fur Physik. — 1916. — Bd. 17. —S. 557-585.

[9] Smoluchowski M. V. Attempt for a mathematical theory of kinetic coagulation of colloid solutions // Z. Phys. Chem. — 1917.—Vol. 92. —P. 129.

[10] Müller H. Zur allgemeinen Theorie ser raschen Koagulation // Fortschrittsberichte über Kolloide und Polymere. — 1928. — Bd. 27, nu. 6. — S. 223-250.

[11] Эффективный метод численного решения математической модели переноса коагулирующих частиц / Загидуллин Р. Р., Смирнов А. П., Матвеев С. А. и Тыртышников Е. Е. // Вестник Московского университета. Серия 15. Вычислительная математика и кибернетика. — 2017. — № 4. — С. 28-34.

[12] Parallel Numerical Algorithm for Solving Advection Equation for Coagulating Particles / Matveev Sergey, Zagidullin Rishat, Smirnov Alexander, and Tyrtyshnikov Eugene // Supercomputing Frontiers and Innovations. — 2018. —Vol. 5, no. 2.

[13] Malyshkin L., Goodman J. The Timescale of Runaway Stochastic Coagulation // Icarus. —2001.—Vol. 150. —P. 314.

[14] Stadnichuk V., Bodrova A., Brilliantov N. V. Smoluchowski aggregation-fragmentation equations: Fast numerical method to find steady-state

solutions // International Journal of Modern Physics B. — 2015. —Vol. 29, no. 29. —P. 1550208.

[15] Lee M. H. On the validity of the coagulation equation and the nature of runaway growth // Icarus. — 2000. — Vol. 143, no. 1.—P. 74-86.

[16] Anderson acceleration method of finding steady-state particle size distribution for a wide class of aggregation-fragmentation models / Matveev S. A., Stadnichuk V. I., Tyrtyshnikov E. E., Smirnov A. P., Ampilogova N. V., and Brilliantov N. V. // Computer Physics Communications. — 2018. — Vol. 224. —P. 154-163.

[17] Stationary mass distribution and non-locality in models of coalescence and shattering / Connaughton C., Dutta A., Rajesh R., Siddharth N., and Zaboronski O. // Phys. Rev. E.— 2018.—Vol. 97. —P. 022137.

[18] Debry Edouard, Sportisse Bruno, Jourdain Benjamin. A stochastic approach for the numerical simulation of the general dynamics equation for aerosols // Journal of Computational Physics. — 2003. —Vol. 184, no. 2. — P. 649-669.

[19] Monte-Carlo Modeling of Aerosol Kinetics / Sorokin A. A., Strizhov V. F., Demin M. N., and Smirnov A. P. // Atomic Energy. — 2015.—Vol. 117, no. 4. —P. 289-293.

[20] Analysis of four Monte Carlo methods for the solution of population balances in dispersed systems / Zhao H., Maisels A., Matsoukas T., and Zheng C. // Powder Technology. — 2007.— Vol. 173, no. 1.—P. 38-50.

[21] Kolodko A., Sabelfeld Karl, Wagner W. A stochastic method for solving Smoluchowski's coagulation equation // Mathematics and Computers in Simulation. —1999. —July. —Vol. 49. —P. 57-79.

[22] Matveev S. A., Smirnov A. P., Tyrtyshnikov E. E. A fast numerical method for the Cauchy problem for the Smoluchowski equation // Journal of Computational Physics. — 2015. — Vol. 282. — P. 23-32.

[23] Newton method for stationary and quasi-stationary problems for Smoluchowski-type equations / Timokhin I. V., Matveev S. A., Siddharth N., Tyrtyshnikov E. E., Smirnov A. P., and Brilliantov N. V. // Journal of Computational Physics. — 2019. — Vol. 382. — P. 124-137.

[24] Timokhin Ivan. Tensorisation in the Solution of Smoluchowski Type Equations // International Conference on Large-Scale Scientific Computing / Springer. —2019. —P. 181-188.

[25] Метод поиска редуцированного базиса для нестационарных задач / Ти-мохин И. В., Матвеев С. А., Тыртышников Е. Е. и Смирнов А. П. // Доклады Российской академии наук. Математика, информатика, процессы управления. —2021.—Т. 497, № 2_21. С. 35-38.

[26] Model reduction for Smoluchowski equations with particle transfer / Timokhin I., Matveev S., Tyrtyshnikov E., and Smirnov A. // Russian Journal of Numerical Analysis and Mathematical Modelling. — 2021. — Vol. 36, no. 3. —P. 177-181.

[27] Âstrom J. A. Statistical models of brittle fragmentation // Advances in Physics. —2006. —Vol. 55, no. 3-4.—P. 247-278.

[28] Krapivsky P. L., Otieno W., Brilliantov N. V. Phase Transitions in Systems with Aggregation and Shattering // Phys. Rev. E. — 2017.—Vol. 96.— P. 042138.

[29] Agoshkov V. I., Dubovski P. B. Solution of the reconstruction problem of a source function in the coagulation-fragmentation equation // Russian

Journal of Numerical Analysis and Mathematical Modelling. — 2002. — Vol. 17, no. 4. —P. 319-330.

[30] Laurenzi I. J., Diamond S. L. Kinetics of random aggregation-fragmentation processes with multiple components // Phys. Rev. E. — 2003. — May. — Vol. 67. —P. 051103.

[31] Iverson Kenneth E. A Programming Language. — John Wiley & Sons, 1962.

[32] Tensor Train versus Monte Carlo for the multicomponent Smoluchowski coagulation equation / Matveev S. A., Zheltkov D. A., Tyrtyshnikov E. E., and Smirnov A. P. // Journal of Computational Physics. — 2016.—Vol. 316. —P. 164-179.

[33] Fast and Accurate Finite-difference Method Solving Multicomponent Smoluchowski Coagulation Equation with Source and Sink Terms / Smirnov A. P., Matveev S. A., Zheltkov D. A., and Tyrtyshnikov E. E. // Procedia Computer Science. — 2016.—Vol. 80. —P. 2141-2146.

[34] Тензорные разложения для решения уравнений математических моделей агрегации, допускающих многочастичные столкновения / Стефо-нишин Д. А., Матвеев С. А., Смирнов А. П. и Тыртышников Е. Е. // Выч. мет. программирование.— 2018. —Т. 19, № 4.— С. 390-404.

[35] Saad Y., Schultz M. H. GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems // SIAM Journal on Scientific and Statistical Computing. — 1986. — July. — Vol. 7, no. 3. — P. 856-869.

[36] Hager William W. Updating the Inverse of a Matrix // SIAM Review. — 1989. —Vol. 31, no. 2. —P. 221-239.

[37] Hayakawa H. Irreversible kinetic coagulations in the presence of a source // Journal of Physics A: Mathematical and General. — 1987. — Vol. 20, no. 12. —P. L801-L805.

[38] Oseledets I. V. Tensor-Train Decomposition // SIAM Journal on Scientific Computing. —2011.—Vol. 33, no. 5. —P. 2295-2317.

[39] Oseledets I. V. Approximation of 2d x 2d matrices using tensor decomposition // SIAM Journal on Matrix Analysis and Applications. — 2010. — Vol. 31, no. 4. —P. 2130-2145.

[40] Oseledets I. V., Tyrtyshnikov E. E. TT-cross approximation for multidimensional arrays // Linear Algebra and its Applications. — 2010.—Vol. 432, no. 1. —P. 70-88.

[41] Dolgov S. V., Savostyanov D. V. Alternating Minimal Energy Methods for Linear Systems in Higher Dimensions // SIAM Journal on Scientific Computing. —2014.—Vol. 36, no. 5. —P. A2248-A2271.

[42] Pinnau R. Model Reduction via Proper Orthogonal Decomposition // Model Order Reduction: Theory, Research Aspects and Applications. Mathematics in Industry / ed. by A. Schilders W. H., van der Vorst H. A., J. Rommes. — Springer, Berlin, Heidelberg, 2008. — Vol. 13.

[43] Gubisch Martin, Volkwein Stefan. Proper Orthogonal Decomposition for Linear-Quadratic Optimal Control // Model Reduction and Approximation / ed. by Benner Peter, Ohlberger Mario, Cohen Albert, Willcox Karen. — Society for Industrial and Applied Mathematics, 2017. — Computational Science & Engineering. — P. 3-64.

[44] Sirovich Lawrence. Turbulence and the Dynamics of Coherent Structures. I-III // Quart. Appl. Math. —1987.—Vol. 45, no. 3.—P. 561-590.

[45] Steady oscillations in aggregation-fragmentation processes / Bril-liantov N. V., Otieno W., Matveev S. A., Smirnov A. P., Tyrtyshnikov E. E., and Krapivsky P. L. // PHYSICAL REVIEW E. — 2018. — Vol. 98, no. 1.

[46] Collective Oscillations in Irreversible Coagulation Driven by Monomer Inputs and Large-Cluster Outputs / Ball R. C., Connaughton C., Jones P. P., Rajesh R., and Zaboronski O. // Physical Review Letters. — 2012.— Окт. —Vol. 109, no. 16.

[47] Периодические режимы в конденсирующейся аэродисперсной системе с источником / Загайнов В. А., Лушников А. А., Бахтырева М. С., Луценко А. О. и Ходжер Т. В. // Доклады Академии Наук. —2007. — Т. 414, № 1. —С. 93-96.

[48] Галкин В. А. Уравнение Смолуховского // Физматлит, Москва. — 2001.

[49] Harten Ami. High resolution schemes for hyperbolic conservation laws // Journal of Computational Physics. — 1983. —Vol. 49, no. 3. — P. 357-393.

[50] Sweby P. K. High Resolution Schemes Using Flux Limiters for Hyperbolic Conservation Laws // SIAM Journal on Numerical Analysis. — 1984. — Vol. 21, no. 5. —P. 995-1011.

[51] Van Leer Bram. Towards the ultimate conservative difference scheme. IV. A new approach to numerical convection // Journal of Computational Physics. — 1977.—Vol. 23, no. 3. —P. 276-299.

[52] Горейнов С. А., Тыртышников Е. Е., Замарашкин Н. Л. Псевдоскелетные аппроксимации матриц // Доклады РАН. — 1995. —Т. 343, № 2. — С. 151-152.

[53] Berenger Jean-Pierre. A perfectly matched layer for the absorption of elec-

tromagnetic waves // Journal of Computational Physics. — 1994.—Vol. 114, no. 2. —P. 185-200.