ЛЕКЦІЯ № 10
ЧИСЛОВЕ РІШЕННЯ ДИФЕРЕНЦІАЛЬНИХ РІВНЯНЬ З ЧАСТИННИМИ ПОХІДНИМИ
Класифікація диференціальних рівнянь в частинних похідних
Найбільш розробленими є методи наближеного розв’язання для диференціальних рівнянь в частинних похідних другого порядку з двома незалежними змінними. Для розв’язання багатьох практичних задач необхідно розглядати так звані лінійні диференціальні рівняння в частинних похідних, тобто диференціальні рівняння першого степеня відносно шуканої функції та всіх її похідних і які не містять їх добутків. Такі рівняння можна подати в наступному вигляді:
(10.1)
У рівнянні (10.1) шуканою є функція z, а х і у – незалежними змінними. Функції А (х, у), В (х, у), С (х, у), а (х, у), b (х, у), с (х, у) – неперервні функції від х і у такі, що мають неперервні частинні похідні.
Проведемо класифікацію диференціальних рівнянь в частинних похідних, засновану на розгляді рівняння (10.1). Уведемо для зручності позначення
і розглянемо спрощену форму рівняння (10.1):
(10.2)
яка відповідає рівнянню (10.1) при 
Рівняння (10.2) завжди може бути приведено до однієї з трьох стандартних канонічних форм. Цими формами є еліптичні, параболічні і гіперболічні диференціальні рівняння в частинних похідних. Тип рівняння визначається значенням коефіцієнтів в рівнянні (10.2) і пов'язаний зі знаком дискримінанта 
в рівнянні (10.2).
В залежності від знака дискримінанта маємо:
Δ < 0 еліптичний тип в точці (х, у);
Δ = 0 – параболічний тип в точці (х, у);
Δ > 0 – гіперболічний тип в точці (х, у).
Якщо коефіцієнти А, В і С є величинами сталими, не залежними від х і у, то канонічні рівняння є повністю еліптичними, параболічними або гіперболічними.
Скінченнорізницеві апроксимації
Скінченнорізницеві апроксимації для частинних похідних є найбільш поширеним підходом до числового інтегрування диференціальних рівнянь у частинних похідних. Частинні похідні замінюються відповідними різницевими співвідношеннями по відповідним незалежним змінним. У загальному випадку розмірність області, у якій необхідно знайти розв’язання диференціального рівняння у частинних похідних, дорівнює кількості незалежних змінних. У випадку двох незалежних змінних х і у область є двовимірною.
Метод, що застосовується для скінченнорізницевої апроксимації, базується на покритті області сіткою прямокутних кліток шириною h (в напрямку осі Ох) і висотою k (в напрямку осі Оу). Величина залежної змінної 
встановлюється в будь-якій точці в межах області. Зокрема, коли задана одна точка прямокутної сітки з координатами 
, то чотири точки, що її оточують мають координати 
; 
; 
; 
. Геометричний спосіб покриття області сіткою показаний на рис. 10.1.
Застосовуючи скінченнорізницеву схему, можна похідні подати наступним чином
 |
Для наочного подання цих похідних зручно застосовувати шаблони, які мають вигляд:
У центрах кіл цих шаблонів указуються коефіцієнти диференціального рівняння. Коло центральної частини шаблону відповідає величині 
. Додатним приростам по горизонтальним лініям (рис. 10.1) відповідає лівий кінець, від'ємним – правий кінець шаблону.
Розглянемо застосування викладеного підходу до скінченнорізницевої апроксимації частинних похідних. Для функції застосуємо наступне позначення у вузлах сітки
Оскільки 
означає похідну від z по х при сталому у, отримуємо
.
При цьому обчислювальний шаблон має вигляд
 |
Цей шаблон містить тільки горизонтальні елементи, оскільки s (або у) є величинами сталими.
Аналогічно можна отримати шаблон для похідної в напрямку s.
 |
Змішана похідна виходить аналогічно:
 |
Множенням окремих елементів рядка на кожний елемент стовпця отримаємо
10.3. Розрахунок поля температур в однорідній штабі за явною схемою
В якості ілюстрації числових методів розв’язання диференціальних рівнянь розглянемо числове розв’язання диференціального рівняння параболічного типу, пов’язаного з поширенням тепла в однорідній штабі.
Процес теплопровідності виникає, якщо тіло є неоднорідно нагрітим. Найпростіша для вивчення теплопровідності система – лінійна однорідна штаба. У простій моделі бічні поверхні штаби вважаються теплоізольованими, тобто через них немає обміну теплом між штабою і навколишнім середовищем.
Отже, визначимо розподіл температури 
в металевій штабі для будь-якого значення часу 
за умови, що фізичні властивості матеріалу штаби є величинами сталими.
У виразі під 
розуміється температура, під х – поздовжня координата, а під 
– час.
Розподіл температури в металевій штабі може бути описаним наступним диференціальним рівнянням з початковими та граничними умовами:
(10.3)
(10.4)
де А – стала величина, яка залежить від теплопровідних властивостей матеріалу, з якого виготовлено штабу, а також від деяких інших чинників.
Початкова умова задає розподіл температури в штабі в початковий момент часу (його вважають рівним нулю). Граничні умови (у нашому випадку їх має бути дві) у цьому простішому варіанті вказують на те, яка температура підтримується на кінцях штаби.
Зазначимо, що супроводження рівняння теплопровідності початковими і граничними умовами, робить постановку задачі фізично однозначною. Нагадаємо, що якщо для диференціального рівняння задані початкові умови (умови в початковий момент часу), то така задача називається задачею Коші, якщо ж задані граничні умови (на границях досліджуваної області), то така задача називається крайовою задачею, якщо задані і початкові, і граничні умови, то маємо змішану крайову задачу.
Поставлену задачу краще розв’язувати в безрозмірних одиницях. На рис. 10.2 показано металеву штабу одиничної довжини.
Для числового розв’язання поставленої задачі апроксимуємо диференціальне рівняння (10.3) разом з початковими і граничними умовами (10.4) скінченними різницями.
 |
Моделювання процесу теплопровідності пов'язано з дискретизацією як зміни температури за часом, так і вздовж штаби.
Для цього побудуємо сітку в області, обмеженій 
і 
(рис. 10.3), з кроком 
і 
, де M і N – деякі цілі числа. Під функцією 
розуміють апроксимацію рівняння (10.3) скінченними різницями.
  |
Для будь-якого вузла сітки (i, j), для якого лічильник і не дорівнює 0 та M, а лічильник j не дорівнює 0, похідні диференціального рівняння (10.3) можна замінити скінченними різницями наступним чином:
Тоді рівняння (10.3) наближено замінюється виразом:
Звідки отримують
(10.5)
Нас буде цікавити температура в точці 
, тому розв'яжемо рівняння (10.5) відносно температури в точці
(10.6)
де
Таким чином, із диференціального рівняння (10.3) отримано зв'язок між температурою наступного (j + 1) і поточного шарів (j).
 |
На рис. 10.4 хрестиками і маленькими колами позначені ті вузли сітки, які задіяні при визначенні похідних у рівнянні (10.6).
Якщо для будь-якого рівня часу 
відомі значення температури 
, то з рівняння (10.6) можна безпосередньо визначити значення температури при 
, тобто в шарі (j + 1). Подібна схема розрахунку має назву явної.
Слід відмітити, що сукупність вузлів у фіксований момент часу називається шаром.
Але виникають питання. Де взяти значення температури в нульовому шарі, тобто при 
? Зазначимо, що для визначення температури в нульовому шарі застосовуються початкові умови, які записуються у вигляді:
(10.7)
Подібні ж питання постають для граничних точок, що мають значення лічильника і, які дорівнюють 0 і М. Виходячи з граничних умов, можна зробити висновки, що
(6.6)
Застосувавши початкові умови у вигляді (6.5) і граничні умови у вигляді (6.6), можна за рівнянням (6.4) розрахувати значення температури у шарі (j = 1) і продовжити розрахунки до шару з номером N, тобто поступово збільшуючи час прогріву штаби.
Розглянемо процес зміни температури у штабі з паралельними сторонами 
за умови, що відомий коефіцієнт α її теплопровідності та в початковий момент часу 
температура 
рівномірно розподілена вздовж штаби. Дві взаємно протилежні сторони знаходяться при постійній температурі 
. Необхідно знайти розподіл температури залежно від часу та координати розташування пробної точки.
Рівняння теплопровідності матиме наступний вигляд:
. (6.7)
Уведемо для розгляду безрозмірні величини
Тоді рівняння (6.7) може бути переписаним і поданим у наступному вигляді:
яке будемо розв’язувати за наявності таких умов:
За цих обставин початкові і граничні умови приймуть доволі просту форму:
Слід зазначити, що модифікацією функціональних залежностей 
для початкових і граничних умов можна розрахувати відповідне поле температур.
 |
Результати розрахунків поля температур у безрозмірному вигляді показані на рис. 6.4. На цьому рисунку маленькі кола відповідають розрахунковим точкам. Криві побудовані в інтервалі безрозмірного часу від 0,05 до 0,35 з кроком 0,05. Із наведених даних випливає, що зі збільшенням часу температура зростає, наближаючись до одиниці.
6.5. Розрахунок поля температур в однорідній штабі за неявною схемою
У вище розглянутому прикладі застосування явного методу розв’язання диференціального рівняння (6.1) апроксимовані температури у вузлах сітки залежали тільки від 
i 
. З розгляду рис. 6.5 випливає, що за цим підходом визначення температури v можливе тільки в зоні, позначеній літерою А. Саме вузли цієї зони задіяні при визначенні розподілу температури в штабі. У той же час відомо, що для визначення поля температур, обумовленого диференціальним рівнянням (6.1) та початковими і граничними умовами (6.2), необхідно проводити розрахунки як зоні А, так і в зоні В.
Явна різницева схема зручна в застосуванні але стійка лише при виконанні умови
Цю обставину треба враховувати при виборі кроку за часом і довжиною штаби.
 |
Розв’яжемо рівняння (6.1) з початковими і граничними умовами (6.2) за дещо іншою схемою апроксимаційного подання похідних, яку називають неявною схемою. Запишемо скінченні різниці похідних у наступному вигляді:
(6.8)
 |
На відміну від попереднього прикладу, коли друга похідна визначалася в j -му шарі з часом
, зараз вона знаходиться із значеннями апроксимованої температури в шарі j +1, якому відповідає час . Це призводить до іншої схеми розташування чотирьох вузлових точок (рис. 6.6), задіяних в обчислювальному процесі.
Розв’яжемо рівняння (6.8) відносно , яке запишемо в правій частині розв’язку. Будемо мати вираз:
. (6.9)
Коефіцієнт λ має той же сенс, що й у формулі 6.4.
Отже, на відміну від вище розглянутої явної схеми кожне різницеве рівняння другої схеми містить на кожному новому шарі три невідомі значення, які неможливо визначити відразу ж, як це має місце в явній схемі. При цьому друга різницева схема складається з лінійних триточкових рівнянь, тобто кожне рівняння містить невідому функцію в трьох точках нового шару.
Граничні та початкові умови залишаються тими ж самими, що й у вище розглянутому прикладі:
,
,
.
Рівняння (6.9) можуть бути записані для кожної вузлової точки з діапазону 
для будь-якого моменту часу, що призведе до системи з М – 1 рівнянь з М – 1 невідомими 
.
Розглянемо розв’язання наступної системи лінійних рівнянь:
(6.10)
Цю систему рівнянь зручніше записати в такому вигляді:
(6.11)
У системах (6.10) і (6.11) для скорочення запису індекси (j + 1) при всіх температурах v пропущені. Праві частини рівнянь (6.10), які є відомими величинами, для спрощення запису позначені як 
, де 
Система рівнянь (6.11) має назву тридіагональної матриці, оскільки заповненими є тільки три діагональні елементи кожного рядка. Цю систему можна розв’язати відомим методом Гауса, але при цьому мають бути заповнені всі нульові члени, тому матриця буде мати значний розмір і буде виконано багато зайвих обчислювальних операцій, що є недоцільним.
Краще для розв’язання тридіагональної матриці застосовувати так званий метод прогонки в наступній формі:
. (6.12)
Коефіцієнти b i і γ i будуть розглянуті трошки нижче.
Підстановка виразу (6.12) до системи (6.11) призведе до наступного запису:
,
з якого можна визначити vi:
.
Це дає можливість застосувати наступну рекурсивну схему:
;
.
З першого рівняння системи (6.11) маємо
,
де 
.
Коефіцієнт γ1 дорівнює 
Із останнього рівняння системи (6.11) можна отримати
.
Звідки
.
Таким чином, алгоритм розв’язання тридіагональної матриці у стислому вигляді запишеться наступним чином:
в якому
Таким чином, у разі неявної схеми, щоб розрахувати значення функції температур в кожен наступний момент часу, тобто, щоб перейти на наступний шар за часом, необхідно кожного разу вирішувати методом прогонки лінійну систему.
Результати розрахунків за неявною схемою в графічному вигляді наведені на рис. 6.7. Їм відповідають криві без додаткових позначень. На цьому ж рисунку для порівняння показані результати розрахунків поля температур в однорідній штабі, отримані із застосуванням явної схеми. Це криві з нанесеними на них точками, вони повністю відповідають графічним даним, наведеним на рис. 6.4.
Як випливає з розгляду рис. 6.7, результати розрахунків проведені за явною і неявною схемами дуже близькі, що і очікувалося.
Зрозуміло, що коло задач, для розв’язання яких застосовується числове математичне моделювання, не обмежується розрахунками температурних полів. Існує безліч задач, рішення яких можна отримати ретельно підібраними числовими методами.
 |
На завершення треба відзначити, що сучасні числові методи і потужні ЕОМ надали можливість розв’язувати такі задачі з числового математичного моделювання, про які раніше можна було тільки мріяти. Але застосовувати числові методи не так просто. Треба пам’ятати, що ЕОМ уміють виконувати тільки арифметичні дії та логічні операції. Тому окрім розробки математичної моделі необхідна ще розробка відповідного алгоритму, який зводить всі розрахунки до послідовності арифметичних і логічних дій. Вибирати числову математичну модель та алгоритм треба з урахуванням швидкодії комп’ютера та об’єму його пам’яті: занадто складна модель може виявитися комп’ютеру не під силу, а дуже проста – не відповідатиме фізичній точності об’єкта, що вивчається засобами числового математичного моделювання.
Безпосередньо алгоритм і програма для ЕОМ мають бути ретельно перевірені. Але перевірка програми справа нелегка, про що свідчить популярне твердження: "В будь-якій як завгодно малій програмі є щонайменше одна помилка". Про це треба обов’язково пам’ятати, розв’язуючи задачі числового математичного моделювання.
