Подякувати Студентському архіву довільною сумою
Admin
26.02.2023 12:38
Дякуємо, що користуєтесь нашим архівом!
Інформація про навчальний заклад
ВУЗ:
Національний університет Львівська політехніка
Інститут:
Не вказано
Факультет:
Прикладна математика
Кафедра:
Не вказано
Інформація про роботу
Рік:
2008
Тип роботи:
Методичні вказівки до лабораторної роботи
Предмет:
Чисельні методи
Група:
МЕ
Частина тексту файла (без зображень, графіків і формул):
МІНІСТЕРСТВО ОСВІТИ І НАУКИ УКРАЇНИ
НАЦІОНАЛЬНИЙ УНІВЕРСИТЕТ "ЛЬВІВСЬКА ПОЛІТЕХНІКА"
ЛІНІЙНІ БАГАТОКРОКОВІ МЕТОДИ ЧИСЕЛЬНОГО РОЗВ’ЯЗУВАННЯ
ЗАДАЧІ КОШІ ДЛЯ ЗВИЧАЙНИХ ДИФЕРЕНЦІАЛЬНИХ РІВНЯНЬ
Методичні вказівки
до лабораторної роботи з курсу
"Чисельні методи"
для студентів базового напрямку
6.0802 "Прикладна математика"
Затверджено
на засіданні кафедри
прикладної математики
Протокол №4 від 9.11.2006
Львів – 2008
Лінійні багатокрокові методи чисельного розв’язування задачі Коші для звичайних диференціальних рівнянь: Методичні вказівки до лабораторної роботи з курсу "Чисельні методи" для студентів базового напрямку 6.0802 "Прикладна математика" / Укл.: М.В. Кутнів, Я.В. Пізюр – Львів: Видавництво Національного університету "Львівська політехніка", 2007. – 14 с.
Укладачі Кутнів М.В., канд.фіз.-мат.наук, доц.
Пізюр Я.В., канд.фіз.-мат.наук, доц.
Відповідальний за випуск Костробій П.П., канд.фіз.-мат.наук, проф.
Рецензент Каленюк П.І., д-р.фіз.-мат.наук, проф.
ТЕОРЕТИЧНА ЧАСТИНА
§1. Лінійні багатокрокові методи
1. Методи Адамса. Введемо на інтервалі рівномірну сітку з кроком . Якщо рівняння (1) §2 проінтегрувати на відрізку , то одержимо
. (1)
Припустимо, що нам відомі наближені значення точного розв’зку задачі (1), (2) §2, тоді можна вважати також, що ми маємо і величини . Замінимо функцію в (1) інтерполяційним многочленом Ньютона, який проходить через точки . Його можна виразити через різниці назад:
(2)
Тоді чисельний аналог (1) буде задаватися формулою
.
Після заміни змінної в останньому інтегралі та підставляння виразу (2) будемо мати
(3)
де коефіцієнти обчислюються за формулами
Формула (3) дозволяє визначити явно, тому її називають явним методом Адамса.
Розглянемо частинні випадки (3). Якщо для , обчислити та виразити різниці назад через , то одержимо такі формули
Зауважимо, що при ми маємо явний метод Ейлера.
Формули (3) одержані при інтегруванні інтерполяційного многочлена від до , тобто зовні інтервалу інтерполяції Добре відомо, що зовні цього інтервалу інтерполяційний многочлен дає досить погане наближення. Тому дослідимо також методи, що грунтуються на інтерполяційному многочлені, який додатково використовує точку , тобто
(4)
Підставляючи цей многочлен у (1), одержимо наступний неявний метод:
(5)
де коефіцієнти визначаються за формулами
Наведемо приклади формул (5). При будемо мати неявний метод Ейлера
при правило трапецій
Насправді ці два методи – однокрокові. При одержимо відповідно такі методи
Формули (5) визначають неявно (на кожному кроці для обчислення необхідно розв’язати нелінійне рівняння), а тому вони називаються неявними методами Адамса.
Неявні формули Адамса мають загальний вигляд:
. (6)
2. Формули диференціювання назад. Багатокрокові формули Адамса основані на чисельному інтегруванні, тобто інтеграл в (1) апроксимується деякою квадратурною формулою. Тепер розглянемо багатокрокові методи, які грунтуються на чисельному диференціюванні.
Припустимо, що відомі значення розв’язку диференціального рівняння (1) §2 . Щоб вивести формулу для , використаємо інтерполяційний многочлен , який проходить через точки Як і многочлен (4) , його можна виразити через різниці назад , а саме
(7)
Визначимо тепер невідоме значення так, щоб многочлен задовольняв диференціальне рівняння хоча б в одному вузлі сітки, тобто
Враховуючи , що , продиференціюємо (7) по змінній
Для одержимо явні формули
,
де
.
При будемо мати явний метод Ейлера, а при правило середньої точки
У випадку формула має вигляд
.
Однак вона нестійка, як і всі решта формул при (див. п.4), а тому непридатна для розрахунків.
Кращі властивості мають формули, які одержуються з (7) при . Це неявні формули
(8)
з коефіцієнтами
які після диференціювання набувають вигляду
при .
Тому (8) зводиться до формули
.
Такі багатокрокові методи називають формулами диференціювання назад. Вони вперше були виведені Кертісом і Хіршфельдером.
Наведемо приклади цих формул, виразивши різниці назад через
Формули диференціювання назад мають вигляд:
.
3. Умови порядку апроксимації. Методи Адамса та формули диференціювання назад є частинними випадками лінійних багатокрокових (різницевих) методів, які описуються загальною формулою
(9)
де дійсні числа. Будемо вважати, що виконані такі умови :
Рівняння (9) треба розглядати як рекурентне співвідношення, яке виражає нове значення через знайдені раніше значення Для початку розрахунку необхідно задати значень . Значення визначається вихідною задачею, а саме Величини можна обчислити, наприклад, за допомогою методів Рунге-Кутта або багатокрокових методів нижчого порядку. Надалі будемо вважати, що початкові значення задані.
Якщо , то формула (9) задає явний метод; якщо , то – неявний.
Зауважимо, коефіцієнти в (9) визначені з точністю до множника. Щоб ліквідувати цю довільність, будемо вважати, що виконується умова
(10)
яка означає, що права частина різницевого рівняння (9) апроксимує праву частину рівняння (1) §2 .
Локальною похибкою багатокрокового методу (9) називається величина при точних значеннях При це означення збігається з означенням локальної похибки однокрокових методів.
Похибкою апроксимації або нев’язкою різницевого методу (9) називається функція
(11)
яка одержується в результаті підставляння точного розв’язку диференціальної задачі (2.1), (2.2) в різницеве рівняння (9) та ділення на .
Лема. Розглянемо диференціальне рівняння (2.1) з неперервно диференційовною функцією і розв’язком . Тоді для локальної похибки лінійного багатокрокового методу виконується рівність
,
де – деяке проміжне значення між і .
Доведення. З означення локальної похибки
За теоремою про скінченні прирости будемо мати
Звідси випливає твердження леми.
Кажуть, що багатокроковий метод має -й порядок апроксимації, якщо виконується одна з умов: або . Зауважимо, що згідно з лемою ці умови еквівалентні.
Дослідимо, як впливає вибір коефіцієнтів на похибку апроксимації багатокрокового методу (9). Розкладаючи функції в точці за формулою Тейлора, одержимо
Підставляючи ці розклади у вираз для похибки апроксимації (11), будемо мати
Звідси випливає, що похибка апроксимації має -й порядок, якщо виконуються умови
. (12)
Разом з умовою нормування (10) ми отримали систему лінійних алгебраїчних рівнянь відносно невідомих
Для того, щоб ця система не була перевизначена, необхідно вимагати, щоб Ця вимога означає, що порядок апроксимації лінійних -крокових методів не може перевищувати . Отже, найвищий досяжний порядок апроксимації неявних -крокових методів дорівнює , а явних – .
Явні методи Адамса дають точний розв’язок тоді, коли і є многочленом, степінь якого менший . Тоді за допомогою явних методів Адамса одержуємо точні розв’язки диференціальних рівнянь
,
які мають вигляд . Це означає, що похибка апроксимації дорівнює нулю, а отже
.
Оскільки ця рівність збігається з рівністю (12) при , то явні методи Адамса мають порядок не нижче .
Неявні методи Адамса дають точний розв’язок, якщо є многочленом на одиницю вищого степеня, ніж для явних, тому вони мають порядок не нижче .
Аналогічно можна показати, що формули диференціювання назад мають порядок не нижче .
Приклад 1. Знайти головний член похибки апроксимації двокрокового методу
.
Показати, що метод має третій порядок апроксимації.
Розв’язування. Враховуючи співвідношення
,
,
похибку апроксимації цього методу розкладемо у ряд Тейлора в околі точки
.
Отже, головний член похибки апроксимації має вигляд , а метод має третій порядок апроксимації.
4. Стійкість багатокрокових методів. Високого порядку апроксимації та малої локальної похибки ще недостатньо для того, щоб багатокроковий метод був придатний для практичних розрахунків. Чисельний метод може бути чутливим до малих збурень початкових даних тобто бути нестійким. Це означає, що як завгодно малі зміни початкових даних можуть викликати як завгодно великі зміни розв’язку навіть при і фіксованому .
Розглянемо явний лінійний двокроковий метод третього порядку апроксимації вигляду
Застосуємо його до задачі
тоді будемо мати
(13)
За початкові візьмемо значення точного розв’зку
За аналогією з лінійними диференціальними рівняннями будемо шукати розв’зок рівняння (13) у вигляді де – деяка невідома константа. У результаті одержимо характеристичне рівняння
(14)
Загальний розв’язок різницевого рівняння (13) обчислюється за формулою
(15)
де – корені рівняння (14), а коефіцієнти визначаються з початкових умов. Оскільки апроксимує , перший член в (15) апроксимує значення точного розв’язку в точці . Збурення в розв’язок різницевого рівняння вносить другий член, який називають "паразитним", тому що при . Цей "паразитний" розв’язок з ростом стає великим і починає переважати у розв’язку .
Звернемося тепер до питання стійкості багатокрокового методу (9). Важливу роль відіграє поведінка розв’язку при . Очевидно, що (9) при зводиться до формули
. (16)
Її можна розглядати як чисельний метод (9), застосований до розв’язання диференціального рівняння
.
Підставляючи в (16) , одержимо характеристичне рівняння
. (17)
Багатокроковий метод (9) називається стійким (нуль-стійким), якщо корені характеристичного рівняння (17) задовольняють кореневу умову, тобто всі корені за модулем не більші за одиницю, а серед коренів модуль яких дорівнює одиниці немає кратних.
Для явного і неявного методів Адамса характеристичне рівняння (17) має вигляд . Крім простого кореня, рівного , характеристичне рівняння має нульовий корінь кратності . Отже, методи Адамса – стійкі.
Можна показати, що -крокова формула диференціювання назад стійка при і нестійка при .
Оскільки методи Рунге-Кутта є однокроковими, то для них характеристичне рівняння має вигляд . Це рівняння має лише один корінь , а отже, методи Рунге-Кутта – стійкі.
Справджується твердження.
Теорема. Якщо багатокроковий метод (9) стійкий і має -й порядок апроксимації, то він збіжний з -м порядком точності.
Доведення див., наприклад, [23].
5. Методи Нордсіка. Усі лінійні багатокрокові методи чисельного інтегрування еквівалентні знаходженню многочлена, який апроксимує розв’язок . Ідея методів Нордсіка полягає у тому, щоб зобразити цей многочлен через похідні від нульового до -го порядку включно, тобто з допомогою вектора
величини мають зміст наближених зачень .
Щоб визначити процедуру чисельного інтегрування, необхідно задати правило знаходження за відомим і диференціальним рівнянням (2.1). При використанні розкладу в ряд Тейлора (наприклад, при ) таке правило має вигляд
(18)
У загальному випадку прогнозоване значення
(19)
де –трикутна матриця Паскаля, елементи якої визначаються формулою
Оскільки ця апроксимація не використовує рівняння (2.1), то рівність , взагалі кажучи, не буде виконуватися, тому коректуємо , додаючи до нього величину, кратну , і будемо вимагати, щоб виконувалася рівність
(20)
де . Друга компонента вектора повинна бути рівна 1 , щоб виконувалася умова
Співвідношення (19) і (20) можуть бути об’єднані в одну формулу
(21)
де . Процедуру (21) називають процедурою Нордсіка.
Дослідимо метод Нордсіка на стійкість. Для цього підставимо в (21) , тоді одержимо
Наприклад, при ця матриця має вигляд
.
Для стійкості методу необхідно, щоб всі власні числа матриці задовольняли кореневу умову. Зауважимо, що і є власними значеннями матриці , а її характеристичний многочлен не залежить від . Величини зручно вибрати таким чином, щоб всі решта власні значення були рівні нулю. У випадку це виконується при . Коефіцієнт можна вибрати з умови перетворення в нуль константи похибки методу. У нашому випадку , а весь метод задається вектором
Зазначимо, що цей метод еквівалентний трикроковому методу Адамса. Дійсно, запишемо розклад (18) для точки
тоді
Розв’яжемо цю систему рівнянь відносно , тоді будемо мати
(22)
З рівності (21) для першої компоненти вектора Нордсіка отримаємо
Підставивши (22) в останню рівність, одержимо трикроковий неявний метод Адамса.
Можна показати, що метод Нордсіка еквівалентний деякій багатокроковій формулі порядку не нижче . З допомогою вектора Нордсіка зручно змінювати крок чисельного інтегрування (див. §7).
§2. Чисельне інтегрування жорстких систем
звичайних диференціальних рівнянь
Чисельні методи розв’язання задачі Коші для звичайних диференціальних рівнянь покоординатно без будь-яких змін переносяться на системи звичайних диференціальних рівнянь.
Наприклад, лінійні багатокрокові методи у випадку систем рівнянь можна записати у вигляді
, (1)
де – вектори наближень до розв’язку і правої частини системи диференціальних рівнянь у відповідних точках. Однак при чисельному розв’язанні систем звичайних диференціальних рівнянь можуть виникнути додаткові труднощі, пов’язані з жорсткістю системи.
1. Поняття жорсткої задачі. Суть явища жорсткості полягає в тому, що розв’язок, який потрібно обчислити, змінюється повільно, однак існують швидкі згасаючі збурення. Наявність таких збурень перешкоджає знаходженню чисельного розв’язку, який повільно змінюється.
Розглянемо приклад жорсткої задачі
точний розв’язок якої має вигляд
Цей розв’язок складається з двох компонент, одна з яких змінюється значно швидше (оскільки , наприклад ), ніж друга . Після дуже малого відрізка часу , де , швидка (жорстка) компонента прямує до нуля і суттєвий вклад у розв’язок робить повільна (гладка) компонента.
Для виявлення деяких особливостей чисельних методів застосуємо до цієї задачі явний та неявний методи Ейлера:
Тоді для явного методу Ейлера будемо мати
а для неявного
Одержані чисельні розв’язки порівняємо з точним
Вираз можна інтерпретувати як збурення повільного розв’язку в момент . Це збурення швидко спадає завдяки , де , тобто при . Чисельний метод повинен бути в змозі подавляти різницю для значень . Різниця у випадку неявного методу Ейлера швидко спадає і для будь-якого фіксованого при . Для явного методу різниця буде зменшуватися тільки у випадку , тобто при . Ця умова накладає сильне обмеження на величину кроку при . З іншого боку, якщо , то апроксимується виразом досить точно для значно більших , ніж ті, що задовольняють умову . Отже, величина кроку обмежується чисельною стійкістю, а не точністю.
Для лінійних систем диференціальних рівнянь зі сталими коефіцієнтами подібні компоненти розв’язку, які сильно відрізняються, виникають, коли матриця системи містить сильно розкидані власні значення. Розглянемо, наприклад, систему
з початковими умовами . Власні значення матриці коефіцієнтів цієї системи
цієї системи дорівнюють: . Тоді розв’язок буде мати вигляд
де – швидка компонента розв’язку, а – повільна.
Розглянемо нелінійну систему звичайних диференціальних рівнянь (1.2). Нехай i – два різні розв’язки системи. Утворимо різницю , яка задовольняє рівняння
. (2)
Проведемо розклад правої частини (2) за формулою Тейлора. Оскільки
то маємо
У результаті розкладу система (2) набуває вигляду
(3)
де через позначено матрицю Якобі з елементами
Знехтувавши в (3) величиною , одержимо так звану систему рівнянь першого наближення
,
яка є системою лінійних диференціальних рівнянь відносно , оскільки функція задана. Припустимо, що в околі точки матриця змінюється достатньо мало, тоді одержимо систему диференціальних рівнянь зі сталими коефіцієнтами
Припустимо, що – власні числа матриці , які в загальному випадку є комплексними числами. Задача (1.1) називається жорсткою в околі точки , якщо
Нехай – власне число з найбільшою за модулем негативною дійсною частиною, тоді величина буде найменшою часовою константою. Мірою жорсткості задачі є величина
(4)
яка називається коефіцієнтом жорсткості. Якщо , то система буде сильно жорсткою, якщо , тоді система нежорстка.
На практиці жорсткі задачі розпізнають, виходячи з фізичних міркувань. Будь-яка фізична система, яка моделюється звичайними диференціальними рівняннямии і має фізичні компоненти з часовими константами, які сильно відрізняються, зводиться до жорсткої задачі. Фізичні компоненти з найменшими часовими сталими дуже швидко змінюються і роблять задачу жорсткою, тоді як дослідник проявляє найбільший інтерес до компонент з великими часовими константами. У жорсткій задачі розв’язки, які змінюються повільно, визначаються саме цими останніми компонентами. Існуючі методи чисельного розв’язання жорстких задач дозволяють інтегрування з розмірами кроків такого ж порядку, як і великі часові константи. Класичні явні методи можуть застосовуватися тільки з розмірами кроків порядку найменших часових констант і тому є надзвичайно неефективними.
Однією з властивостей жорстких задач є наявність великої константи Ліпшиця
(5)
Якщо вважати, що часовий маштаб задачі зв’язаний з розв’язком, який змінюється повільно на інтервалі , то вираз (5) можна замінити на більш конкретну умову .
2. Абсолютна стійкість чисельних методів. Для дослідження стійкості чисельних методів інтегрування жорстких систем диференціальних рівнянь, розглянемо лінійну систему зі сталими коефіцієнтами
. (6)
Застосовуючи лінійний багатокроковий метод (1) до системи (6), одержимо
. (7)
Припустимо, що матриця має різних власних значень , тоді існує невироджена матриця така, що , де – діагональна матриця з елементами . Якщо покласти і помножити систему (6) зліва на , то одержимо
або
(8)
Зробимо аналогічну заміну в (7), тоді
. (9)
Співвдношення (9) збігається зі скінченно-різницевою апроксимацією лінійним багатокроковим методом (1) рівнянь (8). Аналогічні висновки можна зробити відносно методів Рунге - Кутта. Отже, стійкість чисельних методів будемо досліджувати на скалярному модельному рівнянні
(10)
де – комплексне число, оскільки власні значення матриці , взагалі кажучи, комплексні числа.
Багатокроковий метод, застосований до модельного рівняння (10), має вигляд
де – комплексний параметр. Якщо розв’язок цього рівняння шукати у вигляді , то одержимо характеристичне рівняння
. (11)
Зауважимо, що якщо однокроковий метод застосовувати до модельного рівняння (10), то , де називається функцією стійкості однокрокового методу. Характеристичне рівняння для однокрокового методу має вигляд .
Чисельний метод називається абсолютно стійким (стійким при будь-яких ) для , якщо для цього корені характеристичного рівняння (11) задовольняють кореневу умову, тобто всі корені за модулем не більші за одиницю, а серед коренів модуль яких рівний одиниці немає кратних. Областю абсолютної стійкості називається множина всіх точок комплексної площини , для яких даний метод є абсолютно стійким.
Однокроковий метод буде абсолютно стійким, якщо .
Метод називається -стійким, якщо його область абсолютної стійкості включає всю півплощину . Суть наведеного означення полягає в тому, що -стійкий чисельний метод є абсолютно стійким, якщо стійким є розв’язок вихідного диференціального рівняння.
Явний метод Ейлера, застосований до модельного рівняння (10), має вигляд
Вимога абсолютної стійкості приводить до оцінки
або
де . Область абсолютної стійкості в цьому випадку буде всередині круга одиничного радіуса з центром у точці .
Для неявного методу Ейлера
область абсолютної стійкості – зовні круга одиничного рідіуса з центром в точці .
Отже, неявний метод Ейлера є -стійкий, а явний – не -стійкий.
Розглянемо формулу трапецій
.
Для рівняння (10) цей метод має вигляд
Звідси випливає, що тоді і тільки тоді, коли . Отже, формула трапецій теж -стійка.
Для знаходження області стійкості лінійного багатокрокового методу виразимо з рівняння (11) параметр через змінну , тоді дістанемо
. (12)
Покладемо у рівнянні (12) таким, що приймає всі значення на одиничному колі і обчислимо відповідні значення . Оскільки будь-яке на одиничному колі можна записати у вигляді , то можна підставити замість у рівняння (12), що приведе до
. (13)
Множина точок, породжених (13), є геометричним місцем точок , для яких . Оскільки коло є замкненою кривою, то також замкнена. Геометричне місце точок є симетричним відносно осі , тому достатньо обчислити тільки для .
При розв’язуванні жорстких систем диференціальних рівнянь було б бажано використовувати саме -стійкі методи, тому що умови їх стійкості не накладають обмежень на крок . Однак серед явних методів (Рунге - Кутта, багатокрокових) не існує -стійких.
Дальквіст довів, що не існує -стійкого неявного лінійного багатокрокового методу порядку апроксимації вище другого.
У зв’язку з цим було введено ще одне означення стійкості.
Чисельний метод називається -стійким, якщо область його стійкості містить кут (рис.1).
Зокрема, -стійкість збігається з -стійкістю.
Доведено, що ні для якого не існує явного -стійкого методу.
Найбільш поширеним класом методів розв’язування жорстких задач є формули диференціювання назад, області стійкості яких зображені на рис.2 (області стійкості розташовані зовні відповідних кривих ).
Формули диференціювання назад від 1-го до 6-го порядків апроксимації – -стійкі, причому чим більший порядок, тим менше .
Приклад 1. При яких значеннях параметра чисельний метод
буде -стійким?
Розв’язування. Застосуємо цей метод до модельного рівняння (10), тоді одержимо
де функція стійкості
є аналітичною функцією в лівій півплощині комплексної площини . Згідно з теоремою про максимум модуля аналітичної функції, максимум досягається на границі області, тобто на уявній осі . Отже, умова -стійкості буде виконуватися, якщо . Остання умова для заданої функції буде мати вигляд . Звідси . Отже, чисельний метод -стійкий за умови .
§3. Реалізація лінійних неявних багатокрокових методів
Для обчислення за формулою (6.1) необхідно розв’язати систему нелінійних алгебраїчних рівнянь з невідомими вигляду
(1)
де
У найпростішому випадку для розв’язування системи (1) можна використати метод простої ітерації
(2)
де – задане початкове наближення до , для обчислення якого можна використати явні багатокрокові методи, –номер ітерації. Цей ітераційний процес ефективний тільки у випадку нежорстких задач. Для жорстких задач він збіжний при менших, ніж найменша часова константа системи диференціальних рівнянь.
У випадку жорсткої задачі систему (1) розв’язують з допомогою методу Ньютона. Запишемо (1) у вигляді
.
Тоді ітераційний метод Ньютона буде мати вигляд
або
де – одинична матриця порядку – матриця Якобі правої частини системи, . Обчислення за цією формулою вимагає значних затрат, пов’язаних з необхідністю обчислення матриці Якобі, і розв’язанням системи лінійних алгебраїчних рівнянь з невідомими для кожного .
Ці затрати можуть бути зменшені за рахунок застосування модифікованого методу Ньютона і за рахунок врахування структури матриці Якобі. Для обчислення з допомогою модифікованого методу Ньютона необхідно розв’язати систему лінійних алгебраїчних рівнянь
Щоб розв’язати таких систем лінійних рівнянь методом Гаусса необхідно один раз розкласти матрицю
у вигляді ( – ліва трикутна матриця, – права трикутна матриця), що відповідає прямому ходу методу Гаусса і розв’язати систем з трикутними матрицями, що відповідає оберненому ходу методу Гаусса. Матрицю Якобі можна обчислювати як чисельно, так і аналітично. Величину часто вибирають рівною 3.
Ефективна програма чисельного інтегрування повинна передбачати зміну величини кроку. Однак у випадку багатокрових методів процедура зміни кроку викликає ускладнення, оскільки ці методи вимагають значень чисельного розв’язку в точках, які розташовані на однаковій віддалі. Існують дві можливості виходу з цього:
1) визначення з допомогою інтерполяції многочленами значень у точках, які необхідні при зміні кроку;
2) побудова багатокрокових методів для нерівномірної сітки (див., напр., [23])
Розглянемо детальніше перший підхід. Лінійні багатокрокові методи можна використовувати у формі з розділеними різницями (див. §5 п.1,2) або з допомогою зображення Нордсіка , кожна компонента якого апроксимує відповідну похідну точного розв’язку (див. §5 п.5). При такому способі зберігання інформації початкове наближення обчислюють за формулою де –трикутна матриця Паскаля, елементи якої визначаються за формулою
Повний алгоритм прогнозу і корекції при можна записати у вигляді
де матриця-рядок визначає метод (9), . Довільна зміна кроку у векторі Нордсіка досягається просто множенням -го коефіцієнта на , причому , де – новий крок , тобто змінюється так: . Однак, у цьому випадку при зміні порядку необхідно обчислити новий вектор Нордсіка. У формі з розділеними різницями зручніше змінювати порядок методу, але важче змінювати крок.
При застосуванні лінійних багатокрокових методів важливим є питання контролю похибки, вибору кроку чисельного інтегрування, порядку методу, який є найбільш точним і ефективним у даній точці .
Головний член локальної похибки лінійного багатокрокового методу порядку апроксимації має вигляд
де – константа, яка залежить від порядку .
Припустимо, що до точки чисельне інтегрування проходило успішно, а для знаходження ми вибрали крок і порядок . Щоб вирішити, чи придатне , необхідно оцінити локальну похибку. Для практичної оцінки локальної похибки методу порядку можна використати, наприклад, формулу
, (3)
де – розв’язок, одержаний методом порядку . На основі формули (3) обчислюється похибка
,
де –маштабний множник. Для обчислення абсолютної похибки , а для відносної . Можна використати змішане маштабування (див. §4).
Далі необхідно вибрати новий крок і порядок. Ідея вибору кроку полягає в знаходженні найбільшого , при якому похибка не перевищує . Похибку в точці можна наблизити її значенням в точці , тоді крок вибирається з умови
.
Звідси
.
Коефіцієнт дозволяє врахувати вплив величини . Для вибору оптимального порядку вибирають кроки
,
де
.
Тоді новий крок , а порядок , де .
Необхідно зазначити, що в кожній з реалізацій багатокрокових методів описаний алгоритм, удосконалений і доповнений ще рядом процедур. Наприклад, крок зберігається сталим, якщо відношення близьке до , що спрощує обчислення.
§4.Коротка характеристика програми
Систему ЗДР будемо розв’язувати за допомогою програми STIFF [4]. Програма контролює точність чисельного розв’язку і автоматично вибирає крок і порядок методу з метою розв’язати задачу з мінімальними затратами. Написана на мові FORTRAN з подвійною точністю. В програму STIFF включені неявні методи Адамса від 1-го по 12-й порядок точності, формули диференціювання назад від 1-го по 5-й з фіксованим кроком [].
Для використання програми необхідно викликати підпрограму STIFF. Зв’язок з викликаючою (головною) програмою здійснюється за допомогою таких формальних параметрів:
SUBROUTINE STIFF(Y,YMAX,ERROR,PW,FSAVE,IWORK,NYDIM)
IMPLICIT REAL*8 (A-H,O-Z)
LOGICAL EVALJA,CONVER
COMMON/STCOM1/ T,H,HMIN,HMAX,EPS,N,MF,KFLAG,JSTART,MAXDER
COMMON/STCOM2/ HUSED,NQUSED
COMMON/STCOM3/ ML,MU
COMMON/STCOM4/ NSTEP,NFUN,NJAC
DIMENSION Y(NYDIM,1),EL(13),TQ(4),YMAX(1)
DIMENSION ERROR(1),PW(1),FSAVE(1),IWORK(1)
C -------------------------------------------------------------
C ВХІДНІ ПАРАМЕТРИ:
С ------------------
C Y ДВОВИМІРНИЙ МАСИВ РОЗМІРУ NYDIM*LMAX
C (LMAX=13 ДЛЯ МЕТОДІВ АДАМСА, LMAX=6
C ДЛЯ ФОРМУЛ ДИФЕРЕНЦІЮВАННЯ НАЗАД).
C В Y(I,1), I=1,...,N ЗАДАЮТЬСЯ ПОЧАТКОВІ ЗНАЧЕННЯ.
C
C NYDIM ВИКОРИСТОВУЄТЬСЯ ДЛЯ ОПИСУ РОЗМІРІВ МАСИВІВ
C (NYDIM>=N).
С
C N КІЛЬКІСТЬ ДИФЕРЕНЦІАЛЬНИХ РІВНЯНЬ ПЕРШОГО ПОРЯДКУ.
С
C T ПОЧАТКОВЕ ЗНАЧЕННЯ НЕЗАЛЕЖНОЇ ЗМІННОЇ.
С
C H ВЕЛИЧИНА ПОЧАТКОВОГО КРОКУ.
С
C HMIN МІНІМАЛЬНЕ АБСОЛЮТНЕ ЗНАЧЕННЯ КРОКУ ІНТЕГРУВАННЯ.
C ОСКІЛЬКИ СПОЧАТКУ ЗАДАЧА РОЗВ'ЯЗУЄТЬСЯ МЕТОДОМ
C ПЕРШОГО ПОРЯДКУ ТОЧНОСТІ, ТО БАЖАНО HMIN
C ВИБИРАТИ НА РІВНІ ШУМУ ЕОМ.
C
C HMAX МАКСИМАЛЬНЕ АБСОЛЮТНЕ ЗНАЧЕННЯ КРОКУ ІНТЕГРУВАННЯ.
C
C EPS ТОЧНІСТЬ ЧИСЕЛЬНОГО РОЗВ'ЯЗКУ.
C
C YMAX МАСИВ РОЗМІРУ NYDIM.
С ЯКЩО ВИМАГАЄТЬСЯ КОНТРОЛЬ АБСОЛЮТНОЇ ПОХИБКИ, ТО
C ПЕРЕД ПЕРШИМ ВИКЛИКОМ STIFF НЕОБХІДНО ЗАДАТИ
C YMAX(I)=1.D0.
C ЯКЩО ВИМАГАЄТЬСЯ КОНТРОЛЬ ВІДНОСНОЇ ПОХИБКИ, ТО
C ПЕРЕД КОЖНИМ ВИКЛИКОМ STIFF НЕОБХІДНО ЗАДАТИ
C YMAX(I)=DMAX1(DABS(Y(I,1)),FLOOR),
С FLOOR - ДОДАТНЕ ЧИСЛО (ФІЗИЧНИЙ НОЛЬ ЗАДАЧІ),
С ВИКОРИСТОВУЄТЬСЯ ДЛЯ ТОГО, ЩОБ YMAX(I) НЕ
C ДОРІВНЮВАЛО 0.
С ЯКЩО ВИМАГАЄТЬСЯ КОНТРОЛЬ КОМБІНОВАНОЇ ПОХИБКИ:
C АБСОЛЮТНОЇ ПРИ DABS(Y(I,1))<=1 І ВІДНОСНОЇ (ПО
C ВІДНОШЕННЮ ДО НАЙБІЛЬШОГО ЗНАЧЕННЯ DABS(Y(I,1)) ПРИ
C DABS(Y(I,1))>1, ТО ПЕРЕД ПЕРШИМ
C ВИКЛИКОМ STIFF НЕОБХІДНО ЗАДАТИ
C YMAX(I)=DMAX1(DABS(Y(I,1)),1.D0),
C А ПРИ НАСТУПНИХ ЗВЕРТАННЯХ ДО STIFF НЕОБХІДНО
C ВИБРАТИ YMAX(I)=DMAX1(DABS(Y(I,1)),YMAX(I)).
C
C MF ВКАЗІВНИК МЕТОДУ, ДОДАТНЕ ЧИСЛО З ДВОМА ДЕСЯТКОВИМИ
C ЦИФРАМИ MF=METH*10+MITER;
C METH=1 - МЕТОД АДАМСА;
C METH=2 - ФОРМУЛИ ДИФЕРЕНЦІЮВАННЯ НАЗАД;
C MITER=0 - ВИКОРИСТОВУЄТЬСЯ МЕТОД ПРОСТОЇ ІТЕРАЦІЇ;
C MITER=1 - ВИКОРИСТОВУЄТЬСЯ МОДИФІКОВАНИЙ
C ІТЕРАЦІЙНИЙ МЕТОД НЬЮТОНА, МАТРИЦЯ ЯКОБІ
C ОБЧИСЛЮЄТЬСЯ АНАЛИТИЧНО ПІДПРОГРАМОЮ PEDERV;
C MITER=2 - ВИКОРИСТОВУЄТЬСЯ МОДИФІКОВАНИЙ
C ІТЕРАЦІЙНИЙ МЕТОД НЬЮТОНА, МАТРИЦЯ ЯКОБІ
C ОБЧИСЛЮЄТЬСЯ ЧИСЕЛЬНИМ ДИФЕРЕНЦІЮВАННЯМ;
C MITER=3 - ВИКОРИСТОВУЄТЬСЯ МОДИФІКОВАНИЙ
C ІТЕРАЦІЙНИЙ МЕТОД НЬЮТОНА, МАТРИЦЯ ЯКОБІ
C ЗАМІНЮЄТЬСЯ ДІАГОНАЛЬНОЮ МАТРИЦЕЮ;
C MITER=4 - ТЕ САМЕ, ЩО І MITER=1, АЛЕ МАТРИЦЯ ЯКОБІ
C - СТРІЧКОВА;
C MITER=5 - ТЕ САМЕ, ЩО І MITER=2, АЛЕ МАТРИЦЯ ЯКОБІ
C - СТРІЧКОВА.
C
C JSTART ІНДИКАТОР ТИПУ ВИКЛИКУ ПІДПРОГРАМИ.
C JSTTART=0 - ПРОВЕСТИ ПЕРШИЙ КРОК ІНТЕГРУВАННЯ;
C JSTTART>0 - ПРОВЕСТИ НАСТУПНИЙ КРОК ІНТЕГРУВАННЯ ЗА
C УМОВИ, ЩО ПОПЕРЕДНІЙ КРОК БУВ УСПІШНИМ
C І КОРИСТУВАЧ НЕ ЗМІНИВ ПАРАМЕТРИ;
C JSTTART<0 - ПОВТОРИТИ ОСТАННІЙ КРОК З НОВИМ
C ЗНАЧЕННЯМ H АБО EPS, АБО MF
C (АБО З НОВИМИ ЗНАЧЕННЯМИ ЦИХ ПАРАМЕТРІВ
C ОДНОЧАСНО).
C
C MAXORD МАКСИМАЛЬНИЙ ПОРЯДОК МЕТОДУ. ДЛЯ МЕТОДІВ АДАМСА
C MAXORD<=12, А ДЛЯ ФОРМУЛ ДИФЕРЕНЦІЮВАННЯ
C НАЗАД MAXORD<=5.
C
C MU КІЛЬКІСТЬ НЕНУЛЬОВИХ НАДДІАГОНАЛЕЙ СТРІЧКОВОЇ
C МАТРИЦІ ЯКОБІ MU=MAX(j-i, i<=j, J(i,j).NE.0).
C
C ML КІЛЬКІСТЬ НЕНУЛЬОВИХ ПІДДІАГОНАЛЕЙ СТРІЧКОВОЇ
C МАТРИЦІ ЯКОБІ ML=MAX(i-j, i>=j, J(i,j).NE.0).
C
C
C ВИХІДНІ ПАРАМЕТРИ:
С ------------------
C T ВИХІДНЕ ЗНАЧЕННЯ НЕЗАЛЕЖНОЇ ЗМІННОЇ.
C
C Y ДВОВИМІРНИЙ МАСИВ РОЗМІРУ NYDIM*LMAX
C (LMAX=13 ДЛЯ МЕТОДІВ АДАМСА, LMAX=6 ДЛЯ
C ФОРМУЛ ДИФЕРЕНЦІЮВАННЯ НАЗАД) ЗАЛЕЖНИХ ЗМІННИХ
C ТА ЇХ МАШТАБОВАНИХ ПОХІДНИХ ОБЧИСЛЕНИХ В ТОЧЦІ T.
C Y(I,J+1) МІСТИТЬ J - ТУ ПОХІДНИХ Y(I,1) ПОМНОЖЕНУ
C НА H**J/FACTORIAL(J).
C
C H ВЕЛИЧИНА ОСТАННЬОГО ВИКОРИСТАНОГО КРОКУ.
C
C ERROR МАСИВ РОЗМІРУ NYDIM, ЯКИЙ МІСТИТЬ ЗНАЧЕННЯ ПОХИБКИ
C ЧИСЕЛЬНОГО ІНТЕГРУВАННЯ НА ДАНОМУ КРОЦІ.
C
C JSTART ПОРЯДОК ТОЧНОСТІ ВИКОРИСТАНОГО МЕТОДУ (АБО ПОРЯДОК
C МАКСИМАЛЬНОЇ ПОХІДНОЇ).
C
C KFLAG INTEGER
C KFLAG=0 - КРОК УСПІШНИЙ,
C KFLAG=-1 - ТОЧНІСТЬ НЕ ДОСЯГНУТА ДЛЯ DABS(H)=HMIN,
C KFLAG=-2 - ЗБІЖНІСТЬ КОРЕКТОРА НЕ ДОСЯГНУТА ДЛЯ
C DABS(H)=HMIN.
C
C PW МАСИВ РОЗМІРУ LP
C (ЯКЩО MITER=0, ТО LP=1.
C ЯКЩО MITER=1,2, ТО LP=NYDIM*N.
C ЯКЩО MITER=3, ТО LP=N.
C ЯКЩО MITER=4,5, ТО LP=NYDIM*(2*ML+MU+1)).
C ВИКОРИСТОВУЄТЬСЯ ДЛЯ ЗБЕРІГАННЯ МАТРИЦІ ЯКОБІ.
C
C NSTEP КІЛЬКІСТЬ КРОКІВ ЧИСЕЛЬНОГО ІНТЕГРУВАННЯ (ВДАЛИХ І
C НЕВДАЛИХ).
C
C NFUN КІЛЬКІСТЬ ОБЧИСЛЕНЬ ПРАВИХ ЧАСТИН СИСТЕМИ.
C
C NJAC КІЛЬКІСТЬ ОБЧИСЛЕНЬ МАТРИЦІ ЯКОБІ.
C
C HUSED ВЕЛИЧИНА КРОКУ, ЯКИЙ ВИКОРИСТОВУВАВСЯ ДЛЯ
C ОБЧИСЛЕННЯ РОЗВ'ЯЗКУ В ДАНІЙ ТОЧЦІ T.
C
C NQUSED ПОРЯДОК МЕТОДУ, ЯКИЙ ВИКОРИСТОВУВАВСЯ ДЛЯ
C ОБЧИСЛЕННЯ РОЗВ'ЯЗКУ В ДАНІЙ ТОЧЦІ T.
C
C IWORK РОБОЧИЙ МАСИВ РОЗМІРУ NYDIM.
C
C FSAVE РОБОЧИЙ МАСИВ РОЗМІРУ 2*NYDIM.
C
C
C ПІДПРОГРАМА STIFF ВИКОРИСТОВУЄ ТАКІ ПІДПРОГРАМИ:
С ------------------
C MATDEC - ЗДІЙСНЮЄ LU РОЗКЛАД МАТРИЦІ.
C
C MATSOL - РОЗВ'ЗУЄ СИСТЕМИ ЛІНІЙНИХ РІВНЯНЬ З ТРИКУТНИМИ
C МАТРИЦЯМИ (ЯКЩО MITER=1 АБО 2).
C
C PREDIC - ОБЧИСЛЮЄ ПОЧАТКОВЕ НАБЛИЖЕННЯ.
C
C COSET - ВИЗНАЧАЄ ЗНАЧЕННЯ КОЕФІЦІЄНТІВ ВІДПОВІДНОГО
C БАГАТОКРОКОВОГО МЕТОДУ.
C
C RESCAL - ПЕРЕМАШТАБОВУЄ МАСИВ ФУНКЦІЙ ПРИ ЗМІНІ КРОКУ
C ЧИСЕЛЬНОГО ІНТЕГРУВАННЯ.
C
C BANMAT - РОЗВ'ЗУЄ СИСТЕМУ ЛІНІЙНИХ РІВНЯНЬ ЗІ СТРІЧКОВОЮ
C МАТРИЦЕЮ (ЯКЩО MITER=4 АБО 5).
C
C ADDVEC - ДОДАЄ ДВА ВЕКТОРИ.
C
C DOTPRO - ПІДПРОГРАМА-ФУНКЦІЯ, ЯКА СКАЛЯРНО ПЕРЕМНОЖАЄ ДВА
C ВЕКТОРИ.
C
С
C ПІДПРОГРАМА STIFF ВИКОРИСТОВУЄ ТАКОЖ ПІДПРОГРАМИ, ЯКІ
С ПОВИНЕН НАПИСАТИ КОРИСТУВАЧ:
С ------------------
C DIFFUN ОБЧИСЛЮЄ МАСИВ ПРАВИХ ЧАСТИН СИСТЕМИ ЗДР
C SUBROUTINE DIFFUN(N,T,Y,YDOT)
C IMPLICIT REAL*8 (A-H,O-Z)
C REAL*8 Y(N),YDOT(N)
C YDOT(1)=...
C ...
C RETURN
C END
C
C PEDERV ОБЧИСЛЮЄ МАТРИЦЮ ЯКОБІ ПРАВИХ ЧАСТИН СИСТЕМИ ЗДР
C SUBROUTINE PEDERV(N,T,Y,YDOT)
C IMPLICIT REAL*8 (A-H,O-Z)
C REAL*8 Y(N)
С COMMON /STCOM4/ PW(1)
C PW(1)=...
C ...
C RETURN
C END
С ЯКЩО MITER=1, ТО МАТРИЦЯ ЯКОБІ ЗАДАЄТЬСЯ В
C ОДНОВИМІРНОМУ МАСИВІ РОЗМІРУ NYDIM*N І РОЗМІЩУЄТЬСЯ
C ПО СТОВПЦЯХ. У ВИПАДКУ СТРІЧКОВОЇ МАТРИЦІ ЯКОБІ
C (MITER=4) ВОНА ЗАДАЄТЬСЯ У ВИГЛЯДІ МАТРИЦІ РОЗМІРУ
C NYDIM*(MU+LU+1) І РОЗМІЩУЄТЬСЯ В ОДНОВИМІРНОМУ
C МАСИВІ РОЗМІРУ NYDIM*(2*MU+LU+1) ПО СТОВПЦЯХ.
C -------------------------------------------------------------
Приклад 1. Нехай на відрізку потрібно чисельно розв’язати задачу Коші
з точністю та початковим кроком інтегрування . Матриця Якобі правих частин цієї ситеми ЗДР має вигляд
Оскільки власні значення матриці Якобі , то задача є жорсткою.
Головна програма та підпрограми DIFFUN і PEDERV для прикладу 1:
PROGRAM EXAMPLE1
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION Y(10,13),YMAX(10),ERROR(10),PW(100),
1 FSAVE(20),IWORK(10)
COMMON/STCOM1/T,H,HMIN,HMAX,EPS,N,MF,KFLAG,JSTART,MAXORD
COMMON/STCOM2/HUSED,NQUSED
COMMON/STCOM3/ML,MU
COMMON/STCOM4/NSTEP,NFUN,NJAC
NYDIM=10
EPS=1.D-2
KB=0
401 CONTINUE
N=3
T=0.0D0
TEND=1.D0
Y(1,1)=1.D0
Y(2,1)=0.D0
Y(3,1)=0.D0
H=3.3D-8
HMAX=TEND
HMIN=1.D-15
JSTART=0
MF=22
MAXORD=5
WRITE(0,20) MF,EPS
20 FORMAT(//3X,'MF=',I2/,' EPS='D11.3)
NSTEP=0
NFUN=0
NJAC=0
DO 30 I=1,N
30 YMAX(I)=1.D0
40 CONTINUE
CALL STIFF(Y,YMAX,ERROR,PW,FSAVE,IWORK,NYDIM)
IF(KFLAG.EQ.0)GO TO 60
WRITE(0,50) KFLAG
50 FORMAT(/' KFLAG=',I2/)
STOP
60 CONTINUE
IF(DABS(TEND-T).LE.1.D-15) GO TO 90
IF(TEND-T-H) 80,40,40
80 E=TEND-T
S=E/H
DO 85 I=1,N
DO 85 J=1,JSTART
85 Y(I,1)=Y(I,1)+Y(I,J+1)*S**J
T=T+E
GO TO 60
90 CONTINUE
WRITE(0,556) H,T,(Y(I,1),I=1,N)
556 FORMAT(1X,5D16.8)
WRITE(0,95) NSTEP,NFUN,NJAC
95 FORMAT(/' NSTEP=',I4,' NFUN= ',I5,' NJAC=',I4)
KB=KB+1
IF(KB.GE.3) GO TO 402
EPS=EPS*1.D-2
GO TO 401
402 CONTINUE
STOP
END
SUBROUTINE DIFFUN (N,T,Y,YDOT)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION Y(1),YDOT(1)
YDOT(1)=-Y(1)+1.D8*Y(3)*(1.D0-Y(1))
YDOT(2)=-1.D1*Y(2)+3.D7*Y(3)*(1.D0-Y(2))
YDOT(3)=-YDOT(1)-YDOT(2)
RETURN
END
SUBROUTINE PEDERV(N,T,Y,PW,NYDIM)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION Y(1),PW(1)
PW(1)=-1.D0-1.D8*Y(3)
PW(2)=0.D0
PW(3)=-PW(1)
PW(NYDIM+1)=0.D0
PW(NYDIM+2)=-1.D1-3.D7*Y(3)
PW(NYDIM+3)=-PW(NYDIM+2)
N2=NYDIM*2
PW(N2+1)=1.D8*(1.D0-Y(1))
PW(N2+2)=3.D7*(1.D0-Y(2))
PW(N2+3)=-PW(N2+1)-PW(N2+2)
RETURN
END
Отримані результати розв’язку прикладу 1 в кінці відрізка програмою STIFF для такі: .
Таблиця 1. Характеристики програми STIFF
NSTEP
NFUN
NJAC
EPS=
22
50
6
EPS=
40
94
8
EPS=
70
156
11
ПОСЛІДОВНІСТЬ ВИКОНАННЯ ЛАБОРАТОРНОЇ РОБОТИ
1. Одержати варіант завдання.
2. Вивчити відповідний лекційний матеріал.
3. Написати головну програму та підпрограми DIFFUN, PEDERV.
4. Провести відлагодження програми і отримати розв’язок задачі при MF=21,22 з точностями .
5. Вивести на друк в точці TEND значення чисельного розв’язку, а також основні характеристики програми NSTEP,NFUN,NJAC.
ЗМІСТ ЗВІТУ
1. Постановка задачі (конкретний варіант).
2. Текст головної програми та підпрограм DIFFUN, PEDERV.
Ділись своїми роботами та отримуй миттєві бонуси!
0%
Оголошення від адміністратора