Міністерство освіти і науки України
Сумський державний університет
Кафедра інформатики
Курсова робота
на тему:
«Дослідження збіжності рішень для диференціальних рівнянь участинних похідних, отриманих методом сіток»
Суми 2006
Вступ
Актуальність теми. Задачі, що відносяться до диференціальних рівняньу частинних похідних другого порядку, виникають у різних прикладних областях,зокрема в задачах акустики, електродинаміки, динамічної теорії пружності тощо.Найчастіше при їх розв’язуванні використовують метод скінчених різниць(методсіток). За цим методом вихідну диференціальну задачу у частинних похіднихзамінюють відповідною різницевою схемою, що є системою скінченої кількостіалгебраїчних рівнянь. Тобто для розв’язування неперервної задачі будують дискретнумодель, характер поведінки якої описують різницеві рівняння. Очевидно, будь-якадискретна модель не тотожна вихідній неперервній задачі.
Особливістю наближених методів є те, що кожному рівнянню можнапоставити у відповідність велику кількість різницевих апроксимацій, що маютьмайже однакові характеристики. Тому побудова різницевих схем, властивості якихякнайповніше відповідають вихідній диференціальній задачі, — суть і предметметоду скінчених різниць, а розвиток теорії різницевих схем природно шукають упокращенні порядку апроксимації, а також у зменшенні кількості арифметичнихоперацій для знаходження розв’язків. Іншими словами, різницева схема повиннаякомога краще моделювати властивості вихідного диференціального рівняння, дотого ж кількість арифметичних дій, потрібних для знаходження розв’язку, маєбути по можливості пропорційна кількості вузлів сітки.
Побудова різницевих схем для рівнянь у частинних похідних зузагальненими розв’язками, швидкість збіжності яких узгоджена з гладкістю цихрозв’язків, привертає сьогодні особливу теоретичну увагу. Як зазначається абоприймається за очевидне у кожній роботі з чисельних методів, основним питаннямдля теорії та практики наближених методів є питання точності розв’язку.Дослідження задач з негладкими розв’язками для рівнянь гіперболічного типупотребують особливої уваги через те, що негладкості середовища для такихрівнянь не зникають з часом. Проблема узагальнюється таким чином: як покращититочність наближеного методу, не збільшуючи при цьому паразитичних осциляцій,які з’являються при переході на кожний наступний ярус. Це явище виникає, колирозв’язок негладкий, має розриви та особливі точки (наявні сконцентрованізовнішні сили, точкові джерела тощо). Причина таких осциляцій — дисперсіярізницевої схеми по відношенню до диференціальної задачі, тобто відмінність(відставання або випередження) фазової швидкості сіткових гармонік від гармонікдиференціальних. Звідси ясно, якою важливою є побудова таких схем длярозв’язування гіперболічних рівнянь, де враховані дисперсійні властивостінеперервної моделі і, можливо, до мінімуму зведений спотворюючий вплив цихвластивостей.
Стан проблеми. Огляд літератури. Проблема існування дисперсіїрозглядалася багатьма авторами. Першими роботами, у яких було відмічено зв’язокміж дисперсією різницевих схем та втратою точності розв’язків рівняньгіперболічного типу, були роботи К. Роберта, В. Вайса та Дж. Фромма, в якихдисперсія досліджувалася для різних різницевих схем, що апроксимуютьгіперболічне рівняння першого порядку. Надалі на існування дисперсії дляодновимірних рівнянь вказувалося в роботах С. Орзаґа, Р. Чина. М.М. Москальковимпоказаний зв’язок осциляцій сіткових розв’язків з дисперсією гармонікрізницевої схеми для одновимірних гіперболічних рівнянь як першого, так ідругого порядку. Різними авторами проводився дисперсійний аналіз для різнихзадач математичної фізики. Для рівнянь газової динаміки розв’язувалися проблемизв’язку дисперсії та стійкості різницевих схем. Для спектральних методіврозв’язування задач гідродинаміки подібні проблеми ставилися. Огляд робіт,присвячених цьому питанню, можна знайти в роботі Л. Трефетхена, дедосліджувалися дисперсійні властивості різницевих схем, що апроксимуютьдвовимірне хвильове рівняння. За допомогою методу диференціальних наближеньпроблема дисперсії також досліджувалася.
Пропонувалися різні методи боротьби з наслідками дисперсії —паразитичнимиосциляціями розв’язків гіперболічних рівнянь. Один з них — введення в рівняннятак званої штучної в’язкості (див. роботу П.Роуча). Однак цей спосіб не завждизадовільний: втрачається справжній профіль розв’язку та ускладнюютьсяалгоритми, особливо, коли необхідно працювати з великою кількістю вузлів сітки.За іншими методами, що враховують існування дисперсії, пропонується вводити додатковийантидисперсійний ярус, або розглядати не ортогональні сітки як на площині, такі в тривимірному просторі.
В роботах О.С. Макаренка та М.М. Москалькова було вперше доведено,що в двовимірному випадку, виявляється, існує залежність дисперсії різницевоїсхеми не лише від номера сіткових гармонік, а й від напрямку руху хвилі. Длянеявних різницевих схем для двовимірного гіперболічного рівняння було показано,що можна суттєво покращити дисперсію у заданому напрямку руху хвилі задопомогою вибору вагів схеми.
В роботі В.Л. Макарова, С.В. Макарова, М.М. Москалькова булопомічено, що у різницевих схем на правильних трикутних сітках дисперсійнівластивості кращі за дисперсійні властивості схем на звичайних шаблонахпрямокутної сітки. Пояснюється це тим, що існує залежність між виглядом шаблонута дисперсією схеми. Відмічалося, також, що на правильній трикутній сітці можнапобудувати схеми четвертого порядку точності для рівняння Пуассона. Впершеподібне покращення апроксимації для не ортогональних сіток розглядав В.І. Лебедев.
Чисельні методи рішення диференціальних рівнянь у частиннихпохідних 2-го порядку
Загальний вигляд диференціальних рівнянь у часткових похідних 2-гопорядку:
/> (1)
Розв’язком рівняння (1) називається функція u=u(x,y), щоперетворює це рівняння на тотожність. Графіком розв’язку є поверхня в просторіOxyu.
Рівняння (1) називається лінійним, якщо воно першого степеня щодошуканої функції і всіх її похідних і не містить їхніх добутків, тобто церівняння може бути записане у вигляді
/> (2)
У загальному випадку A,B,C,a,b,c — це коефіцієнти, що можутьзалежати тільки від х, у. Уводиться визначення дискримінанта:
/>
У залежності від знака дискримінанта D лінійне диференціальнерівняння (2) відноситься до одного з наступних типів:
Якщо D>0, то рівняння еліптичне.
Якщо D=0, то рівняння параболічне.
Якщо D
Якщо D не зберігає знак – рівняння змішаного типу.
Приклади:
/> рівняння еліптичного типу.
/> - рівняння параболічного типу.
/> рівняння гіперболічного типу.
Рівняння вигляду /> — називається характеристичнимрівнянням. Розв’язки цього рівняння називають характеристиками.
Початкові і крайові умови
Диференціальні рівняння в часткових похідних мають незліченумножину розв’язків, тому для однозначності розв’язку необхідно до вихідногорівняння приєднати додаткові умови. Для диференціальних рівнянь у частковихпохідних 2-го порядку вони можуть бути початковими і граничними. По суті,розрізнити ці умови можна лише в тому випадку, коли одна з незалежних зміннихдиференціального рівняння відіграє роль часу, а інша – роль координати. Тодіякщо умови задані для початкового моменту часу, то це початкові умови, а умови,що відносяться до фіксованих значень координат – граничні або крайові.
Представлення часткових похідних у скінчено-різницевому вигляді
Розглянемо диференціальне рівняння />;
Апроксимуємо часткові похідні відповідними різницями:
/>
/>
/>
/>.
Аналогічно можна записати
/>
/>
тобто різницю зміщаємо до центра
/>
Формула для змішаної похідної
/>
/>
/>
Ці формули переходу до різницевих схем можна записати,використовуючи позначення:
/>,
/>,
/>.
Скористаємося розкладанням у ряд Тейлора:
/>
Для того, щоб оцінити похибку 2-ої похідної заміняємо /> в першійзаміні, у другій заміні />.
Отже, для />
/>
а для />
/>
Користуючись цими розкладаннями, одержимо
/>
Аналогічну формулу можна записати для похідної по у.
/>.
Аналогічно можна одержати оцінки для інших похідних.
Метод сіток
Ідея методу сіток відома давно ще за часів Ейлера. Однак,практичне використання цього методу наштовхувалося на серйозні труднощі, томущо одержання з його допомогою досить точного розв’язку крайової задачі звичайноприводило до колосальних систем алгебраїчних рівнянь, на розв’язок яких приручному розрахунку були потрібні роки. Положення різко змінилося з появоюелектронних обчислювальних машин. Метод сіток допускає зручну реалізацію наЕОМ, тому що застосування його зазвичай зводиться до масової повторюваностіоднорідних циклів. В даний час метод сіток є одним з найбільш ефективнихметодів розв’язку лінійних диференціальних рівнянь. Метод сіток (інакше методскінчених різниць) для наближеного розв’язку крайових задач двовимірнихдиференціальних рівнянь полягає в наступному:
У плоскій області G, у якій розшукується розв’язок, будуєтьсясіткова область Gh, що складається з однакових осередків і наближає дануобласть G;
Задане диференціальне рівняння заміняється у вузлах побудованоїсітки відповідним скінчено-різницевим рівнянням;
На підставі граничних умов установлюються значення шуканогорозв’язку в граничних вузлах області Gh
Розв’язавши отриману систему скінчено-різницевих рівнянь, мизнайдемо значення шуканої функції у вузлах сітки, тобто будемо мати чисельнийрозв’язок нашої задачі. Вибір сіткової області здійснюється в залежності відконкретної задачі, але у всіх випадках контур Гh сіткової області Gh вартовибирати так, щоб він якнайкраще апроксимував контур Г заданої області G. Сіткабудується таким чином, щоб вузли (xi,yi) сітки Sh або належали області G, абовідстояли від її границі Г на відстань меншому, ніж h. Точки (вузли) сітки Shназиваються сусідніми, якщо вони розміщені один від одного в напрямку осі Охабо осі Оу на відстань, що дорівнює кроку сітки h. Вузол Ah сітки Shназивається внутрішнім, якщо він належить області G, а всі чотири сусідніх зним вузла – множині Sh; інакше він називається граничним. Граничний вузол сіткиSh називається вузлом першого роду, якщо він має сусідній внутрішній вузол цієїсітки, інакше граничний вузол називається вузлом другого роду. Внутрішні вузлиі граничні вузли першого роду сітки Sh називаються розрахунковими точками.Граничні вузли другого роду не входять в обчислення і можуть бути вилучені ізсітки.
На перший погляд процедура застосування методу сіток, щоскладається з трьох етапів, може здатися простою і легко реалізованою. Однакнасправді це не так. Через велику розмаїтість типів і розмірів сіток, видіврівнянь у часткових похідних, граничних і початкових умов, можливихкінцево-різницевих апроксимацій цих рівнянь і методів їхнього розв’язку, чисельнерозв’язку рівнянь у часткових похідних вимагає модифікацій алгоритму прирозгляді кожного конкретного приклада.
Стійкість скінчено-різницевої схеми для розв’язку рівняньпараболічного типу (рівняння теплопровідності)
Як приклад рівняння параболічного типу розглянемо рівняннятеплопровідності:
/>,
де u=u(x,t) – температура, t – час, /> - довжина стрижня.
Для простоти покладемо, а=1.
Початкові і крайові умови:
/>,
/> и./>/>.
При використанні скінчено-різницевої схеми для розв’язку крайовоїзадачі виникає питання про стійкість такої схеми. Під цим розуміють наступне:скінчено-різницева схема називається стійкою, якщо малі похибки в процесірозв’язку загасають або у всякому разі залишаються малими при необмеженомузбільшенні номера поточного шару.
Для рівняння
/> (3)
скінчено-різницева схема матиме вигляд:
/>. (4)
З'ясуємо умови стійкості з граничними і початковими умовами
/> (5)
Маємо:
/> і />,
де />, />.
Переходячи до скінчених різниць у рівнянні (4), будемо мати:
/>=0 (6)
У граничних вузлах сітки />Г виконані такі умови:
/>, />, />.
Припустимо, що в точках початкового шару t=0 допущена помилка />, тобто
/>,
і нехай /> - розв’язок рівняння (6):
/>. (7)
яке задовольняє граничним умовам, що містять помилку:
/>, />, />.
Нас цікавить, як зміниться похибка /> при необмеженому зростанні номераj. Віднімаючи з рівняння (7) рівняння (6), для похибки /> одержимо скінчено-різницеверівняння.
/>=0. (8)
На границі Г області маємо:
/> (8а)
Частковий розв’язок рівняння (8) будемо шукати у вигляді
/>, (9)
де числа /> і p (р>0) підберемо так, щобвираз (9) задовольняв рівнянню (8) і однорідним крайовим умовам.
/>.
Користуючись ними маємо:
/>/>,
звідки випливає, що pl=m/>і /> (m=1,2,3……).
Отже,
/>.
Підставляючи цей вираз в рівняння (8), будемо мати:
/> (10)
Після перетворень рівняння (10) набуде вигляду:
/>.
Звідси
/>. (11)
Зауважимо, що />не залежить від точки (/>). Таким чином,для однорідного рівняння (8) одержуємо лінійно незалежні розв’язки вигляду:
/> (m=1,2,…....,n-1),
причому кожен розв’язок задовольняє однорідним крайовим умовам.
Лінійна комбінація цих розв’язків
/> (12)
також є розв’язком рівняння (8), що задовольняє при будь-якихзначеннях коефіцієнтів /> однорідним крайовим умовам. Цікоефіцієнти підбираються так, щоб виконувалася перша умова (8а), тобто щоб /> (і=1,2,…...,n-1).
Для стійкості розглянутої скінчено-різницевої схеми (6) необхідно,щоб при будь-яких значеннях постійних /> функція />, обумовлена рівністю (12),залишалася обмеженою при />.
Для цього досить, щоб для всіх m була виконана рівність:
/>.
Звідси
/> і /> (m=1,2,…...,n-1).
Остання нерівність буде виконана, якщо виконується умова:
/>. (13)
Отримані нерівності дають достатні умови стійкості розглянутоїскінчено-різницевої схеми для змішаної задачі у випадку рівняннятеплопровідності (параболічного типу).
Гіперболічний тип
Як приклад рівняння гіперболічного типу розглянемо рівнянняколивань однорідної обмеженої струни:
/>,
де u=u(x,t) – зсув струни, t – час, /> - координата довільної точкиструни.
Для простоти покладемо, а=1.
Початкові і крайові умови:
/>.
/> и./>/>
Схема стійка, якщо виконано умову Куранта k
Недолік схеми в тім, що як тільки обрана величина кроку сітки h унапрямку x, з'являються обмеження на величину кроку /> за змінною t. Якщо необхіднозробити обчислення для великого значення величини T, то може знадобитися великакількість кроків по змінній t. Зазначений недолік характерний для всіх явнихрізницевих схем.
Еліптичний тип
Як приклад рівняння еліптичного типу розглянемо задачу Дирихле(перша крайова задача для рівняння Лапласа />):
/>,
Крайова умова:
на колі /> (Г) виконується />.
Перепишемо систему рівнянь у вигляді, зручному для застосуванняметоду простої ітерації:
для внутрішніх вузлів
ui,k= -/>(14)
для граничных вузлів
/> (15)
Тут для внутрішніх вузлів використовувався п’ятиточковий шаблон,зображений. Припустимо, що gi,k
для внутрішніх вузлів
/>
для граничних вузлів
/>
р=0,1,2,…,/>задане.
Доведено, що якщо gi,k збігаються доточного розв’язку різницевої схеми ui,k або системи рівнянь (14), (15) і маємісце оцінка
max />
i,k i,k
де q=max />.
i,k
Доведення цього твердження полягає в перевірці умови збіжності методу простої ітерації для системи лінійних рівнянь, при цьому мається на увазі, що невідомий вектор />утворює елементи ui,k. Наприклад, компоненти вектора />можна перенумерувати таким чином: нехай /> тоді
x1= u1.1, />x 2 =u 2.1,…, x N1 = u N1.1;
x N1+1=u 2.1 x N1+2 =u 2.2,…, x 2N1 =u N1.2;
…………………………x N1N2 =u N1N2.
Відносно вектора />= /> різницева схема є системоюлінійних рівнянь в матричному записі /> де матриця А має в кожному рядкуне більше п’яти елементів
/>/> /> />…
/> /> />…
А=… />… /> /> />… />…
… />… /> /> />… />...
… /> />
Це пов’язано з тим, що похідні в кожному внутрішньому вузлі (i,k)апроксимувались за п’ятьма сусідніми вузлами.
/> Розв’язання різницевихрівнянь при h 0 збігається до точного розв’язання крайової задачі зі швидкістю,яка визначається порядком апроксимації рівнянь та крайових умов. Таким чином,для точного розв’язання (u(x,y)/>) оцінки похибки
/> max />O(h2), h 0 (16)
i, k
Оцінка похибки (16) є справедливою, якщо точний розв’язокнеперервно диференційований чотири рази в області G. Для областей з кутовимиточками, наприклад прямокутника, взагалі кажучи, u(x,y) />. Але якщо граничнафункція, тобто /> задовольняє в кутах спеціальніумови узгодження, то точний розв’язок u(x,y) /> і є вірною оцінка (16).
Для прямокутної області G=/>такими умовами узгодження можутьбути:
достатня гладкість />;
функція />повинна задовольняти в кутахпрямокутника диференціальне рівняння.
Оцінка похибки (8.96) має в основному теоретичне значення,оскільки містить константу С, яку практично важко визначити
/>
max />ch2+ O(h2), h 0
i, k
Тому в реальних розрахунках використовується правило Рунге оцінкипохибки, аналогічне тому, яке використовується в чисельному розв’язанні задачіКоші і розв’язанні звичайних диференціальних рівнянь. Робиться два варіантирозрахунку /> зкроком h та />;тоді похибка має вигляд
max />/> max />+О(h2)
і головна частина похибки визначається на вузлах, що збігаються.
Потрібно зазначити, що рівномірними прямокутними сітками найбільшзручно користуватиcя при розв’язанні задач у прямокутних областях. Якщо областьмає форму паралелограма(скошена система), то користуються координатами, осі якихпаралельні сторонам цього паралелограма. Декартові прямокутні координатипов’язані з косокутними координатами /> співвідношеннями />, де а — кут між />. Удиференціальних виразах похідні за х та у замінюються похідними за />. Усі похідніапроксимуються за допомогою центральних різниць. Якщо область має форму кола,зручно користуватись полярними координатами />
Наведемо деякі загальні зауваження. При чисельному розв’язаннікрайових задач для диференціальних рівнянь в частинних похідних методом сітокможуть бути використані тільки різницеві схеми, які збігаються, оскільки вцьому разі можна розраховувати на отримання наближеного розв’язку задачі,достатньо близького до точного. Але й різницеві схеми, що збігаються не завждиможуть бути використані при практичному розв’язанні задачі, оскільки привикористанні методу сіток при обчисленні значень граничних функцій та правоїчастини виникають похибки. Щоб ці похибки не спотворили істинного розв’язаннярізницевої схеми, остання повинна бути стійкою за граничними умовами і заправою частиною. При використанні нестійкої різницевої схеми спотворення істинногорозв’язку тим сильніше, чим дрібніша сітка; при використанні ж великої сітки неможна розраховувати на те, що розв’язок різницевої схеми буде близький доточного розв’язку крайової задачі для диференціального рівняння в силу поганоїрізницевої апроксимації рівняння.
Крім того, під час розв’язання різницевої задачі в процесірозрахунків нам обов’язково доведеться округляти значення розв’язків у вузлахсітки. Ці помилки можуть значно спотворити картину розв’язання, тому необхідноювимогою є стійкість різницевої схеми що до помилок, які виникають в результатіокруглення значень розв’язку у вузлах сітки. Оскільки помилки округленнязначень розв’язку в вузлах сітки, принаймні, в найпростіших випадках можнакомпенсувати зміною правої частини різницевого рівняння, то особливо суттєвою євимога до стійкості правої частини. Необхідно взяти до уваги й числовийалгоритм, який використовується для розв’язання різницевої схеми. Навіть у випадку,коли різницева схема стійка за граничними умовами і за правою частиною, приневдалому виборі алгоритму для розрахунку розв’язання цієї різницевої схемиможе відбутися сильне накопичення обчислювальної похибки, у цьому разінестійким буде сам процес розрахунку. Нестійкі алгоритми розрахунку практичнонепридатні у випадку дрібної сітки.
Вибір оптимального кроку
Припустимо, що межа абсолютної погрішності при обчисленні функції /> в кожній точцізадовольняє нерівність
/> (17)
Хай в деякій околиці крапки /> похідні, через які виражаютьсязалишкові члени, безперервні і задовольняють нерівностям
/> (18)
де /> - деякі числа. Тоді повнапогрішність (без урахування погрішностей округлення) не перевершує відповідновеличин
/>
Мінімізація по /> цих величин приводить донаступних значень:
/> (19)
при цьому
/> (19а)
Якщо при вибраному для якої-небудь значенні /> відрізок не виходить замежі околиці точки, в якій виконується відповідна нерівність (17), то знайдене єоптимальним і повна погрішність чисельного диференціювання оцінюєтьсявідповідною величиною (19).
Дослідження точності
Дослідження точності одержаних виразів при чисельних розрахункахзручно робити за допомогою апостеріорної оцінки, по швидкості убивання членіввідповідного ряду Тейлора. Якщо крок сітки достатньо малий, то погрішністьблизька до першого відкинутого члена.
Таким чином порядок точності результату по відношенню до крокусітку рівний числу залишених в ній членів, або іншими словами, він рівний числувузлів інтерполяції мінус порядок похідної. Тому мінімальне число вузлівнеобхідне для обчислення m-ой похідної, рівне m+1; воно забезпечує першийпорядок точності.
Ці висновки відповідають принципу: при почленном диференціюванніряду швидкість його збіжності зменшується.
Якщо врахувати погіршення збіжності ряду при диференціюванні, томожна зробити висновок: навіть якщо функція задана добре складеною таблицею надосить докладній сітці, то практично чисельним диференціюванням можна визначитипершу і другу похідні, а третиною і четвертую – лише задовільно. Більш високіпохідні рідко вдається обчислити із задовільною точністю.
Структура похибки розв'язку задачі
Побудувавши математичну модель, намагаються знайти її розв'язок.Для складних прикладних задач, як правило, не існує точного розв'язку у виглядіявних формул або скінченої послідовності арифметичних операцій, кожна з якихвиконується точно. Тоді вдаються до чисельних методів — могутнього математичногозасобу розв'язування задач. Найпростіші чисельні методи виникли і широковикористовувалися задовго до появи ЕОМ. Але є багато прикладних задач, для якихзнайти розв'язок без застосування ЕОМ практично неможливо. Сучасні швидкодіючіЕОМ стали стимулом для розробки нових чисельних методів.
Застосовувати чисельні методи для розв'язування прикладних задачна базі ЕОМ треба обережно, оскільки точність знайденого розв'язку залежить відбагатьох факторів. При цьому слід уміти оцінити похибку обчисленого розв'язку.
Похибка розв'язку задачі складається з похибки математичноїмоделі, неусувної похибки, похибки методу і обчислювальної похибки.
Похибка математичної моделі пов'язана з тим, що модель описуєявище наближено, з припущеннями і спрощеннями. Тому треба мати уявлення проточність кінцевого результату, щоб спростити побудову математичної моделі.
Неусувна похибка зумовлена похибками у вхідних даних задачі. Воназалежить від методу розв'язування задачі. Але, щоб правильно обрати метод івизначити точність обчислень, важливо знати межі неусувної похибки.
Похибка методу пов'язана з необхідністю заміни неперервної моделідискретною або з обривом нескінченного ітераційного процесу після скінченоїкількості ітерацій.
Знайти чисельно розв'язок у(х) в усіх точках відрізка [а;b]неможливо, оскільки їх безліч. Тому на відрізку [а;b] беруть скінчену кількістьточок (хі = a+ih, i = 0,1… п, а + пh
в них знаходять значення у(х). Початкове значення y(a)=y0 намвідоме. Для знаходження інших значень уi= у(хі) (i = 1,2,.., n) диференціальнерівняння розглядають не на всьому відрізку [a;b], а тільки в зазначених точках y’(xi)=f(xi,y(xi)).
Замінивши похідну у'(х) її наближеним значенням (уi+1 — уі)/h,дістають систему рівнянь
уi+1 — уі = пf(хі, уі), i = 0,1,2,...,n-1. (20)
Звідси послідовно знаходять y1, y2,… yn.
Якщо в рівняння (20) замість уi i уi+1 підставити точні значеннярозв'язку у(хi) і y(xi+i)> то рівності задовольняться лише наближено.
Похибку, яку дістають від заміни неперервної моделі дискретною,називають похибкою дискретизації (або похибкою апроксимації),
Крім похибки дискретизації, існує інший тип похибки чисельнихметодів. В основі багатьох методів лежить ідея ітераційного процесу, в ходіякого будується за певним правилом послідовність наближень до розв'язку задачі.Якщо ця послідовність має границю, коли кількість членів послідовності прямуєдо нескінченності, тоді ця границя буде розв'язком даної задачі. Але на ЕОМможна обчислити тільки скінчену кількість членів послідовності. Похибку,спричинену обривом ітераційного процесу, називають похибкою збіжності. Похибкуметоду намагаються звести до величини, яка в кілька разів менша від похибкивхідних даних.
Отже, похибку чисельного методу можуть утворювати похибкидискретизації або похибки збіжності, або ж для деяких методів обидва типипохибок одночасно. Всі ці похибки, а також методи їх аналізу і регулюваннярозглядаються при побудові конкретних чисельних методів.
Обчислювальні похибки пов'язані з похибками округлення чисел. Обчислення,як ручні, так і на ЕОМ, виконують з певною кількістю значущих цифр. Це вноситьу, результат похибку округлення, яка нагромаджується в ході обчислень. Похибкиокруглення можуть по-різному впливати на кінцевий результат. У результатівиконання мільйонів операцій, кожна з яких вносить невелику похибку, сумарнапохибка округлень може значно перевищити шуканий результат обчислень. Але вокремих операціях похибки округлень можуть мати різні знаки і частковокомпенсувати одна одну. Тому, якщо немає систематичних причин, випадковенагромадження похибок округлення незначне.
Систематичною причиною нагромадження похибок є, наприклад,віднімання близьких за величиною чисел, оскільки при малій абсолютній похибцічисел х1 і х2 відносна похибка (∆x1+∆x2)/|x1-x2| результату можестати великою.
Обчислювальні похибки виникають і під час перетворення чисел зоднієї системи числення в іншу, якщо основа однієї системи числення не єстепенем основи іншої. Це може призвести до того, що в новій системі численнячисло стане ірраціональним.
Втрата точності може статися і при додаванні до великого числадуже малих чисел. Для зменшення похибки додавати числа варто в порядку їхзростання. У машинній арифметиці комутативний і дистрибутивний закони алгебрине завжди виконуються. Обчислювальний алгоритм треба будувати так, щоб похибкаокруглень була значно меншою від усіх інших похибок.
Поняття стійкості та коректності
Похибки у вхідних даних задачі — неусувні. Обчислювач не може їхзменшити, але мусить знати, як вони впливають на точність кінцевого результату.Одні задачі мають похибку результату такого самого порядку, як і порядокпохибки вхідних даних, в інших задачах похибка результату може на кількапорядків перевищувати похибку вхідних даних. Чутливість задачі до неточностей увхідних даних характеризується поняттям стійкості.
Задача називається стійкою за вхідними даними, якщо її розвязок неперервнозалежить від вхідних даних, тобто малому приросту ∆х вхідної величинивідповідає малий приріст ∆у шуканого розв'язку. Іншими словами, малі похибкивхідних даних спричинюють малі похибки розв'язку задачі. Якщо ця умова невиконується, то задача вважається нестійкою за вхідними даними. Це означає, щонавіть незначні похибки вхідних даних можуть привести до як завгодно великихпохибок розв'язку, тобто розв'язок може бути зовсім спотворений. Томузастосовувати безпосередньо до таких задач чисельні методи не можна, оскількипохибки округлень при застосуванні методу будуть катастрофічно нагромаджуватисьу ході обчислень. Наведемо приклад нестійкої задачі, який належить Уілкінсону.
Введемо тепер поняття коректності задачі.
Задача називається коректно поставленою, якщо для будь-якихвхідних даних з деякого класу існує єдиний і стійкий за вхідними даними їїрозв'язок.
Для розв'язування некоректно поставлених задач застосовуватикласичні чисельні методи не варто, оскільки похибки округлень при розрахункахможуть катастрофічне зростати і призвести до результату, далекого від шуканогорозв'язку. Для розв'язування некоректно поставлених задач використовують такзвані методи регуляризацїї, які замінюють дану задачу коректно поставленою.
Програмна реалізація
Розглянемо приклад рівняння в частинних похідних гіперболічноготипу
/>
Область: />
Початкові умови:
/>
Граничні умови:
/>
Крок по осі Х: hx=0,1
Крок по осі Y: hy=0,05
Реалізуємо в Pascal
Const
n = 10;
m = 10;
hx = 1/n;
hy = 0.05;
l = hy/hx;
Var
u:Array [0..n,0..m] of Real;
i,j:Integer;
fo:Text;
Function f(x:Real):Real;
begin
f:=(1-x)*cos(Pi*x/2);
end;
Function g(x:Real):Real;
begin
g:=2*x+1;
end;
Function fi(y:Real):Real;
begin
fi:=2*y+1;
end;
Function psi(t:Real):Real;
begin
psi:=0;
end;
Begin
For i:=1 to n-1 do begin
u[i,0]:=f(i*hx);
u[i,1]:=f(i*hx)+hy*g(i*hx);
end;
For j:=0 to m do begin
u[0,j]:=fi(j*hy);
u[n,j]:=psi(j*hy);
end;
For j:=1 to m-1 do
For i:=1 to n-1 do
u[i,j+1]:=2*(1-sqr(l))*u[i,j]+sqr(l)*(u[i+1,j]+u[i-1,j])-u[i,j-1];
Assign(fo,'result.txt');
ReWrite(fo);
For j:=m downto 0 do begin
Write(fo,j*hy:4:2,' |');
For i:=0 to n do Write(fo,u[i,j]:8:4);
WriteLn(fo);
end;
For j:=1 to 94 do Write(fo,'-');
WriteLn(fo);
Write(fo,' y/x |');
For j:=0 to n do Write(fo,j*hx:8:4);
Close(fo);
End.
RESULT.txt
0.50 | 2.0000 1.8693 1.7537 1.6381 1.5294 1.4450 1.27240.9448 0.6245 0.3283 0.0000
0.45 | 1.9000 1.7551 1.6251 1.4873 1.3910 1.3473 1.24990.9784 0.6246 0.3280 0.0000
0.40 | 1.8000 1.6447 1.4945 1.3467 1.2658 1.2361 1.18390.9914 0.6389 0.3198 0.0000
0.35 | 1.7000 1.5355 1.3646 1.2229 1.1533 1.1193 1.08280.9644 0.6616 0.3114 0.0000
0.30 | 1.6000 1.4248 1.2419 1.1171 1.0498 1.0019 0.96120.8913 0.6724 0.3127 0.0000
0.25 | 1.5000 1.3121 1.1338 1.0257 0.9511 0.8863 0.83230.7809 0.6480 0.3258 0.0000
0.20 | 1.4000 1.2018 1.0433 0.9427 0.8549 0.7734 0.70410.6502 0.5763 0.3380 0.0000
0.15 | 1.3000 1.1015 0.9672 0.8628 0.7602 0.6636 0.57970.5145 0.4635 0.3252 0.0000
0.10 | 1.2000 1.0172 0.8986 0.7834 0.6670 0.5569 0.46000.3823 0.3289 0.2658 0.0000
0.05 | 1.1000 0.9489 0.8308 0.7037 0.5754 0.4536 0.34510.2562 0.1918 0.1556 0.0000
0.00 | 1.0000 0.8889 0.7608 0.6237 0.4854 0.3536 0.23510.1362 0.0618 0.0156 0.0000
y/x | 0.0000 0.1000 0.2000 0.3000 0.4000 0.5000 0.60000.7000 0.8000 0.9000 1.0000
Для дослідження збіжності розв’язку даної задачі, змінимо крок розбиття сітки. Наприклад, візьмемо крок,в двічі менший, як по осі X, так і осі Y і знайдемо розв’язок за допомогою цієїж програми.
Отримаємо:
0.50 | 2.0000 1.7347 1.5044 1.2505 0.6479 0.0000
0.40 | 1.8000 1.4766 1.2682 1.1311 0.6793 0.0000
0.30 | 1.6000 1.2472 1.0499 0.9330 0.6538 0.0000
0.20 | 1.4000 1.0568 0.8517 0.6944 0.5347 0.0000
0.10 | 1.2000 0.9008 0.6654 0.4551 0.3218 0.0000
0.00 | 1.0000 0.7608 0.4854 0.2351 0.0618 0.0000
-----------------------------------------------------------------
y/x | 0.0000 0.2000 0.4000 0.6000 0.8000 1.0000
Крок був взятий так, щоб можна було порівняти розв’язок в окремихвузлах, як в першому, так і в другому випадку. Методом порівняння можнапобачити, що розв’язок в вузлах сітки дещо відрізняється. З цього можна зробитивисновок, що на збіжність розв’язку також впливає такий фактор, як крокрозбиття сітки.
Виконаємо також обчислення в середовищі Excel і порівняємоотримані дані.
З отриманого результату, що наведений дещо нижче, можна зробитивисновок, що отримані дані співпадають, що свідчить о вірності розв’язкузадачі, хоча і мають різну точність, що пов’язана з точністю обчислення (вPascal точність залежить від типу змінної).
/>
/>
Узгодженість і збіжність
Різницева схема, що сходиться, повинна давати рішення, яке прагнеістинного рішення ДР при подрібненні сітки. Порушення збіжності може бутиобумовлено наступними причинами:
Неузгодженість PC з початковим ДР. При подрібненні сітки значенняшуканої функції при заданих незалежних змінних прагне певної межі. але ця межане співпадає з істинним рішенням.
Нестійкість PC. При подрібненні сітки значення шуканої функції призаданих незалежних змінних не прагне певної межі. У будь-якому випадку в ходіобчислювального експерименту нестійкість проявляє себе як швидке, катастрофічненаростання чисел.
Математичні основи питань збіжності добре розвинені тільки длялінійних ДР. Результати лінійної теорії використовуються у вигляді навіднихміркувань для нелінійних задач, а їх застосовність перевіряється потім в ходіобчислювального експерименту. Для лінійних ДУ в ЧП існує т.н, теорема Лакса, якаговорить, що за наявності стійкості апроксимація є необхідною і достатньоюумовою збіжності PC.
апроксимація + стійкість => збіжність PC
У вживанні до нелінійних ДР в ЧП використовування теореми Лаксаприпускає використовування принципу «заморожених» коефіцієнтів. Підметодом «заморожених» коефіцієнтів розуміють простий прийом, коликоефіцієнти в ДР замінюють деякими константами і аналізують одержане ДУ зпостійними коефіцієнтами. Якщо аналіз показує нестійкість PC для«заморожених» коефіцієнтів, то таку схему виключають з розгляду. Якщож схема для рівняння з «замороженими» коефіцієнтами буде злагодженоюі стійкою, тобто надія, що ця схема буде тією, що сходиться і для нелінійногоДР.
Так, щоб бути упевненим, що PC, що використовується, дає рішення,що сходиться при подрібненні сітки до істинного рішення початкового ДР,необхідно довести її узгодженість і стійкість.
Погрішність апроксимації
PC апроксимації ДР злагоджена з початковим ДР (або апроксимує ДР),якщо в межі, коли розміри осередків сітки прагнуть нуля, PC еквівалентна даномуДР в кожній з вузлових точок, тобто
PC і ДР при цьому записують, перекидаючи всі члени в одну сторону(щоб з другого боку залишився нуль).
Кількісною характеристикою узгодженості є погрішність (помилка)апроксимації ДР даної PC
Е=(РС-ДУ).
В термінах погрішності апроксимації: PC злагоджена з своїм ДР,якщо погрішність апроксимації у всіх вузлах сітки прагне нуля при будь-якомуподрібненні кроків сітки.
Провідні члени погрішності називаються порядком апроксимації ДРданої PC.
Як правило:
• для тих PC, що сходяться помилка чисельного рішення зменшуєтьсяподібно погрішності апроксимації;
• порядок апроксимації ДР (провідні члени погрішностіапроксимації) визначається порядком точності формул, використаних дляапроксимації похідних;
• для злагоджених явних PC можна встановити співвідношення міжкроками сітки, при виконанні якого досягається більш високий порядокапроксимації.
Висновок
Диференціальні рівняння в частинних похідних є широко вживанимматематичним апаратом при розробці моделей в самих різних областях науки ітехніки. На жаль, явне рішення цих рівнянь в аналітичному вигляді виявляєтьсяможливим тільки в окремих простих випадках, і, як результат, можливість аналізуматематичних моделей, побудованих на основі диференціальних рівнянь,забезпечується за допомогою наближених чисельних методів рішення. Об'ємвиконуваних при цьому обчислень звичайно є значним і використовуваннявисокопродуктивних обчислювальних систем є традиційним для даної областіобчислювальної математики. Проблематика чисельного рішення диференціальнихрівнянь в частинних похідних є областю інтенсивних досліджень.
Стосовно самої теми курсової роботи, треба відзначити, щозбіжність рішень ДРЧП досягається, в першу чергу через теорему Лакса –апроксимація + стійкість породжує збіжність. Також треба зазначити, що назбіжність також впливає, хоча і в меншій мірі, вибір кроку розбиття сітки, атакож різноманітні похибки, хоча вони і не значно впливають на розв’язок, протедещо спотворюють його. Тому, щоб досягти найбільш точного і оптимальногорішення потрібно враховувати всі фактори, що можуть впливати на збіжність іточність даного розв’язку.
Література
Ляшенко М.Я.,Головань М.С. Чисельні методи.-К.: Либідь, 1996.-288с.
www.software.unn.ac.ru/ccam/files/HTML_Version/index.html
www.arptek.ru/i18n
Лекції почисельним методам Л.Д. Назаренко.
www.ict.nsc.ru/rus/texbooks/akhmerov/matmodel/2-4-7.html
www.nsc.ru/rus/textbooks/akhmerov/ode/s-33/s-33.html