Пермский государственный гуманитарно-педагогический университет
К вопросу о численном расчете пространственно-распределенных динамических систем с запаздыванием по времени
Д.А. Брацун
Введение
В последнее время все больший интерес вызывают динамические системы, имеющие наследственные слагаемые. Состояние такой динамической системы в какой-то фиксированный момент времени может зависеть в общем случае от всей предыдущей эволюции системы или определяться только конкретным моментом в прошлом. В последнем случае возникает динамическая система, зависящая от запаздывающего аргумента. Диапазон приложений теории таких систем включает в себя популяционную динамику и социальные процессы [1-3], нелинейные химические реакции [4] и процессы генной регуляции [5], поведение систем с автоматическим управлением [6, 7], механику жидкости [8-10] и т.д. Запаздывание может быть обусловлено самыми различными причинами, например, растянутостью изучаемого процесса во времени (воспроизводство популяции в биологии), ограниченностью скорости распространения сигнала (электромагнитная волна в релятивистской электродинамике), наличием инерционности некоторых элементов (в теории управлении с обратной связью) или существованием цепочек многоэтапных последовательных реакций с известным результатом в конце (процессы транскрипции-трансляции в математической генетике). В ряде случаев исключение запаздывания из рассмотрения эволюции системы позволяет относительно адекватно описывать реальные процессы. Однако зачастую такое упрощение приводит к парадоксальным выводам [11]. Необходимо отметить, что уравнения, содержащие члены с запаздывающими аргументами, анализировать непросто. В подавляющем числе случаев подобные задачи не удается решить аналитически, а численное решение также вызывает определенные трудности.
Среди моделей с запаздыванием особое место занимают пространственно распределенные системы, так как они дают возможность изучать влияние запаздывания на структурообразование в системе. Вообще говоря, наиболее распространённым методом расчета нелокальных систем является конечно-разностный метод. Простейшая явная схема позволяет организовать расчет с привлечением пространственных полей всего для двух моментов времени. Особенность нелокальной системы с запаздыванием заключается в том, что для нахождения значения функций с запаздывающим аргументом необходимо хранить данные полей всех временных слоев в пределах диапазона запаздывания. На практике такой подход реализуем только для сравнительно малоразмерных моделей с небольшим временем запаздывания. В противном случае объем оперируемых данных, а значит, и время расчета модели значительно возрастают. Оценки показывают, что даже для двумерных эволюционных моделей прямолинейное хранение всех временных слоёв требует настолько значительных объемов оперативной памяти, что это сильно затрудняет использование персональных компьютеров. Применение конечно-разностных схем с переменным шагом по времени или расчёт эволюции трёхмерных пространственных моделей только усугубляют ситуацию.
В данной работе предлагается простой алгоритм, который позволяет запоминать в рамках диапазона запаздывания только некоторые пространственные поля, а недостающую информацию интерполировать в ходе проведения расчета. Эффективность метода демонстрируется на примере численного моделирования процессов пространственно-временной динамики белков, ответственных за поддержание биоритмов в клетках.
1. Метод расчета систем реакции-диффузии с запаздыванием по времени
Наиболее распространенным способом построения моделей реакции-диффузии, включающих одновременно запаздывание и диффузию, является простое добавление соответствующего диффузионного слагаемого к обыкновенному дифференциальному уравнению. Этот подход получил свою критику в литературе [4], где указывается, что диффузия и запаздывание не являются независимыми друг от друга явлениями. Однако эта проблема не является принципиальной в рамках данной статьи, так как предложенные способы корректировки моделей не меняют общей конструкции уравнений и не устраняют технические трудности, возникающие при изучении этих уравнений.
Таким образом, мы построим модель реакции-диффузии с запаздыванием традиционным способом. Рассмотрим систему нелинейных кинетических уравнений диффузионного типа
,
где - набор физических (или химических) переменных, которые определяют нелинейную кинетику системы в отсутствие диффузии; - матрица коэффициентов пространственной диффузии; Д - оператор Лапласа в декартовых координатах , и . Нелинейная функция содержит как обычные, так и запаздывающие во времени величины. Время запаздывания ф для простоты фиксировано, хотя для предлагаемого метода расчёта это не принципиально. Система уравнений (1) должна быть дополнена граничными и начальными условиями.
Численное исследование динамических систем вида (1), как правило, осуществляется с помощью метода конечных разностей, который является одним из наиболее распространенных численных методов решения уравнений в частных производных. При построении разностных схем заметным преимуществом обладают явные схемы, отличающиеся логической простотой в сочетании с их надёжностью. Однако в сравнении с неявными схемами они обладают одним существенным недостатком, а именно - жестким ограничением на шаг по времени для обеспечения устойчивости схем. Тем не менее, популярность явных схем среди исследователей со временем только растет, так как целью современного научного анализа изучаемой динамической системы все больше становится не просто поиск ее стационарных решений, но и изучение нелинейной динамики системы. Кроме того, появление высокопроизводительных вычислительных систем, позволило серьезно компенсировать ограничение явных методов. И все-таки одними из ключевых проблем, которые возникают при компьютерном моделировании задач математической физики, являются несовершенство архитектуры вычислительной техники, выражающееся в задержке ("лаге") при передаче данных от процессора к оперативной памяти, и небольшой размер кэш-памяти процессора. Указанные проблемы в полной мере проявляют себя в задачах с запаздыванием, в которых, во-первых, на каждом шаге по времени используются запаздывающие значения функции, а во-вторых, для этого необходимо хранить большой объем данных внутри диапазона времени запаздывания.
Введем в рассмотрение равномерную пространственную сетку
, ; ;
, ,
где - шаг сетки; , , - количество узлов вдоль координат , , соответственно. Запишем систему уравнений (1) в конечно-разностной форме, используя центральные разности для пространственных производных и разность вперед для производной по времени
где и - шаг по времени, - число шагов по времени, выполненных в рамках диапазона времени запаздывания ф. Как видно из выражения (3), для вычисления поля на шаге по времени требуется информация о поле не только в предыдущий момент времени , но и в момент в прошлом. Если время запаздывания ф велико, то велико и число . Необходимо сохранять в памяти данные для всех временных слоев в диапазоне времени запаздывания, так как они участвуют в построении решения. Ситуация усугубляется тем, что для явной схемы очень важно контролировать шаг по времени для устойчивости схемы он должен быть меньше определенного значения (критерий Куранта)
.
Допустим, что размер интегрируемой области l=50 и рассматривается однородная сетка . Если время запаздывания в системе , то согласно критерию Куранта должно выполняться условие (при условии , что очень часто встречается в задачах реакции-диффузии). В этом случае в диапазоне времени запаздывания укладывается не менее 80 фиксированных шагов по времени. Если каждый элемент сеточной функции определяется вещественным числом, занимающим 8 байт памяти, то для хранения данных потребуется выделить порядка 610,4 Мб памяти. Учитывая обстоятельство, что запаздывающие значения используются на каждом шаге по времени, организация хранения данных на внешнем носителе в данном случае является не допустимой, так как это сопровождается значительными затратами времени на доступ к этим устройствам. Выгрузка данных в кэш-память процессора также невозможна, так как эта память имеет ограниченный объем. Так, на момент написания статьи максимальный размер кэш-памяти даже третьего уровня (L3-cache) составляет немногим более 24 Мбайт. Таким образом, ключевым фактором, влияющим на быстроту проведения расчетов, является выгрузка в оперативную память больших объемов данных на каждом шаге по времени. Считывание из оперативной памяти также занимает определенное время.
Если компоненты вектора в ходе эволюции испытывают значительные колебания, то шаг по времени согласно критерию (4) также меняется. Здесь возникает следующая проблема если на промежутке времени запаздывания целое число сделанных шагов не укладывается, то в момент просто не существует поля , которое могло бы быть использовано для построения текущего поля. В этом случае необходимо либо сгенерировать это недостающее поле путем интерполяции по двум ближайшим соседям, либо отказаться от плавающего шага по времени.
Для преодоления всех этих трудностей мы предлагаем следующее решение данные хранить не на каждом шаге по времени, а только в определенные моменты, которые будем называть опорными. Значения переменных для промежуточных слоев времени восстанавливаются путем интерполяции данных опорных слоёв, хранимых в памяти (рис. 1). Предполагается, что в зависимости от параметров конкретной модели, гладкости функций и размерности накладываемой сетки могут быть использованы различные виды интерполяции.
Рис. 1 Схема выделения опорных временных слоев (отмечены аббревиатурой BL) - текущее время расчета, - время запаздывание, - шаг по времени
Вторая идея алгоритма заключается в использовании оптимального количества опорных слоев, т. е. предполагается, что в процессе расчета плотность расположения таких слоев в пределах диапазона запаздывания может меняться. Частота выделения опорных слоев определяется гладкостью функции - это позволит фиксировать точность вычислений. Данный подход также позволяет снизить погрешность в системах с динамически меняющимся шагом по времени.
Для оценки гладкости функции введем функцию , которая определяет относительный характер максимального изменения значения поля. промежуточный поле пространственный опорный
,
где - компонента вектора .
Далее введем в рассмотрение важный параметр алгоритма - число , определяющее сколько и каких слоев на интервале запаздывания должны стать опорными и заноситься в оперативную память. Оно означает, что в потоке данных каждый -й слой является опорным. В случае, когда моделируемая функция достаточно быстро меняет свое значение, для максимально точного расчета запаздывающих значений необходимо использовать все данные в пределах диапазона запаздывания . Фактически это значение параметра соответствует отказу от применения алгоритма и прямолинейному запоминанию всего потока данных. Для определенного класса задач, один пример из которых будет представлен ниже, эмпирически установлено, что этот случай наступает при . При меньших значениях алгоритм можно применять без существенной потери точности. Для минимизации погрешности интерполяции и сохранения достаточно высокой скорости расчета параметр должен принимать значения из отрезка [1, 20]. Наибольший положительный эффект применения данного алгоритма проявляется, если функция меняется достаточно медленно () для вычисления запаздывающего значения в этом случае можно ограничиться каждым двадцатым временным слоем (), сокращая тем самым объем хранимых данных и количество операций с этими данными. При больших значениях параметра удовлетворить заданной точности расчетов становится сложнее. Таким образом, для нахождения оптимального числа опорных слоев будем пользоваться следующим эмпирическим правилом
.
Соотношения (6) позволяют задать степень детализации данных полей, хранимых в памяти компьютера, в зависимости от скорости изменения функции , которая определяется переменной .
Рис. 2. Зависимость количества опорных слоев от гладкости функции
При этом параметр не зависит от шага по времени, позволяя избавиться от лишних данных, когда функция меняется слабо (рис. 2).
Очевидно, что в рамках предлагаемого метода возникает проблема восстановления значений полей для всех временных слоёв, данные для которых не запоминаются в ходе расчета. Учитывая, что число опорных слоев в потоке данных может варьироваться, для получения запаздывающих значений применялся интерполяционный полином Ньютона (7). Процедура производится для каждого узла сетки, при этом узлами интерполирования являются значения всех хранимых опорных слоев (рис. 3).
Рис. 3. Схема интерполяции требуемых данных по хранимым в памяти данным опорных слоев
Необходимо заметить, что восстанавливаемый временной слой всегда находится между первым и вторым опорными слоями (см. рис. 3). Если функция достаточно гладкая, при построении полинома можно использовать только эти опорные слои, сокращая время расчета. Но в общем случае для N опорных слоев многочлен Ньютона имеет вид
,
где R - разделенные разности соответствующего порядка, - время опорного слоя N, - компонента вектора для слоя N.
2. Пример реализации метода расчет эволюции циркадианных колебаний
В работе [12] нами была предложена модель, описывающая пространственно-временную динамику белков, ответственных за механизм циркадианных колебаний в организме Neurospora crassa. Главным элементом механизма колебаний является эффект запаздывания в биохимических реакциях транскрипции генов. Молекулярные процессы, которые приняты во внимание в нашей модели, схематически представлены на рис. 4.
Рис. 4. Схема процессов транскрипции генов frq и wc-1
Изучение генома нейроспоры выявило несколько генов, которые участвуют в регулировании суточных ритмов этого организма. В связи с этим мы выделяем две основные компоненты механизма колебаний. Одной из них является ген frq, который кодирует белок FRQ (отмечен на рис. 4 белым цветом). Как известно, в середине темного времени суток уровень FRQ низок, зато к середине светлого времени суток он достигает своего максимума [13]. Второй компонент формируют два гена wc-1 и wc-2, которые кодируют соответствующие белки WC-1 и WC-2. Изучение поведения этих белков показало, что они формируют гетеродимерный комплекс WCC [14], который мы рассматриваем в качестве партнера для белка FRQ (отмечен на рис. 4 черным цветом).