Журнал вычислительной математики и математической физики, 2022, T. 62, № 4, стр. 553-563

Общая структура редуцированных базисов для задач агрегационной кинетики различной размерности

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

¹ МГУ им. М.В. Ломоносова, факультет ВМК
119991 Москва, Ленинские горы, стр. 52, Россия

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

³ МГУ, Московский центр фундаментальной и прикладной математики
119991 Москва, Ленинские горы, Россия

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

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

EDN: DYSZEV
DOI: 10.31857/S004446692204010X

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

Аннотация

Данная работа посвящена исследованию структуры линейных пространств, позволяющих приближать решения систем дифференциальных уравнений, возникающих в задачах агрегационной кинетики с источником и стоком частиц, на больших интервалах времени с высокой точностью. Основным параметром для семейства рассматриваемых моделей агрегационной кинетики является их размерность N, соответствующая максимальному допустимому размеру частицы в системе. При больших значениях предельного допустимого размера частицы $N$ исследование таких задач оказывается затруднительным, так как начинает требовать значительных вычислительных ресурсов. Для ускорения расчетов мы используем идеологию редукции размерности при помощи поиска базиса искомого маломерного пространства методом снимков. В результате использования такого подхода можно установить существование такого пространства, построить его базис и свести исходную нелинейную задачу высокой размерности к задаче существенно меньшей размерности, позволяющей организовать вычисления более эффективно без существенной потери качества численных решений. Главным результатом данной работы является найденная возможность использования базиса задач большей размерности $N$ для решения задач меньшей размерности $M < N$ без подготовки и построения базиса для размерности $M$. В численных экспериментах мы демонстрируем, что достаточно просто взять первые $M < N$ компонент векторов базиса задачи размерности $N$ и использовать полученную систему в качестве базиса для новой задачи размерности $M$. С использованием такого базиса нам удается получать превосходную точность численных решений редуцированной задачи размерности $M$, хотя сами по себе решения задач размерности $M$ и $N$ не обязаны быть близки. Дополнительно в работе предлагается несколько оценок для норм решений уравнений необратимой агрегации с источником. Из оценок следует, что прямая подстановка векторов базиса $N$ в задачи размерности $M > N$ не приводит к успеху, а расширение базиса для задач большей размерности требует дополнительных усилий, что согласуется с результатами экспериментов. Библ. 28. Фиг. 3.

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

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

Процессы агрегации и фрагментации играют важную роль в огромном наборе природных явлений от роста полимерных структур [1]–[3], коагуляции тромбоцитов в крови [4], динамики аэрозолей и примесей в атмосфере [5] до образования звезд и планет [6]–[8]. В случае наличия источника мономеров классические уравнения Смолуховского принимают вид формально бесконечной системы дифференциальных уравнений:

(1.1)

$\frac{{d{{n}_{k}}}}{{dt}} = J{{\delta }_{{k1}}} + \frac{1}{2}\sum\limits_{i + j = k} {{{C}_{{ij}}}} {{n}_{i}}{{n}_{j}} - {{n}_{k}}\sum\limits_{i = 1} {{{C}_{{ik}}}} {{n}_{i}},\quad k = 1,\; \ldots ,\;\infty ,$

где

• ${{n}_{k}}$ – концентрация частиц массы $k$;

• ${{C}_{{ij}}}$ – “ядро”, характеризующее вероятность столкновения частиц размеров $i$ и $j$;

• $J$ – мощность источников частиц размера $1$.

Искомыми величинами здесь выступают концентрации ${{n}_{k}}$ частиц размера $k$ на единицу объема среды. При данной записи мы предполагаем, что рассматривается пространственно однородная система, частицы сталкиваются парно и в результате столкновений частиц размеров $i$ и $j$ возникает новая частица размера $i + j$, т.е. выполняется закон сохранения массы. Достаточные условия на коэффициенты ${{C}_{{ij}}}$ для выполнения закона сохранения массы можно найти, например, в работах [9], [10]. Асимптотические свойства стационарных решений уравнений Смолуховского с источником мономеров хорошо изучены в случае коэффициентов

(1.2)

${{C}_{{ij}}} = {{i}^{\nu }}{{j}^{\mu }} + {{i}^{\mu }}{{j}^{\nu }}$

в работе [11]. В данной работе мы сконцентрируемся на параметрическом семействе моделей агрегации с использованием конечного количества уравнений $N$, считая, что образующиеся агрегаты большего размера выпадают в осадок под действием внешних сил (например, гравитации), а в качестве основного ядра коагуляции возьмем частный случай (1.2):

${{C}_{{ij}}} = {{i}^{a}}{{j}^{{ - a}}} + {{i}^{{ - a}}}{{j}^{a}}$

с монодисперсными начальными условиями ${{n}_{k}}(0) = {{\delta }_{{k1}}}$.

В таких условиях исходную систему (1.1) естественно заменить схожей конечной системой

(1.3)

$\frac{{d{{n}_{k}}}}{{dt}} = J{{\delta }_{{k1}}} + \frac{1}{2}\sum\limits_{i + j = k} {{{C}_{{ij}}}} {{n}_{i}}{{n}_{j}} - {{n}_{k}}\sum\limits_{i = 1}^N {{{C}_{{ik}}}} {{n}_{i}},\quad k = 1,\; \ldots ,\;N.$

Данное семейство моделей вызывает высокий интерес исследователей и позволяет наблюдать ряд интересных эффектов при вариациях параметра размерности системы дифференциальных уравнений $N$: из асимптотически устойчивого стационарного положения система с ростом $N$ переходит в режим стационарных колебаний по времени [12], [13], что может иметь значение для анализа популяционной динамики в биологии [14]. Свойства смены характера решений на данный момент окончательно не исследованы, но ряд результатов для более простых моделей [15]–[17] позволяет сформировать гипотезу, что колебания в данных системах возникают в результате бифуркации Андронова–Хопфа.

Основной целью данной работы является исследование взаимосвязи решений системы (1.3) при различных $N$, имея в виду возможность использования информации об одном из них для быстрого построения других. Этот вопрос мы исследуем в контексте идеи редукции модели – аппроксимации многомерной задачи в пространстве меньшей размерности с соответствующим сокращением вычислительных расходов. Конкретное воплощение этой общей идеи следует из работ [18], [19], в которых для уравнения Смолуховского типа (1.3) показана возможность построения маломерной аппроксимации целиком в содержащем решение пространстве малой размерности, без необходимости обращаться к полной системе даже для промежуточных вычислений.

В общем случае редукция модели охватывает широкий спектр различных вычислительных техник с разной степенью привязки к конкретным задачам [20]–[24]; предложенный нами метод (см. алгоритм 1) основан на идее POD (Proper Orthogonal Decomposition) [25], состоящей в простейшем варианте в построении базиса искомого подпространства малой размерности из старших сингулярных векторов матрицы, составленной из некоторых предполагаемых представителей искомого маломерного пространства, которых, следуя [22], мы набираем из “снимков” решения $n({{t}_{k}})$ в равноотставленные моменты времени. Построение базиса таким методом требует, в общем случае, знания решения заранее (в классическом POD система предполагается зависящей от дополнительных параметров, и “снимки” берутся из решений при нескольких характерных значениях параметра; мы же заинтересованы в применении редукции к отдельно взятой системе). Для обхода этого ограничения в [18] предложен адаптивный алгоритм, строящий по начальному отрезку времени базис для решения в последующие моменты. Эффективность полученного метода проверена на задачах, имеющих циклическое по времени решение – на них показано ускорение вычислений в разы [19] – однако базовый факт наличия подпространства малой размерности, приближенно содержащего в себе решение, экспериментально подтверждается и в других случаях.

В то же время само построение базиса остается достаточно дорогостоящей процедурой, требующей к тому же предварительного построения значительной части искомого решения. Для пространственно неоднородных задач эту проблему можно частично обойти, если ограничиться редукцией по пространству и искать базис одновременно с решением; это было проделано в [26], и показало ускорение в десятки раз. Для пространственно однородной задачи мы предлагаем рассмотреть возможность переноса базиса между системами вида (1.3) с различными $N$. Для этого естественно задаться вопросом о связи решений систем разного размера, в том числе на достаточно больших промежутках времени (искомое линейное пространство определяется поведением решения в целом, не ограничиваясь каким-либо конечным временным интервалом). В этой связи мы

• экспериментально продемонстрируем существование существенных различий в долгосрочном поведении систем с различным $N$;

• получим оценку на величину этого расхождения;

• покажем, что, тем не менее, решения задач с меньшим $N$ с большой точностью лежат в линейной оболочке урезанных решений с большим $N$.

Алгоритм 1. Общая схема используемого алгоритма построения базиса редукции для (1.3). Подробное описание алгоритма см. в [18], [19]

Require: “Ширина” для “окон” $\tau > 0$, параметры точности $\varepsilon {\text{'}} > \varepsilon > 0$, $\delta > 0$

Ensure: $V$ – искомый базис.

$V \leftarrow 0 \in {{\mathbb{R}}^{{N \times 0}}}$

$k \leftarrow 0$

repeat

Решить исходную систему(1.3) полностью на интервале $[k\tau ,(k + 1)\tau ]$ любым способом.

$\hat {V} \leftarrow {\text{базис}}\;{\text{POD}}\;{\text{для}}$ n(t) на [kτ, (k + 1)τ]

$\Delta \leftarrow {{\left\| {(I - V{{V}^{{\text{т}}}})\hat {V}} \right\|}_{2}}$

if $\Delta > \varepsilon {\text{'}}$ then

$U\Sigma {{W}^{{\text{т}}}} \leftarrow (V\,|\,\hat {V})$ {сингулярное разложение с сингулярными числами ${{\sigma }_{i}}$}

Найти $i:{{\sigma }_{i}} > \delta ,\quad {{\sigma }_{{i + 1}}} \leqslant \delta $

$V \leftarrow {\text{первые}}$ i столбцов U

end if

$k \leftarrow k + 1$

until $\Delta \leqslant \varepsilon $

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

Кроме того, мы получим некоторые общие оценки на отличие решений уравнений (1.3) и (1.1). Полученные оценки и используемая техника отчасти схожи с первыми главами [9], с тем отличием, что мы исследуем именно дискретную формулировку уравнения Смолуховского и фокусируемся больше на конкретных численных оценках, а не качественных результатах.

2. ТЕОРЕТИЧЕСКИЕ ОЦЕНКИ

В данном разделе мы сконцентрируемся на общих оценках отличия решений систем уравнений (1.3) в терминах нормы ${{\left\| {\, \cdot \,} \right\|}_{1}}$.

2.1. Решения конечных систем

Прежде всего получим априорную оценку на степень расхождения решений (1.3) при различных $N$. Для этого введем ряд обозначений:

• ${\mathbf{n}}(t)$ – решение системы вида (1.3) с $M$ уравнениями; $n(t):{{\mathbb{R}}_{{ \geqslant 0}}} \to {{\mathbb{R}}^{M}}$;

• ${\mathbf{\tilde {n}}}(t)$ – решение такой же системы, но уже с $N > M$ уравнениями; $\tilde {n}(t):{{\mathbb{R}}_{{ \geqslant 0}}} \to {{\mathbb{R}}^{N}}$;

• ${{{\mathbf{\tilde {n}}}}_{{(1)}}}(t)$, ${{{\mathbf{\tilde {n}}}}_{{(2)}}}(t)$ – первые $M$ и последние $N - M$ компонент $\tilde {n}(t)$ соответственно; ${{{\mathbf{\tilde {n}}}}_{{(1)}}}(t):{{\mathbb{R}}_{{ \geqslant 0}}} \to {{\mathbb{R}}^{M}}$, ${{{\mathbf{\tilde {n}}}}_{{(2)}}}(t):{{\mathbb{R}}_{{ \geqslant 0}}} \to {{\mathbb{R}}^{{N - M}}}$;

• $C$ – матрица ядра $M \times M$;

• $\tilde {C}$ – матрица ядра $N \times N$; $\tilde {C} = \left( {\begin{array}{*{20}{c}} C&{{{{\tilde {C}}}_{{12}}}} \\ {\tilde {C}_{{12}}^{{\text{т}}}}&{{{{\tilde {C}}}_{{22}}}} \end{array}} \right)$.

Кроме того, нам будет полезно ввести в рассмотрение билинейный аналог квадратичного оператора из правой части (1.3):

${{\mathcal{S}}_{M}}(x,y) = \frac{1}{2}\sum\limits_{i + j = k} {{{C}_{{ij}}}} {{x}_{i}}{{y}_{j}} - \frac{1}{2}{{x}_{k}}\sum\limits_{i = 1}^M {{{C}_{{ik}}}} {{y}_{i}} - \frac{1}{2}{{y}_{k}}\sum\limits_{i = 1}^M {{{C}_{{ik}}}} {{x}_{i}}.$

Тогда уравнения на ${\mathbf{n}}$ и ${{{\mathbf{\tilde {n}}}}_{{(1)}}}$ можно переписать в векторном виде:

$\frac{{d{\mathbf{n}}}}{{dt}} = \mathcal{J} + {{\mathcal{S}}_{M}}({\mathbf{n}},{\mathbf{n}}),$

$\frac{{d{{{{\mathbf{\tilde {n}}}}}_{{(1)}}}}}{{dt}} = \mathcal{J} + {{\mathcal{S}}_{M}}({{{\mathbf{\tilde {n}}}}_{{(1)}}},{{{\mathbf{\tilde {n}}}}_{{(1)}}}) - \left( {{{{\tilde {n}}}_{k}}\sum\limits_{i = M + 1}^N {{{C}_{{ik}}}} {{{\tilde {n}}}_{i}}} \right)_{{k = 1}}^{M},$

где под $\left( {f(k)} \right)_{{k = 1}}^{M}$ понимается вектор из значений $f(k)$ при $k$ от 1 до $M$.

Вычитая теперь одно уравнение из другого, с учетом билинейности и симметрии ${{\mathcal{S}}_{M}}$, получаем

$\frac{{d{\mathbf{n}}}}{{dt}} - \frac{{d{{{{\mathbf{\tilde {n}}}}}_{{(1)}}}}}{{dt}} = {{\mathcal{S}}_{M}}({\mathbf{n}} - {{{\mathbf{\tilde {n}}}}_{{(1)}}},{\mathbf{n}} + {{{\mathbf{\tilde {n}}}}_{{(1)}}}) + \left( {{{{\tilde {n}}}_{k}}\sum\limits_{i = M + 1}^N {{{C}_{{ik}}}} {{{\tilde {n}}}_{i}}} \right)_{{k = 1}}^{M},$

откуда уже легко оценить сверху модуль разности

$\begin{gathered} \left| {\frac{{d{{n}_{k}}}}{{dt}} - \frac{{d{{{\tilde {n}}}_{k}}}}{{dt}}} \right| \leqslant \frac{1}{2}\sum\limits_{i + j = k} {{{C}_{{ij}}}\left| {{{n}_{i}} - {{{\tilde {n}}}_{i}}} \right|} ({{n}_{j}} + {{{\tilde {n}}}_{j}}) + {{{\tilde {n}}}_{k}}\sum\limits_{k = M + 1}^N {{{C}_{{ik}}}} {{{\tilde {n}}}_{i}} + \\ + \;\frac{1}{2}\left| {{{n}_{k}} - {{{\tilde {n}}}_{k}}} \right|\sum\limits_{k = 1}^M {{{C}_{{ik}}}} ({{n}_{i}} + {{{\tilde {n}}}_{i}}) + \frac{1}{2}({{n}_{k}} + {{{\tilde {n}}}_{k}})\sum\limits_{k = 1}^M {{{C}_{{ik}}}} \left| {{{n}_{i}} - {{{\tilde {n}}}_{i}}} \right|. \\ \end{gathered} $

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

$\begin{gathered} {{\left\| {\frac{{d{\mathbf{n}}}}{{dt}} - \frac{{d{{{{\mathbf{\tilde {n}}}}}_{{(1)}}}}}{{dt}}} \right\|}_{1}} \leqslant \frac{1}{2}\sum\limits_{i + j \leqslant M} {{{C}_{{ij}}}} \left| {{{n}_{i}} - {{{\tilde {n}}}_{i}}} \right|({{n}_{j}} + {{{\tilde {n}}}_{j}}) + \sum\limits_{i,j = 1}^M {{{C}_{{ij}}}} \left| {{{n}_{i}} - {{{\tilde {n}}}_{i}}} \right|({{n}_{j}} + {{{\tilde {n}}}_{j}}) + \\ + \;\sum\limits_{i = 1}^N {\sum\limits_{j = M + 1}^N {{{C}_{{ij}}}} } {{{\tilde {n}}}_{i}}{{{\tilde {n}}}_{j}} \leqslant \frac{3}{2}{{\left\| {{\mathbf{n}} - {{{{\mathbf{\tilde {n}}}}}_{{(1)}}}} \right\|}_{1}}\left\| {C\left( {{\mathbf{n}} + {{{{\mathbf{\tilde {n}}}}}_{{(1)}}}} \right)} \right\| + {\mathbf{\tilde {n}}}_{{(1)}}^{{\text{т}}}{{{\tilde {C}}}_{{12}}}{{{{\mathbf{\tilde {n}}}}}_{{(2)}}}. \\ \end{gathered} $

Последнее слагаемое в оценке представляет собой скалярное произведение и также легко оценивается:

${\mathbf{\tilde {n}}}_{{(1)}}^{{\text{т}}}{{\tilde {C}}_{{12}}}{{{\mathbf{\tilde {n}}}}_{{(2)}}} \leqslant {{\left\| {{{{{\mathbf{\tilde {n}}}}}_{{(1)}}}} \right\|}_{1}}{{\left\| {{{{\tilde {C}}}_{{12}}}{{{{\mathbf{\tilde {n}}}}}_{{(2)}}}} \right\|}_{\infty }}.$

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

$\frac{d}{{dt}}{{\left\| {{\mathbf{n}} - {{{{\mathbf{\tilde {n}}}}}_{{(1)}}}} \right\|}_{1}} \leqslant \frac{3}{2}{{\left\| {C\left( {{\mathbf{n}} + {{{{\mathbf{\tilde {n}}}}}_{{(1)}}}} \right)} \right\|}_{\infty }}{{\left\| {{\mathbf{n}} - {{{{\mathbf{\tilde {n}}}}}_{{(1)}}}} \right\|}_{1}} + {{\left\| {{{{{\mathbf{\tilde {n}}}}}_{{(1)}}}} \right\|}_{1}}{{\left\| {{{{\tilde {C}}}_{{12}}}{{{{\mathbf{\tilde {n}}}}}_{{(2)}}}} \right\|}_{\infty }},$

откуда, если элементы разности на отрезке $[{{t}_{0}},t]$меняют знак лишь конечное число раз,

${{\left\| {{\mathbf{n}} - {{{{\mathbf{\tilde {n}}}}}_{{(1)}}}} \right\|}_{1}} \leqslant \left\| {{\mathbf{n}}({{t}_{0}}) - {{{{\mathbf{\tilde {n}}}}}_{{(1)}}}({{t}_{0}})} \right\| + \int\limits_{{{t}_{0}}}^t {{{{\left\| {{{{{\mathbf{\tilde {n}}}}}_{{(1)}}}(\tau )} \right\|}}_{1}}} {{\left\| {{{{\tilde {C}}}_{{12}}}{{{{\mathbf{\tilde {n}}}}}_{{(2)}}}} \right\|}_{\infty }}d\tau + \frac{3}{2}\int\limits_{{{t}_{0}}}^t {{{{\left\| {C({\mathbf{n}} + {{{{\mathbf{\tilde {n}}}}}_{{(1)}}})} \right\|}}_{\infty }}} {{\left\| {{\mathbf{n}} - {{{{\mathbf{\tilde {n}}}}}_{{(1)}}}} \right\|}_{1}}d\tau ,$

а отсюда уже по лемме Гронуолла-Беллмана получаем

(2.1)

${{\left\| {{\mathbf{n}} - {{{{\mathbf{\tilde {n}}}}}_{{(1)}}}} \right\|}_{1}} \leqslant \exp \left( {{{{\left\| {{\mathbf{n}}({{t}_{0}}) - {{{{\mathbf{\tilde {n}}}}}_{{(1)}}}({{t}_{0}})} \right\|}}_{1}} + \frac{3}{2}\int\limits_{{{t}_{0}}}^t {{{{\left\| {C({\mathbf{n}} + {{{{\mathbf{\tilde {n}}}}}_{{(1)}}})} \right\|}}_{\infty }}} d\tau } \right)\int\limits_{{{t}_{0}}}^t {{{{\left\| {{{{{\mathbf{\tilde {n}}}}}_{{(1)}}}} \right\|}}_{1}}} {{\left\| {{{{\tilde {C}}}_{{12}}}{{{{\mathbf{\tilde {n}}}}}_{{(2)}}}} \right\|}_{\infty }}d\tau .$

Для дальнейшего оценивания норм в правой части нам пригодится следующая

Лемма 1. В ранее введенных обозначениях,

(2.2)

$\min \left( {{{{\left\| {1{\mathbf{n}}(0)} \right\|}}_{1}},\sqrt {\frac{J}{{\mathop {\max }\limits_{1 \leqslant i,j \leqslant N} {{C}_{{ij}}}}}} } \right) \leqslant {{\left\| {{\mathbf{n}}(t)} \right\|}_{1}} \leqslant \max \left( {{{{\left\| {{\mathbf{n}}(0)} \right\|}}_{1}},\sqrt {\frac{{2J}}{{\mathop {\min }\limits_{1 \leqslant i,j \leqslant N} {{C}_{{ij}}}}}} } \right).$

Доказательство. В силу неотрицательности ${{n}_{k}}(t)$, ${{\left\| {{\mathbf{n}}(t)} \right\|}_{1}} = \sum\nolimits_{k = 1}^N {{n}_{k}}(t)$, откуда

$\frac{d}{{dt}}{{\left\| {{\mathbf{n}}(t)} \right\|}_{1}} = \sum\limits_{k = 1}^N {\frac{{d{{n}_{k}}}}{{dt}}} = \frac{1}{2}\sum\limits_{k = 1}^N {\sum\limits_{i + j = k} {} {{C}_{{ij}}}} {{n}_{i}}{{n}_{j}} - \sum\limits_{k = 1}^N {\sum\limits_{i = 1}^N {{{C}_{{ik}}}} } {{n}_{i}}{{n}_{k}} + J.$

Первая из двух сумм равна просто $\sum\nolimits_{i + j \leqslant N} {{C}_{{ij}}}{{n}_{i}}{{n}_{j}}$; упрощая, получаем

(2.3)

$\frac{d}{{dt}}{{\left\| {{\mathbf{n}}(t)} \right\|}_{1}} = J - \frac{1}{2}\sum\limits_{\substack{ 1 \leqslant i,j \leqslant N \\ i + j \leqslant N } } {{C}_{{ij}}}{{n}_{i}}{{n}_{j}} - \sum\limits_{\substack{ 1 \leqslant i,j \leqslant N \\ i + j > N } } {{C}_{{ij}}}{{n}_{i}}{{n}_{j}}.$

Отсюда уже легко получить двустороннюю оценку на производную ${{\left\| {n(t)} \right\|}_{1}}$:

(2.4)

$J - \bar {C}\left\| {\mathbf{n}} \right\|_{1}^{2} \leqslant J - \sum\limits_{1 \leqslant i,j \leqslant N} {{{C}_{{ij}}}} {{n}_{i}}{{n}_{j}} \leqslant \frac{d}{{dt}}{{\left\| {\mathbf{n}} \right\|}_{1}} \leqslant J - \frac{1}{2}\sum\limits_{1 \leqslant i,j \leqslant N} {{{C}_{{ij}}}} {{n}_{i}}{{n}_{j}} \leqslant J - \frac{1}{2}\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{C} \left\| {\mathbf{n}} \right\|_{1}^{2},$

где $\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{C} \leqslant {{C}_{{ij}}} \leqslant \bar {C}$.

Таким образом, получаем оценки непосредственно на ${{\left\| {n(t)} \right\|}_{1}}$:

$f(t) \leqslant {{\left\| {{\mathbf{n}}(t)} \right\|}_{1}} \leqslant g(t),$

$f{\kern 1pt} {\text{'}}(t) = J - {{f}^{2}}(t)\bar {C},\quad f(0) = {{\left\| {{\mathbf{n}}(0)} \right\|}_{1}},$

$g{\kern 1pt} {\text{'}}(t) = J - {{g}^{2}}(t)\frac{1}{2}\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{C} ,\quad g(0) = {{\left\| {{\mathbf{n}}(0)} \right\|}_{1}}.$

Оба уравнения являются частным случаем $h{\text{'}}(t) = a - b{{h}^{2}}(t)$, решением которого является

$h(t) = \left( \begin{gathered} \sqrt {a{\text{/}}b} \operatorname{th} (t\sqrt {ab} + K),\quad h(0) < \sqrt {a{\text{/}}b} , \hfill \\ \sqrt {a{\text{/}}b} \operatorname{cth} (t\sqrt {ab} + K),\quad h(0) > \sqrt {a{\text{/}}b} , \hfill \\ \sqrt {a{\text{/}}b} ,\quad h(0) = \sqrt {a{\text{/}}b} \hfill \\ \end{gathered} \right.$

для некоторого $K$, так что $h(t)$ монотонно стремится к асимптоте $\sqrt {a{\text{/}}b} $.

Подставляя полученные оценки на $f(t)$ и $g(t)$, для $\left\| {{\mathbf{n}}(t)} \right\|_{1}^{{}}$ получаем

$\min \left( {{{{\left\| {{\mathbf{n}}(0)} \right\|}}_{1}},\sqrt {\frac{J}{{\overline C }}} } \right) \leqslant {{\left\| {{\mathbf{n}}(t)} \right\|}_{1}} \leqslant \max \left( {{{{\left\| {{\mathbf{n}}(0)} \right\|}}_{1}},\sqrt {\frac{{2J}}{{\underline C }}} } \right).$

Обозначим правую часть (2.2) $L$ и вернемся к уравнению (2.1):

${{\left\| {{\mathbf{n}} - {{{{\mathbf{\tilde {n}}}}}_{{(1)}}}} \right\|}_{1}} \leqslant \exp \left( {{{{\left\| {{\mathbf{n}}({{t}_{0}}) - {{{{\mathbf{\tilde {n}}}}}_{{(1)}}}({{t}_{0}})} \right\|}}_{1}} + \frac{3}{2}\int\limits_{{{t}_{0}}}^t {{{{\left\| C \right\|}}_{C}}} ({{{\left\| {\mathbf{n}} \right\|}}_{1}} + {{{\left\| {{{{{\mathbf{\tilde {n}}}}}_{{(1)}}})} \right\|}}_{1}}d\tau } \right){{\left\| {{{{\tilde {C}}}_{{12}}}} \right\|}_{C}}\int\limits_{{{t}_{0}}}^t {{{{\left\| {{{{{\mathbf{\tilde {n}}}}}_{{(1)}}}} \right\|}}_{1}}{{{\left\| {{{{{\mathbf{\tilde {n}}}}}_{{(2)}}}} \right\|}}_{1}}d\tau } ,$

где ${{\left\| {\, \cdot \,} \right\|}_{С}}$ – чебышёвская норма матрицы (максимум модуля элементов).

Поскольку ${{\left\| {{{{{\mathbf{\tilde {n}}}}}_{{(1)}}}} \right\|}_{1}} + {{\left\| {{{{{\mathbf{\tilde {n}}}}}_{{(2)}}}} \right\|}_{1}} \leqslant L$, ${{\left\| {{{{{\mathbf{\tilde {n}}}}}_{{(1)}}}} \right\|}_{1}}{{\left\| {{{{{\mathbf{\tilde {n}}}}}_{{(2)}}}} \right\|}_{1}} \leqslant {{L}^{2}}{\text{/}}4$, и окончательно получаем (предполагая, что минимум ${{C}_{{ij}}}$, а значит, и $L$, одинаковы для обеих систем) следующее утверждение.

Лемма 2. В предыдущих обозначениях

(2.5)

${{\left\| {{\mathbf{n}} - {{{{\mathbf{\tilde {n}}}}}_{{(1)}}}} \right\|}_{1}} \leqslant \frac{{{{{\left\| {{{{\tilde {C}}}_{{12}}}} \right\|}}_{С}}{{L}^{2}}}}{4}(t - {{t}_{0}})\exp \left( {{{{\left\| {{\mathbf{n}}({{t}_{0}}) - {{{{\mathbf{\tilde {n}}}}}_{{(1)}}}({{t}_{0}})} \right\|}}_{1}} + 3L\left\| C \right\|_{C}^{{}}(t - {{t}_{0}})} \right).$

В частном случае ядра ${{C}_{{ij}}} = {{i}^{\mu }}{{j}^{\nu }} + {{i}^{\nu }}{{j}^{\mu }}$, ${{C}_{{ij}}} \geqslant 2(ij{{)}^{{\frac{{\mu + \nu }}{2}}}}$; при $\mu + \nu \geqslant 0$ получаем ${{\left\| {n(t)} \right\|}_{1}} \leqslant \max \left( {\left\| {n(0)} \right\|,\sqrt J } \right)$, независимо от $N$. (Полученная в лемме 1 оценка снизу относительно менее полезна в этом случае: если $\nu $ или $\mu $ положительны, оценка стремится к 0 при $N \to \infty $.) Если, кроме того, $J = 1$ и $\sum\nolimits_k {{n}_{k}}(0) \leqslant 1$, получаем $L = 1$ и упрощенную оценку:

${{\left\| {{\mathbf{n}} - {{{{\mathbf{\tilde {n}}}}}_{{(1)}}}} \right\|}_{1}} \leqslant \frac{{{{{\left\| {{{{\tilde {C}}}_{{12}}}} \right\|}}_{С}}}}{4}(t - {{t}_{0}})\exp \left( {{{{\left\| {{\mathbf{n}}({{t}_{0}}) - {{{{\mathbf{\tilde {n}}}}}_{{(1)}}}({{t}_{0}})} \right\|}}_{1}} + 3{{{\left\| C \right\|}}_{C}}(t - {{t}_{0}})} \right).$

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

2.2. Решения конечной и бесконечной систем

Формально заменяя в выкладках предыдущего пункта $N$ на $\infty $ (переходить к пределу по $N$ пока, вообще говоря, нельзя, так как не установлена сходимость $n$), можно перенести полученные оценки на решение исходной системы (1.1), при условии равномерной по $t$ сходимости двойного ряда $\sum\nolimits_{i,j = 1}^\infty {{{C}_{{ij}}}} {{n}_{i}}(t){{n}_{j}}(t)$ и ограниченности ${{C}_{{ij}}}$:

$0 < \underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{C} \leqslant {{C}_{{ij}}} \leqslant \bar {C},\quad 1 \leqslant i,j.$

А именно, полагая теперь ${{\tilde {n}}_{k}}(t)$ решением системы (1.1), имеем следующие леммы.

Лемма 3. В вышеприведенных обозначениях,

$\min \left( {{{{\left\| {{\mathbf{\tilde {n}}}(0)} \right\|}}_{1}},\sqrt {\frac{{2J}}{{\mathop {\sup }\limits_{i,j} {{C}_{{ij}}}}}} } \right) \leqslant {{\left\| {{\mathbf{n}}(t)} \right\|}_{1}} \leqslant \max \left( {{{{\left\| {{\mathbf{\tilde {n}}}(0)} \right\|}}_{1}},\sqrt {\frac{{2J}}{{\mathop {\inf }\limits_{i,j} {{C}_{{ij}}}}}} } \right).$

Доказательство. Доказательство точно следует лемме 1, с тем лишь отличием, что в уравнении (2.3) пропадает последний член, а потому (2.4) принимает вид

$J - \frac{1}{2}\bar {C}\left\| {{\mathbf{\tilde {n}}}} \right\|_{1}^{2} \leqslant \left\{ {J - \frac{1}{2}\sum\limits_{i,j} {{{C}_{{ij}}}} {{{\tilde {n}}}_{i}}{{{\tilde {n}}}_{j}} = \frac{d}{{dt}}{{{\left\| {{\mathbf{\tilde {n}}}} \right\|}}_{1}}} \right\} \leqslant J - \frac{1}{2}\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{C} \left\| {{\mathbf{\tilde {n}}}} \right\|_{1}^{2},$

что и приводит к отличиям в результате.

Лемма 4. В вышеприведенных обозначениях

$\begin{gathered} {{\left\| {{\mathbf{n}} - {{{{\mathbf{\tilde {n}}}}}_{{(1)}}}} \right\|}_{1}} \leqslant \bar {C}\exp \left( {{{{\left\| {{\mathbf{n}}({{t}_{0}}) - {{{{\mathbf{\tilde {n}}}}}_{{(1)}}}({{t}_{0}})} \right\|}}_{1}} + 3L\bar {C}(t - {{t}_{0}})} \right)\int\limits_{{{t}_{0}}}^t {{{{\left\| {{{{{\mathbf{\tilde {n}}}}}_{{(1)}}}} \right\|}}_{1}}} {{\left\| {{{{{\mathbf{\tilde {n}}}}}_{{(2)}}}} \right\|}_{1}}d\tau \leqslant \\ \leqslant \;\frac{{\bar {C}{{L}^{2}}}}{4}(t - {{t}_{0}})\exp \left( {{{{\left\| {{\mathbf{n}}({{t}_{0}}) - {{{{\mathbf{\tilde {n}}}}}_{{(1)}}}({{t}_{0}})} \right\|}}_{1}} + 3L\bar {C}(t - {{t}_{0}})} \right). \\ \end{gathered} $

Доказательство. Доказательство (в оговоренном выше предположении о равномерной сходимости двойного ряда) в точности следует доказательству леммы 2, с заменой $N$ на $\infty $, так что ${\mathbf{\tilde {n}}}$ и ${{{\mathbf{\tilde {n}}}}_{{(2)}}}$ теперь последовательности; соответственно, ${{\left\| {{\mathbf{\tilde {n}}}} \right\|}_{1}} = \sum\nolimits_{i = 1}^\infty {{{{\tilde {n}}}_{i}}} $.

Если при этом известно поведение ${\mathbf{\tilde {n}}}(t)$, отсюда можно получить оценки на качество приближения полного решения конечным на конечном интервале времени; например, при ${{C}_{{ij}}} \equiv 1$ известно стационарное решение [11], в котором ${{\tilde {n}}_{k}} = O({{k}^{{ - 3/2}}})$. Если это же соотношение сохраняется для всех $t$, то ${{\left\| {{{{{\mathbf{\tilde {n}}}}}_{{(2)}}}} \right\|}_{1}} \leqslant K{{M}^{{ - 1/2}}}$ для некоторого $K > 0$, и из леммы 4 сразу получаем (с учетом $n(0) = {{{\mathbf{\tilde {n}}}}_{{(1)}}}(0)$)

${{\left\| {{\mathbf{n}} - {{{{\mathbf{\tilde {n}}}}}_{{(1)}}}} \right\|}_{1}} \leqslant \frac{1}{{\sqrt M }}t\exp \left( {3t\sqrt J } \right)LK,$

так что при любом фиксированном $t$ норма разности стремится к нулю с ростом $M$.

3. ЧИСЛЕННЫЕ ПРИМЕРЫ

В качестве примера рассмотрим систему (1.3) с ${{C}_{{ij}}} = (i{\text{/}}j{{)}^{a}} + {{(j{\text{/}}i)}^{a}}$. При $0 \leqslant a \leqslant 1{\text{/}}2$ решение исходной системы (1.1) стремится к стационарному, известному аналитически и воспроизводимому в конечной системе [27], [28]. При $a > 1{\text{/}}2$, однако, в решении конечной системы наблюдаются циклы; как мы вскоре покажем, в этом случае решения при разных $N$ могут существенно отличаться. В качестве конкретного примера такого значения параметра в численных экспериментах использовалось $a = 0.8$.

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

Фиг. 1.

Первая норма и общая масса решения (1.3) для разных $N$. В то время как общая масса решения весьма существенно растет (для бесконечной системы (1.1) она составляет $Jt$ в момент времени $t$), первая норма (суммарная концентрация частиц всех размеров) остается ограничена в соответствии с леммой 1.

На фиг. 2 показана относительная погрешность приближения ${{{\mathbf{\tilde {n}}}}_{{(1)}}}$ с помощью $n$ для различных значений $M$ и $N$ от ${{2}^{{12}}}$ до ${{2}^{{18}}}$. Погрешность демонстрирует циклический характер, обусловленный циклическим поведением решений, и экспоненциальный рост по обе стороны от пиков, в соответствии с общей формой теоретической оценки.

Фиг. 2.

Относительная погрешность приближения решения системы размера $N$ решением системы размера $M$ в первой норме, для различных пар $M$ и $N$; для каждой пары $M$ и $N$ сравниваются первые $\min (M,N)$ компонент решения (1.3) в один и тот же момент времени. Циклическое поведение решений приводит к циклическому поведению ошибки; рост ошибки вблизи пиков экспоненциальный, что качественно соответствует оценке (2.5).

Несмотря на экспоненциальный рост ошибки вблизи пиков, и общий рост их высоты (связанный, вероятно, с ростом периода и соответствующим ростом максимума $\left| {t - {{t}_{0}}} \right|$), стоит отметить, что при фиксированных $N$ и $t$ до первого пика ошибка падает с ростом $M$.

В [18] была показана возможность построения по начальному отрезку решения базиса линейного пространства, в котором лежит решение, и дальнейшего поиска решения в этом пространстве с меньшими вычислительными затратами. Мы можем рассмотреть возможность решения редуцированной таким образом системы размера $M$ в базисе, полученном для системы размера $N$ (должным образом усеченным или дополненным нулями и реортогонализованном) как непрямое сравнение общего поведения решений систем разного размера, не привязанное к конкретным моментам времени.

Результаты такого эксперимента приведены на фиг. 3 для $M$ и $N$ от ${{2}^{{11}}}$ до ${{2}^{{14}}}$. Как видно на графике, базис, полученный для системы большего размера, позволяет решать меньшие системы с хорошей точностью, т.е. решение меньшего размера с высокой точностью лежит в линейном пространстве, описываемом соответствующими компонентами большего решения, несмотря на высокие поточечные различия. В обратную сторону – от меньшего базиса к большей системе – переход существенно менее успешный; причина здесь, очевидно, в том, что такое редуцированное решение в действительности воспроизводит решение меньшей системы, дополненное нулями — а оно, как мы уже видели, существенно отличается от решения большей.

Фиг. 3.

Погрешность редуцированного решения системы размера $M$, построенного по базису, полученному для системы размера $N$. При $N > M$ векторы базиса обрезаются до размера $M$ и реортогонализуются; при $M > N$ они дополняются нулями до требуемого размера. Использование большего базиса для меньшей системы дает превосходную точность решения, несмотря на существенные отличия между самими решениями.

4. ЗАКЛЮЧЕНИЕ

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

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

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

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

Eigen M. Prionics or the kinetic basis of prion diseases // Biophysical chemistry. 1996. V. 63. № 1. P. A1–A18.
Krapivsky P.L., Redner S., Ben-Naim E. A Kinetic View of Statistical Physics. UK: Cambridge University Press, 2010.
Pöschel T., Brilliantov N.V., Frömmel C. Kinetics of prion growth // Biophysical Journal. 2003. V. 85. № 6. P. 3460–3474.
Garzon Dasgupta A.K., Martyanov A.A., Filkova A.A., Panteleev M.A., Sveshnikova A.N. Development of a Simple Kinetic Mathematical Model of Aggregation of Particles or Clustering of Receptors // Life. 2020. V. 10. № 6. P. 97.
Aloyan A.E., Arutyunyan V.O., Lushnikov A.A., Zagaynov V.A. Transport of coagulating aerosol in the atmosphere // J. of Aerosol Sci. 1997. V. 28. № 1. P. 67–85.
Stoyanovskaya O.P., Snytnikov V.N. Numerical modeling of formation of high density solitary vortices in a circumstellar disk // Vychislitel’nye Metody i Programmirovanie. 2012. V. 13. № 3. P. 377–383.
Esposito L. Planetary Rings. UK: Cambridge University Press, 2006.
Brilliantov N.V., Krapivsky P.L., Bodrova A., Spahn F., Hayakawa H., Stadnichuk V., Schmidt J. Size distribution of particles in Saturn’s rings from aggregation and fragmentation // PNAS. 2015. V. 112. № 31. P. 9536–9541.
Галкин В.А. Уравнение Смолуховского. М.: Физматлит, 2001.
Melzak Z.A. A scalar transport equation // Transactions of the American Mathematical Society. 1957. V. 85. № 2. P. 547–560.
Hayakawa H. Irreversible kinetic coagulations in the presence of a source // J. of Physics A: Math. and General. 1987. V. 20. № 12. P. L801–L805.
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 // Physical Review Letters. 2012. V. 109. № 16. P. 168304.
Esenturk E., Connaughton C. Role of zero clusters in exchange-driven growth with and without input // Physical Review E. 2020. V. 101. № 5. P. 052134.
Somka J., Stocker R. Bursts Characterize Coagulation of Rods in a Quiescent Fluid // Physical Review Letters. 2020. V. 124. № 25. P. 258001.
Pego R.L., Velázquez J.JL. Temporal oscillations in Becker–Doring equations with atomization // Nonlinearity. 2020. V. 33. № 4. P. 1812.
Budzinskiy S.S., Matveev S.A., Krapivsky P.L. Hopf bifurcation in addition-shattering kinetics // Physical Review E. 2021. V. 103. № 4. P. L040101.
Niethammer B., Pego R.L., Schlichting A., Velázquez J.JL. Oscillations in a Becker-D$\backslash $" oring model with injection and depletion. 2021. arXiv preprint arXiv:2102.06751.
Тимохин И.В., Матвеев С.А., Тыртышников Е.Е., Смирнов А.П. Метод поиска редуцированного базиса для нестационарных задач // Докл. АН. Математика, информатика, процессы управления. 2021. Т. 497. № 2_21. С. 35–38.
Timokhin I.V., Matveev S.A., Tyrtyshnikov E.E., Smirnov A.P. Model reduction in Smoluchowski-type equations. 2020. arXiv preprint arXiv:2008.09439.
Moore B. Principal component analysis in linear systems: Controllability, observability, and model reduction // IEEE Transactions on Automatic Control. 1981. T. 26. № 1. C. 17–32.
Glover K. Optimal Hankel-norm approximations of linear multivariable systems and their ${{l}_{\infty }}$-error bounds // Int. J. Control. 1984. T. 39. № 6. C. 115–193.
Sirovich L. Turbulence and the Dynamics of Coherent Structures. I–III // Quart. Appl. Math. 1987. T. 45. № 3. C. 561–590.
Petrov S. Model Order Reduction Algorithms in the Design of Electric Machines // Internat. Conference on Large-Scale Sci. Comput. 2019. P. 140–147.
Gusak J., Daulbaev T., Ponomarev E., Cichocki A., Oseledets I. Reduced-order modeling of deep neural networks // Comput. Math. and Math. Phys. 2021. V. 61. № 5. P. 774–785.
Pinnau R. Model Reduction via Proper Orthogonal Decomposition. Berlin: Springer, 2008.
Timokhin I.V., Matveev S.A., Tyrtyshnikov E.E., Smirnov A.P. Model reduction for Smoluchowski equations with particle transfer // Russian Journal of Numerical Analys. and Math. Modelling. 2021. V. 36. № 3. P. 177–181.
Matveev S.A., Stadnichuk V.I., Tyrtyshnikov E.E., Smirnov A.P., Ampilogova N.V., Brilliantovy N.V. Anderson acceleration method of finding steady-state particle size distribution for a wide class of aggregation–fragmentation models // Comput. Phys. Communicat. 2018. V. 224. P. 154–163.
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. of Comput. Phys. 2019. V. 382. P. 124–137.

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

Инструменты

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

Журнал вычислительной математики и математической физики

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