Ідентифікація параметрів осцилюючих процесів у живій природі модельованих диференціальними

Санкт-Петербурзький Державний Університет

Реферат

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

Виконала студентка 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) отримуємо:

(2)

де - Градієнт функції , А це дає:

(3)

(4)

(5)

Таким чином, вектор є напрямом найшвидшого ріс-та функції в точці , А вектор - Це напрямок найшвидшого її убування в цій точці.

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

Це означає, що задовольняє диференціальному рівнянню:

(6)

або в координатній формі:

(7)

До рівнянь (6) або (7) додаємо початкові умови:

(8)

або в координатній формі:

(9)

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

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

Якщо рішення задачі (6), (8) існує при , То при будь-якому та-ком отримуємо, що:

при (11)

при (12)

і ми маємо право очікувати, що

(13)

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

2. Рівняння у варіаціях

Розглянемо задачу Коші:

(14)

(15)

де - Параметри. Надалі ми розглянемо функціонали, які залежать від параметрів через рішення задачі Коші (14), (15). Тоді градієнтні рівняння будуть залежати від похідних по рішення задачі (14), (15), і ми повинні вміти їх обчислювати. Диференціюючи рівняння (14), (15) за отримуємо, що функції

(16)

задовольняють наступної задачі Коші:

(17)

(18)

Рівняння (17) щодо похідних (16) називають рівняннями в варіаціях для рівнянь (14).

3. Функціонали методу найменших квадратів

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

, (19)

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

У методі найменших квадратів знаходження (ідентифікації) параметрів розглядають функціонал

(20)

де - Фіксовані вагові коефіцієнти, а - Значення перших компонент розв'язання задачі (14), (15) в точці при заданих

У методі найменших квадратів вважають, що значення , Що доставляє мінімум цієї функції , Є адекватним наближенням до реального значення параметра для прийнятої моделі процесу.

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

(21)

Ці градієнтні рівняння треба доповнити початковими умовами:

(22)

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. LK Babadzanjanz, JA Boyle, DR 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. AJ Lotka, "Undamped oscillations derived from the law of mass action", Jour. Amer. Chem. Soc. 42 (1920), 1595-1599.

8. AJ Lotka, "Elements of Physical Biology", Williams and Wilkins, Baltimore, 1925.