Погрішність вирішення й збіжність

Як відомо, збіжність вирішення задач, отриманих сітковими різницевими методами, визначається формулою:

Збіжність = Апроксимація + Стійкість.

Як було показано вище, крайова задача апроксимується різницевою схемою із другим порядком точності щодо кроків сітки. Стійкість різницевої схеми для рівнянь еліптичного типу випливає із принципу максимуму, який виконується при певних обмеженнях на коефіцієнти рівняння [1]. Отже, буде мати місце збіжність, що означає:

1) система лінійних алгебраїчних рівнянь, що отримана при апроксимації крайової задачі різницевою задачею, буде мати єдиний розв'язок;

2) при зменшенні кроку сітки розв'язок різницевої задачі буде наближуватись до точного розв'язку диференціальної крайової задачі.

Практична оцінка погрішності розв'язку крайової задачі може бути отримана за правилом Рунге.

Нехай у загальному випадку є формула для обчислення за значеннями на рівномірній сітці із кроком , а залишковий член цієї формули має таку структуру:

.

Проведемо тепер розрахунок по тій же наближеній формулі для тієї ж точки , але використовуючи рівномірну сітку з іншим кроком . Тоді одержимо:

.

Ясно, що .

Віднімаємо ці дві формули й одержуємо першу формулу Рунге:

(2.7)

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

(2.8)

Застосовуючи наведені формули до розв'язку задачі Діріхле для рівняння Лапласа методом сіток, бачимо, що необхідно розв'язати задачу двічі – використовуючи сітку із кроком , а потім, використовуючи сітку із кроком в 2 рази менше (помітимо, що точок сітки стане в 4 рази більше). Порядок апроксимації , отже, у знаменнику дроби (2.7) буде число 3. У чисельник підставляємо різницю розв'язків, що отримані у відповідних точках при розрахунку з різними кроками.

Задача для самостійного вирішення на практиці

Задача 2.1. Розв'язати рівняння в області із граничними умовами .

Завдання до лабораторної роботи

Завдання 2.1. Розв'язати задачу Діріхле для рівняння Лапласа у квадраті з вершинами , , , з граничними умовами, наведеними в таблиці 2.2. Обрати сітку з кількістю точок розбивки по кожній зі змінних. Знайти наближене вирішення задачі із заданою точністю розв'язку системи алгебраїчних рівнянь .

Завдання 2.2. Виконати завдання 2.1 для .

Завдання 2.3. Оцінити погрішність вирішення задачі для обраних сіток, використовуючи правило Рунге.

Таблиця 2.2 – Варіанти індивідуальних завдань

Номер варіанта
			40 y
			30
			30

2.7. Контрольні питання

1. Які види сіток використовуються в методі кінцевих різниць? Яким чином будують на цих сітках різницеві апроксимації рівнянь і відповідні їм шаблони?

2. Які прямі й ітераційні методи використовують для розв'язку систем алгебраїчних рівнянь у задачах з диференційними рівняннями із частинними похідними (ДРЧП)?

3. Дайте характеристику ітераційних методів, що використовуються для розв'язку систем алгебраїчних рівнянь у задачах з ДРЧП.

4. Як задаються граничні умови? Яким чином задається початкове наближення при вирішенні ДРЧП із використанням ітераційних методів? Відповідь поясніть на прикладі.

5. З яких міркувань вибирають крок сітки в методі кінцевих різниць?

МЕТОД СІТОК

РОЗВ'ЯЗКУ крайових задач для РІВНЯНЬ ПАРАБОЛІЧНОГО ТИПУ

Постановка задачі

Розглянемо розв'язок змішаної крайової задачі для диференційних рівнянь з частинними похідними (ДРЧП):

, , (3.1)

с початковою умовою:

,

і граничними умовами:

, , .

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

Для розв'язку задачі методом кінцевих різниць побудуємо прямокутну сітку (рис. 3.1), вузли якої визначаються формулами:

, , , .

Значення на лівій і нижній границях сітки відомі з початкових і граничних умов. На правій границі відоме значення частинної похідної розв'язку рівняння по змінній .

Замінимо частинні похідні в рівнянні теплопровідності їх кінцево-різницевими апроксимаціями в кожному внутрішньому вузлі:

, (3.2)

(3.3)

 

 

Рисунок 3.1 – Просторово - часова сітка й шаблони, що використовуються для розв'язку рівнянь параболічного типу

Вираз (3.2) можна вважати наближенням похідної як у точці так і в точці порядку .

Вираз (3.3) апроксимує похідну з порядком . Таким чином, порядок апроксимації диференційного рівняння

Для розв'язку змішаної крайової задачі необхідно апроксимувати похідну в граничній умові на правому кінці:

.

Використовуючи кінцево-різницеву апроксимацію, одержуємо:

(3.4).

Порядок апроксимації останньої формули .

Явна різницева схема

Підставимо вирази (3.2) і (3.3) у рівняння (3.1) і розв'яжемо його щодо значень функції на верхньому часовому шарі:

. (3.5)

Формула (3.5) вирішує поставлену задачу, оскільки вона виражає розв'язок у момент часу i+1 через розв'язок у момент часу i.

З (3.4) знаходимо: . (3.6)

Алгоритм обчислень за явною схемою реалізується наступною послідовністю дій.

1. Обчислюємо значення сіткової функції на першому часовому шарі з початкових умов: .

2. Знаходимо розв'язок на сітковому шарі , використовуючи явну формулу:

, , .

3. Обчислюємо величину за формулою (3.6) .

Завершивши кроки 1-3 одержуємо розв'язок при . Для отримання розв'язку при повторюють кроки 2,3, збільшивши на одиницю й використовуючи з попереднього рядка. Явна схема відповідає верхньому шаблону, наведеному на рисунку 3.1.

Недолік явної схеми: якщо крок за часом виявляється досить великим у порівнянні із кроком по , то погрішності обчислень можуть стати настільки великими, що отриманий розв'язок втрачає сенс, тобто розв'язок стає нестійким. Для стійкості явної схеми повинна виконуватися умова ,

яка накладає досить жорсткі обмеження на крок за часом () і веде до значного збільшення часу для отримання рішення. Така схема називається умовно стійкою.

Приклад 3.1. Обчислити за допомогою явної схеми наближений розв'язок змішаної задачі

, , ,

с початковою умовою й граничними умовами .

Розв'язок. Формула (3.5) з урахуванням кроків по - і по - після перетворень прийме вигляд:

.

Схема стійка, якщо . Перевіримо: , тобто обрана різницева схема є стійкою.

Рисунок 3.2 – Просторово - часова сітка для прикладу 3.1

З початкових умов одержимо:

, , , , , ;

із граничних умов одержимо

, ; , ; , .

Підрахуємо значення для :

;

Аналогічно обчислюємо для інших значень . Отримуємо такі значення роподілу температури:

u=1.0000 0.6000 0.2000 0.2000 0.6000 1.0000

1.0000 0.6000 0.4000 0.4000 0.6000 1.0000

1.0000 0.7000 0.5000 0.5000 0.7000 1.0000

1.0000 0.7500 0.6000 0.6000 0.7500 1.0000

1.0000 0.8000 0.6750 0.6750 0.8000 1.0000

1.0000 0.8375 0.7375 0.7375 0.8375 1.0000

Графіки розв'язку для всіх моментів часу наведені на рис. 3.3. Початковий розподіл температури позначений на графіку неперервною лінією.

Щоб отримати оцінку погрішности рішення, проведемо розрахунки з кроком по координаті в два рази менше, тобто , та з кроком за часом в чотири рази менше, тобто . (Рішення з кроком за часом виявилось нестійким, дуже швидко зростала погрішність!)

Програма, що виконує розрахунки:

 

% приклад 5.1

n=11;hx=0.1; m=2*(n-1)+1; ht=0.005;

u=zeros(m,n); a=1; sig=a^2*ht/hx^2;

x=0:hx:1;

u(1,:)=abs(2*x-1);

t=0:ht:1;

v=abs(2*t-1);

plot(t,v, 'LineWidth',2); hold on; grid on

for i=1:m-1

u(i+1,1)=1; u(i+1,n)=1;

for j=2:n-1

u(i+1,j)=u(i,j)+sig*(u(i,j-1)-
2*u(i,j)+u(i,j+1));

end

end

plot(x,u(5,:),'--o','LineWidth',2)

plot(x,u(9,:),':s','LineWidth',2)

plot(x,u(13,:),'-.*','LineWidth',2)

plot(x,u(17,:),'--<','LineWidth',2)

plot(x,u(21,:),'-->','LineWidth',2)

% виводимо значення температури для кожного четвертого % шару за часом та кожне друге значення по координаті х

u(1:4:m,1:2:n)

Рисунок 3.3 – Графіки розв'язку прикладу 3.1 ().

Одержуємо такі значення температури:

u = 1.0000 0.6000 0.2000 0.2000 0.6000 1.0000

1.0000 0.6250 0.3750 0.3750 0.6250 1.0000

1.0000 0.6875 0.4922 0.4922 0.6875 1.0000

1.0000 0.7437 0.5850 0.5850 0.7437 1.0000

1.0000 0.7902 0.6605 0.6605 0.7902 1.0000

1.0000 0.8283 0.7222 0.7222 0.8283 1.0000,

та графік розвязку на рис. 3.4.

Рисунок 3.4 – Графіки розв'язку прикладу 3.1 з кроками

Оцінимо погрішність для часу t=0.1 (остання строка матриці).

,

Відносна погрішність становить трохи більше 2% від значення функції. У більшості випадків таку погрішність можна вважати припустимою.

Неявна різницева схема

Для апроксимації диференціального рівняння використовуємо нижній шаблон рисунка 5.1. Вважаємо, що вираз (5.2) апроксимує похідну за часом у момент . Тоді й другу похідну по х в (5.3) беремо в тому ж часовому шарі. Схема відповідає нижньому шаблону, наведеному на рисунку 5.1.

Маємо таку апроксимацію диференціального рівняння:

(5.7)

Вираз (5.7) записується для кожного внутрішнього вузла стрижня і являє собою систему лінійних алгебраїчних рівнянь для невідомих значень температури в момент часу . Кожне рівняння системи містить три невідомі, тобто система є трьохдіагональною, і може бути вирішена методом прогону або LU - розкладання. Додамо рівняння, що відображають граничні умови й запишемо систему в канонічному виді (позначимо ):

Домінування діагональних елементів системи забезпечує стійкість методу прогону для її розв’язку.

Така схема стійка без будь яких додаткових обмежень на узгодження кроків за часом і координаті x. Розмір кроків визначається тільки вимогами точності отриманого розв'язку.

Оскільки має місце апроксимація крайової задачі порядку й стійкість різницевої схеми, то буде мати місце збіжність розв'язку різницевої задачі до розв'язку даної диференціальної задачі:

Порядок апроксимації по х крайової задачі зменшився за рахунок апроксимації похідної у крайовій умові . Якщо на правому кінці задаються значення функції, то порядок апроксимації крайової задачі буде .

Апостеріорна оцінка погрішності може бути отримана за правилом Рунге на основі подвійного прорахунку за формулою (2.7). Використовуючи екстраполяцію Ричардсона можна уточнити отриманий розв'язок, підвищивши точність розв'язку на порядок (формула (2.8)).