Санкт-Петербургский Государственный Университет
Реферат
Идентификация параметров осциллирующих процессов в живой природе,моделируемых дифференциальными уравнениями
Выполниластудентка 312гр.
Варламова А.А.
Проверил ТокинИ.Б
Санкт-Петербург
2007
Оглавление
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 Осциллирующиехимические реакции
В некоторых химическихреакциях концентрации реагентов осциллируют в следующем смысле. Соединениекаких-то начальных веществ приводит к их химическому взаимодействию, врезультате чего образуются новые вещества, которые также начинаютвзаимодействовать с другими реагента-ми. В течении всех этих реакцийконцентрации реагентов колеблются и, на-конец, все химические преобразованиязавершаются и в качестве результата остаются какие-то определенные вещества,которые уже не реагируют между собой. Первая математическая модельосциллирующих химических реакций была предложена в работе Лотки [7].
Рассматриваетсяматематическая модель взаимодействия на молекулярном уровне веществ /> на основеследующих предположений:
1. При взаимодействии смолекулой вещества /> молекула вещества /> превращается в молекулувещества />.Это описывают в форме молекулярной ре-акции:
/> (1)
Такую реакцию относят кклассу автокаталитических, так как наличие вещества /> обеспечивает превращение другоговещества в />.
2. При взаимодействии смолекулой вещества /> молекула вещества /> пре-вращается вмолекулу вещества />, то есть происходитавтокаталитическая молекулярная реакция:
/> (2)
3. Вещество /> в то же времянеобратимо распадается, превращаясь в вещество />, то есть происходит молекулярнаяреакция
/> (3)
4. Скорости протеканияреакций (1), (2), (3) пропорциональны концентрациям веществ в левых частях этихреакций, то есть равны соответственно:
/>, />, />, (4)
где символами />, />, /> обозначеныконцентрации веществ />, />, /> со-ответственно, а коэффициенты /> -положительные числа.
5. Скорость измененияконцентрации каждого вещества равна сумме скоростей изменения концентрацийэтого вещества во всех реакциях, в которых оно участвует.
Из условий 1-5 следуютравенства:
/>,
/>,
/>,
/>, (5)
где /> - концентрация вещества />. Это системаОДУ Лотки.
2.1.2 Осцилляция популяцийв системе «хищник-жертва»
Первая экологическая модельтипа «хищник – жертва» была предложена в книге Лотки [8]. Она основана на техже уравнениях (5).
Пусть на острове живут жертвы/> (зайцы) ихищники /> (волки).Рассматривается математическая модель изменения величин /> (растительная пища длязайцев), />,/>, /> (умершиеволки) на основе следующих предположений:
1. Наличие зайцев /> и еды для них /> приводит кувеличению количества зайцев, что можно записать формулой:
/> (6)
2. Наличие волков /> и еды для них /> приводит кувеличению количества волков:
/> (7)
3. Волки умирают от болезнейили старости:
/> (8)
4. Скорость измененияколичества зайцев по формуле (6), скорость изменения количества волков поформуле (7) и скорость увеличения количеств умерших волков по формуле (8) равнысоответственно:
/>, />, />, (9)
где символами />, />, /> обозначены количестварастительной пищи, зайцев и волков, а /> - положительные коэффициенты.
5. Скорость изменения каждогоиз количеств /> (количество умерших волков) равнасумме скоростей изменения этих количеств в каждом из процессов (6), (7), (8), вкотором соответствующая величина /> участвует.
Из условий 1-5 следуютуравнения Лотки (5), только символы имеют другой смысл.
Более общие модели поведения /> хищников и /> жертв в различныхэко-логических ситуациях были предложены в лекциях Вольтерры [1]. В связи сэтим, уравнения Лотки (5) называют часто уравнениями Лотки-Вольтерра.
И все же большая часть работпо этой тематике посвящена даже более упрощенному по сравнению с моделью Лоткидвумерному случаю, так как это позволяет применять методы фазовой плоскости длядинамических систем.
Сведение модели (5) кдвумерной основано на предположении, что вели-чина /> постоянна. В случае моделиосциллирующих химических реакций это означает, что вещества /> достаточно много, а вслучае модели «хищник — жертва» это означает, что еды у зайцев достаточномного. Из этого предполо-жения следует, что />. Так как величина /> входит только в послед-нееиз уравнений (5), то второе и третье уравнения отделяются:
/>,
/>, (10)
где />.
2.2 Другие модели
Они излагаются вмногочисленных статьях и книгах. Кроме уже предложенных ранее, дадим здесьссылку еще на одну книгу [6].
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 Численный эксперимент
Мы опишем здесь постановку и результаты одногоиз численных экспери-ментов, проведенных в полном соответствии с рассмотреннойвыше схемой градиентного метода. Эти результаты опубликованы в работе [4].
Обратимся к дифференциальным уравнениям для модели Лотки в п. 3.1и в численном эксперименте будем действовать по следующей схеме:
1. Фиксируемначальные данные
/>, />, />, /> (21)
и параметры
/>, />, /> (22)
2. Приэтих значениях начальных данных /> и параметров
/> численным интегрированием задачиКоши (1),(2) находим значения концентрации реактанта /> в моменты времени />, />, то есть находим /> при />.Теперьможно имитировать «измерения» величин /> по формуле
/>, />, (23)
где /> - независимые случайные величины,равномерно распределенные меж-ду /> и />. Считаем, что /> - измерения, полученныев некотором реаль-ном эксперименте.
3. Фиксируемначальное приближение:
/> (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 ChemicalReactions Modelled by Systems of ODE”, Journal of Computational Methods forSciences and Engineering, 2002.
5. Bert W. Rust, ACMD, Robert W. Ashton,Chemical Science and Technology Laboratory, “Parameter Identifications”,7/15/2001: math.nist.gov/mcsd/Reports/95/yearly/node28.html
6. R.Haberman, “MathematicalModels. Mechanical Vibrations, Population Dynamics, and Traffic Flow. Classicsin 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.