МІНІСТЕРСТВО ОСВІТИ І НАУКИ Р.Ф.
Курганський державний університет
Кафедра прикладної та вищої математики
Лабораторна робота № 43
на тему:
Рішення змішаної задачі для рівняння
гіперболічного типу методом сіток
Група М-2136
Курган 1998
Розглянемо змішану задачу для хвильового рівняння ( В¶ 2 u/ В¶ t 2 ) = c 2 * ( В¶ 2 u/ В¶ x 2 ) (1). Завдання полягає у відшуканні функції u (x, t) задовольняє даним рівнянню при 0 ВЈ T, початковим умовам u (x, 0) = f (x), В¶ u (x, 0)/ В¶ t = g (x), 0 ВЈ x ВЈ a і нульовими крайовими умовами u (0, t) = u (1, t) = 0.
Так як заміна змінних t В® ct приводить рівняння (1) до виду ( В¶ 2 u/ В¶ t 2 ) = ( В¶ < sup> 2 u/ В¶ x 2 ), то в подальшому будемо вважати з = 1.
Для побудови різницевої схеми рішення задачі будуємо в області D = {(x, t) | 0 ВЈ x ВЈ a, 0 ВЈ t ВЈ T } сітку x i = ih, i = 0,1 ... n, a = h * n, t j = j * ttt, j = 0,1 ... , M, t m ​​= T і апроксимуємо рівняння (1) в кожному внутрішньому вузлі сітки на шаблоні типу "хрест".
t
T
j +1
j
j-1
0 i-1 i i +1
Використовуючи для апроксимації приватних похідних центральні різницеві похідні, отримуємо наступну різницеву апроксимацію рівняння (1).
u i, j +1 - 2u ij + u i, j-1 u i +1,, j - 2u ij + u i-1, j
t 2 h 2
(4)
Тут u ij - наближене значення функції u (x, t) у вузлі (x i , t j ).
Вважаючи, що l = t/h, отримуємо тришарову різницеву схему
u i, j +1 = 2 (1 - l 2 ) u i, j + l 2 (u i +1, j - u i-1, j ) - U i, j-1 , i = 1,2 ... n. (5)
Для простоти в даній лабораторній роботі задані нульові граничні умови, тобто m 1 (t) Вє 0, m 2 (t) Вє 0. Значить, в схемі (5) u 0, j = 0, u nj = 0 для всіх j. Схема (5) називається тришарової на трьох часових шарах з номерами j-1, j, j +1. Схема (5) явна, тобто дозволяє в явному вигляді висловити u i, j через значення u з попередніх двох шарів.
Чисельне рішення задачі полягає в обчисленні наближених значень u i, j рішення u (x, t) у вузлах (x i , t j ) при i = 1, ... n, j = 1,2, ... , M. Алгоритм рішення заснований на тому, що рішення на кожному наступному шарі (j = 2,3,4, ... n) можна отримати перерахунком рішень з двох попередніх шарів (j = 0,1,2, ..., n-1) за формулою (5). На нульовому часовому шарі (j = 0) рішення відомо з початкової умови u i0 = f (x i ).
Для обчислення рішення на першому шарі (j = 1) в даній лабораторній роботі прийнятий найпростіший спосіб, що складається в тому, що якщо покласти В¶ u (x, 0)/ В¶ t В»(u (x, t) - u (x, 0))/t (6), то u i1 = u i0 + + t (x i ), i = 1,2, ... n. Тепер для обчислення рішень на наступних шарах можна застосовувати формулу (5). Рішення на кожному наступному шарі виходить перерахунком рішень з двох попередніх шарів за формулою (5).
Описана вище схема апроксимує задачу з точністю до О (t + h 2 ). Невисокий порядок апроксимації по t пояснюється використанням занадто грубою апроксимації для похідної по е в формулі (6).
Схема стійка, якщо виконано умову Куранта t
Недолік схеми в тому, що як тільки вибраних величина кроку сітки h у напрямку x, з'являється обмеження на величину кроку t по змінній t. Якщо необхідно провести обчислення для великого значення величини T, то може знадобитися велика кількість кроків по змінній t. Зазначений гнедостаток характерний для всіх явних різницевих схем.
Для оцінки похибки рішення зазвичай вдаються до методам згущення сітки.
Для вирішення змішаної задачі для хвильового рівняння по явної різницевої схемою (5) призначена частину програми, позначена Subroutine GIP3 Begn ... End. Ця підпрограма обчислює рішення на кожному шарі за значеннями рішення з двох попередніх шарів.
Вхідні параметри:
hx - крок сітки h по змінній х;
ht - крок сітки t по змінній t;
k - кількість вузлів сітки по x, a = hn;
u1 - масив з k дійсних чисел, що містить значення рішень на (j - 1) часовому шарі, j = 1, 2, ...
u2 - масив з n дійсних чисел, що містить значення рішень на j - му часовому шарі, j = 1, 2, ...
u3 - робочий масив з k дійсних чисел.
Вихідні параметри:
u1 - масив з n дійсних чисел, що містить значення рішення з j - м часовому шарі, j = 1, 2, ...
u2 - масив з n дійсних чисел, що містить значення рішення з (j +1) - м часовому шарі, j = 1, 2, ... .
До частини програми, позначеної як Subroutine GIP3 Begin ... End відбувається циклічне звернення, пеоред першим зверненням до програмі елементів масиву u2 присвоюються початкові значення, а елементам масиву u1 - значення на вирішення на першому шарі, вичіслінние за формулами (6). При виході з підпрограми GIP3 в масиві u2 знаходиться значення рішення на новому часовому шарі, а в масиві u1 - значення рішення на попередньому шарі.
Порядок роботи програми:
1) опис масивів u1, u2, u3;
2) привласнення фактичних значень параметрів n, hx, ht, облюдая умова Куранта;
3) привласнення початкового значення рішення елементам масиву і обчислене за формулами (6) значення рішення на першому шарі;
4) звернення до GIP3 в циклі k-1 раз, якщо потрібно знайти рішення на k-му шарі (k Ві 2).
Приклад:
1
0.5 0.5
Вирішити завдання про коливання струни одиничної довжини з закріпленими кінцями, початкове положення якої зображено на малюнку. Початкові швидкості дорівнюють нулю. Обчислення виконати з кроком h по x, рівним 0.1, з кроком t по t, рівним 0.05, провести обчислення для 16 часових шарів з друком результатів на кожному шарі. Таким чином, завдання має вигляд
(В¶ 2 u/В¶ t 2 ) = (В¶ 2 u/В¶ x 2 ) , x ГЋ [0, 1], t ГЋ [0 , T],
u (x, 0) = f (x), x ГЋ [0, a], В¶ u (x, 0)/В¶ t = g (x), x ГЋ [0, a],
u (0, t) = 0, u (1, t) = 0, t ГЋ [0, 0.8],
Г¦ 2x, x ГЋ [0, 0.5],
f (x) = Г g (x ) = 0
Г® 2 - 2x, x ГЋ [0.5, 1 ],
Будуємо сітку з 11 вузлів по x і виконуємо обчислення для 16 шарів по t. Програма, та результати обчислення наведені далі.