Идентификация параметров осциллирующих процессов в живой природе, моделируемых дифференциальными уравнениями

Санкт-Петербургский Государственный Университет Реферат

Идентификация параметров осциллирующих процессов в живой природе, моделируемых дифференциальными уравнениями

Выполнила студентка 312гр.

Варламова А.А.

Проверил Токин И. Б Санкт-Петербург

1. Идентификация параметров в системах описываемых ОДУ

1.1 Градиентные уравнения

1.2 Уравнения в вариациях

1.3 Функционалы метода наименьших квадратов

1.4 Численное решение градиентных уравнений

1.4.1 Полиномиальные системы

1.4.2 Метод рядов Тейлора

1.4.3 Метод Рунге-Кутта

2. Модели осциллирующих процессов в живой природе

2.1 Модель Лотки

2.1.1 Осциллирующие химические реакции

2.1.2 Осцилляция популяций в системе «хищник-жертва»

2.2 Другие модели

3. Идентификация параметров модели Лотки

3.1 Дифференциальные уравнения

3.2 Постановки задачи идентификации и функционалы МНК

3.3 Как ускорить вычисления

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

4. О других методах идентификации Литература

1. Идентификация параметров в системах, описываемых ОДУ

1.1 Градиентные уравнения

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

Рассмотрим вещественнозначную функцию аргумента, и пусть и. Тогда величина

(1)

то есть производная функции по направлению характеризует скорость изменения при изменении в направлении вектора .

Из формулы (1) получаем:

(2)

где — градиент функции, а это дает:

(3)

(4)

(5)

Таким образом, вектор является направлением наискорейшего рос-та функции в точке, а вектор — это направление наискорейшего ее убывания в этой точке.

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

Это означает, что удовлетворяет дифференциальному уравнению:

(6)

или в координатной форме:

(7)

К уравнениям (6) или (7) добавляем начальные условия:

(8)

или в координатной форме:

(9)

Решение задачи Коши (6),(8) (или (7),(9)) определяет градиентную кривую проходящую через точку. Будем рассматривать это решение как век-тор-функцию аргументов и .

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

Если решение задачи (6),(8) существует при, то при любом та-ком получаем, что:

при (11)

при (12)

и мы вправе ожидать, что

(13)

Метод градиентных уравнений нахождения локального минимума функции заключается в численном интегрировании задачи Коши (6),(8) вдоль оси до достижения точки, достаточно близкой к .

1.2 Уравнения в вариациях

Рассмотрим задачу Коши:

(14)

(15)

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

(16)

удовлетворяют следующей задаче Коши:

(17)

(18)

Уравнения (17) относительно производных (16) называют уравнениями в вариациях для уравнений (14).

1.3 Функционалы метода наименьших квадратов

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

(19)

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

В методе наименьших квадратов нахождения (идентификации) параметров рассматривают функционал

(20)

где — фиксированные весовые коэффициенты, а — значения первых компонент решения задачи (14),(15) в точке при заданных

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

Для того, чтобы воспользоваться методом градиентных уравнений, необходимо выписать уравнения (7) для функционала (20):

(21)

Эти градиентные уравнения надо дополнить начальными условиями:

(22)

1.4 Численное решение градиентных уравнений

Обратимся к функционалу, , определенному в п. 1.3. Пря-мой способ нахождения приближенного значения точки, определенной по формуле (17) (то есть точки предполагаемого минимума функционала), — это численное интегрирование градиентных уравнений (21) при начальных условиях (22).

Правые части уравнений (21) зависят от неизвестных через значения функций в точках при, ,. При фиксированных значениях величины могут быть получены численным интегрированием уравнений (14),(17) при начальных условиях (15),(18).

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

(23)

(24)

отправляясь от значения, последовательно получать приближенные значения решения в точках

Числа называют шагами интегрирования, а числа ,…- узлами таблицы или сетки численного интегрирования. Совокупность узлов называют сет-кой, а величины называют значениями решения на узлах сетки. Если то говорят о равномерной сетке или об интегрировании с постоянным шагом.

Численное интегрирование градиентных уравнений, как правило, требует частой смены величины шага интегрирования. Хорошо к быстрой смене шага приспособлены явные методы Рунге-Кутта и метод рядов Тейлора.

Пошаговые методы численного интегрирования обыкновенных дифференциальных уравнений хорошо освещены в литературе по численному анализу (см., например, [2,3]).

1.4.1 Полиномиальные системы

Полиномиальной системой мы будем называть автономную систему ОДУ

(25)

где — алгебраические полиномы по .

Какие системы ОДУ можно свести к полиномиальным и как это делается? Начнем с примера. Рассмотрим задачу Коши:

(26)

(27)

Вводя дополнительные переменные

(28)

получаем следующую квадратичную задачу Коши:

(29)

(30)

Теперь рассмотрим достаточно общий случай. Рассмотрим класс сис-тем ОДУ (23), правые части которых можно представить в виде:

(31)

где все функции, а также все функции

(32)

являются алгебраическими полиномами по .

Любая система из сводится к полиномиальной. Действительно, если в (23),(24) ввести дополнительные переменные то:

(33)

(34)

где все правые части

(35)

— алгебраические полиномы по с постоянными коэффициентами.

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

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

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

1.4.2 Метод рядов Тейлора

Введем в рассмотрение оператор, сопоставляющий решению задачи Коши (23), (24) его полином Тейлора

(36)

порядка. Радиус сходимости ряда обозначим .

Метод рядов Тейлора решения задачи Коши (23), (24) заключается в построении таблицы приближенных значений по формулам:

, (37)

где — натуральные,, , а удовлетворяют неравенствам .

Для программной реализации метода рядов Тейлора необходимы алгоритмы нахождения коэффициентов Тейлора и автоматического выбора величины шага интегрирования.

Нахождение коэффициентов Тейлора

Рассмотрим квадратичную задачу Коши

(38)

(39)

где — вещественные или комплексные постоянные, а — вещественная или комплексная переменная.

Подставляя в (38) разложение Тейлора

(40)

получаем:

(41)

Приводя подобные члены и приравнивая все коэффициенты полученного степенного ряда нулю, получаем искомые формулы:

;

, , (42)

где, .

Аналогичные формулы легко вывести и для общего случая полиномиальной системы степени .

Оценка погрешности и выбор шага

Рассмотрим полиномиальную задачу Коши:

(43)

(44)

где, ,, а максимальная степень полиномов (степень системы (43)) равна .

Введем обозначения:

, (45)

и будем предполагать, что .

Теорема.

Решение задачи (43), (44) голоморфно в круге и удовлетворяет там неравенствам:

(46)

где

, (47)

Используя эту теорему несложно построить алгоритм автоматического выбора шага в методе рядов Тейлора по заданной пользователем границе абсолютной (или относительной) погрешности.

1.4.3 Метод Рунге-Кутта

Этим методам посвящено много работ, и они хорошо изложены в много-численных учебниках (см., например, [2,3]).

2. Модели осциллирующих процессов в живой природе

2.1 Модель Лотки

2.1.1 Осциллирующие химические реакции

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

Рассматривается математическая модель взаимодействия на молекулярном уровне веществ на основе следующих предположений:

1. При взаимодействии с молекулой вещества молекула вещества превращается в молекулу вещества. Это описывают в форме молекулярной ре-акции:

(1)

Такую реакцию относят к классу автокаталитических, так как наличие вещества обеспечивает превращение другого вещества в .

2. При взаимодействии с молекулой вещества молекула вещества пре-вращается в молекулу вещества, то есть происходит автокаталитическая молекулярная реакция:

(2)

3. Вещество в то же время необратимо распадается, превращаясь в вещество, то есть происходит молекулярная реакция

(3)

4. Скорости протекания реакций (1), (2), (3) пропорциональны концентрациям веществ в левых частях этих реакций, то есть равны соответственно:

, , (4)

где символами, , обозначены концентрации веществ, , со-ответственно, а коэффициенты — положительные числа.

5. Скорость изменения концентрации каждого вещества равна сумме скоростей изменения концентраций этого вещества во всех реакциях, в которых оно участвует.

Из условий 1−5 следуют равенства:

(5)

где — концентрация вещества. Это система ОДУ Лотки.

2.1.2 Осцилляция популяций в системе «хищник-жертва»

Первая экологическая модель типа «хищник — жертва» была предложена в книге Лотки. Она основана на тех же уравнениях (5).

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

1. Наличие зайцев и еды для них приводит к увеличению количества зайцев, что можно записать формулой:

(6)

2. Наличие волков и еды для них приводит к увеличению количества волков:

(7)

3. Волки умирают от болезней или старости:

(8)

4. Скорость изменения количества зайцев по формуле (6), скорость изменения количества волков по формуле (7) и скорость увеличения количеств умерших волков по формуле (8) равны соответственно:

, , (9)

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

5. Скорость изменения каждого из количеств (количество умерших волков) равна сумме скоростей изменения этих количеств в каждом из процессов (6), (7), (8), в котором соответствующая величина участвует.

Из условий 1−5 следуют уравнения Лотки (5), только символы имеют другой смысл.

Более общие модели поведения хищников и жертв в различных эко-логических ситуациях были предложены в лекциях Вольтерры. В связи с этим, уравнения Лотки (5) называют часто уравнениями Лотки-Вольтерра.

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

Сведение модели (5) к двумерной основано на предположении, что вели-чина постоянна. В случае модели осциллирующих химических реакций это означает, что вещества достаточно много, а в случае модели «хищник — жертва» это означает, что еды у зайцев достаточно много. Из этого предполо-жения следует, что. Так как величина входит только в послед-нее из уравнений (5), то второе и третье уравнения отделяются:

(10)

где .

2.2 Другие модели

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

3. Идентификация параметров модели Лотки

3.1 Дифференциальные уравнения

Задачу Коши для уравнений Лотки (5) п. 2 запишем, используя более стан-дартные математические обозначения:

(1)

(2)

Задача Коши (17), (18) п. 1 будет следующей:

(3)

, (4)

Как видим, задача Коши (1), (2), (3), (4) полиномиальная, и для ее численного интегрирования можно применять метод рядов Тейлора.

3.2 Постановки задачи идентификации и функционалы МНК

Для конкретных биологических или иных моделей проводят реальные эксперименты по определению величин, от которых зависят функционалы типа (20) п. 1.3. Каждый реальный эксперимент имеет и свои возможности (часто весьма ограниченные) и свою цену (возможно высокую) определения каждой величины .

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

В эксперименте первого типа, при одном и том же начальном данном измеряются значения

(5)

одной из переменных в различные моменты, .

В эксперименте второго типа, при начальных данных, из-меряются значения

(6)

величин, в один и тот же момент времени .

В эксперименте третьего типа, при начальных данных, из-меряются значения

(7)

величин, в моменты времени, ,.

Соответствующие функционалы равны:

(8)

(9)

(10)

где — фиксированные весовые коэффициенты.

Градиентные уравнения и соответствующие начальные условия для этих функционалов следующие:

(11)

(12)

(13)

(14)

3.3 Как ускорить вычисления

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

Для того, чтобы ускорить вычисления на втором этапе, необходимо ускорить численное интегрирование исходных уравнений, уравнений в вариациях и градиентных уравнений. Исходные уравнения и уравнения в вариациях, как правило, полиномиальные и для их численного интегрирования можно использовать метод рядов Тейлора.

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

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

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

Перейдем к формулам. Уравнения точной градиентной задачи Коши

(15)

, (16)

где, мы хотим заменить на приближенные градиентные уравнения:

, (17)

где — полином по, а — набор его коэффициентов.

При этом мы хотим, чтобы величины

(18)

были достаточно малыми при

(19)

где — некоторое фиксированное число. Коэффициенты поли-нома можно найти методом наименьших квадратов с функционалом:

(20)

где, , а — весовые коэффициенты.

Отметим, что при малых в качестве можно рассмотреть полином степени 3 или 4, а при больших и/или — полином степени 2.

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

Мы опишем здесь постановку и результаты одного из численных экспери-ментов, проведенных в полном соответствии с рассмотренной выше схемой градиентного метода. Эти результаты опубликованы в работе.

Обратимся к дифференциальным уравнениям для модели Лотки в п. 3.1 и в численном эксперименте будем действовать по следующей схеме:

Фиксируем начальные данные

, , (21)

и параметры

, (22)

При этих значениях начальных данных и параметров

численным интегрированием задачи Коши (1),(2) находим значения концентрации реактанта в моменты времени, , то есть находим при .

Теперь можно имитировать «измерения» величин по формуле

, (23)

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

Фиксируем начальное приближение:

(24)

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

Об эффективности метода можно судить по затраченному процессорному времени и по величине относительной погрешности:

(25)

Результаты этого численного эксперимента приведены на рисунках 1, 2.

4. О других методах идентификации

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

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

Для минимизации функционала используется программа VARPRO Стэнфордского университета, а численное интегрирование исходных уравнений (для вычисления функционала) проводится при помощи интегратора SDRIV1 Дэвида Кахане.

1. В. Вольтерра, «Математическая теория борьбы за существование». Москва. «Наука», 1976.

2. Э. Хайрер, С. Нёрсетт, Г. Ваннер, «Решение обыкновенных дифференциальных уравнений», I. Нежесткие задачи. Москва. «Мир», 1990.

3. Э. Хайрер, Г. Ваннер, «Решение обыкновенных дифференциальных уравнений», II. Жесткие и дифференциально — алгебраические задачи. Москва. «Мир», 1999.

4. L.K. Babadzanjanz, J.A. Boyle, D.R. Sarkissian, and J. Zhu, «Parameter Identification for Oscillating Chemical Reactions Modelled by Systems of ODE», Journal of Computational Methods for Sciences and Engineering, 2002.

5. Bert W. Rust, ACMD, Robert W. Ashton, Chemical Science and Technology Laboratory, «Parameter Identifications», 7/15/2001: http://math.nist.gov/mcsd/Reports/95/yearly/node28.html

6. R. Haberman, «Mathematical Models. Mechanical Vibrations, Population Dynamics, and Traffic Flow. Classics in Applied Mathematics, 21», SIAM, Philadelphia, 1977.

7. A.J. Lotka, «Undamped oscillations derived from the law of mass action», Jour. Amer. Chem. Soc. 42 (1920), 1595−1599.

8. A.J. Lotka, «Elements of Physical Biology», Williams and Wilkins, Baltimore, 1925.