Вважаємо різницеві схеми в Mathcad Express

Продовжую замальовки про використання безкоштовного математичного редактора Mathcad Express. На цей раз пропоную звернутися до чисельного рішення диференціальних рівнянь (в сьогоднішньому прикладі з приватними похідними). Файл з подальшими розрахунками у форматах Mathcad і XPS ви знайдете тут.

У цій статті розглянемо, як можна порахувати в Mathcad Express диференціальне рівняння дифузії тепла за допомогою найпростішої явної різницевої схеми і зупинимося на властивості стійкості різницевих схем. Мова піде про те, як отримати ось таке рішення рівняння, що моделює остигання одномірного об'єкта:



На графіку показано початковий розподіл температури вздовж осі Х (червона лінія) і два розрахункові профілю – після першого кроку і після кількох кроків за часом.

Рівняння теплопровідності
В якості прикладу диференціального рівняння в приватних похідних ми розглядаємо рівняння теплопровідності (або, по-іншому, дифузії тепла).



Це рівняння параболічного типу, що містить першу похідну по часу t і другу по просторовій координаті x. Воно описує динаміку температури u(x,t) в присутності джерел тепла, наприклад, при нагріванні металевого стрижня. Таким чином, невідомою функцією, що підлягає визначенню, є функція u(x,t). Це рівняння лінійно, якщо джерело і коефіцієнт дифузії D або не залежать від температури, або джерело залежить від неї лінійно, оскільки в цьому випадку невідома функція u(x,t) (і всі її похідні) входить в рівняння першого ступеня.

Лінійне рівняння теплопровідності має аналітичне рішення, в той час, як переважна більшість нелінійних рівнянь в приватних похідних можна вирішити тільки чисельно. Для того, щоб рівняння в приватних похідних мала єдиний розв'язок, необхідно поставити потрібну кількість початкових і граничних умов, тобто, в даному випадку, співвідношення типу u(x,0)=Init(x) (початкова умова) разом з u(0,t)=f1(t) і u(L,t)=f2(t) (граничні умови, якщо рівняння вирішується на інтервалі від 0 до L).

В наведеній інтерпретації одномірне рівняння дифузії тепла описує динаміку температури деякого подовженого тіла, наприклад, металевого стрижня:



Ми будемо розглядати випадок нульового джерела тепла, а в якості початкової умови виберемо деякий реалістичний профіль температури (стрижень, сильно нагріте в центрі). В якості граничних умов задамо підтримання країв стрижня при постійній температурі. Зрозуміло, динаміка рішення описує очікувану картину вирівнювання температури стрижня з часом (за рахунок перетікання тепла від центру до периферії стрижня завдяки механізму теплопровідності). Швидкість передачі тепла визначається значенням коефіцієнта дифузії D. (Зацікавленому читачеві я пропоную поекспериментувати з доданими розрахунками в Mathcad Express, вибираючи різні значення коефіцієнта дифузії).

Різницева схема
Покажемо, як реалізовувати в Mathcad основний обчислювальний підхід до вирішення рівнянь в приватних похідних – метод сіток. Будемо робити це «вручну», щоб не вийти за рамки безкоштовної версії Mathcad Express.

Суть методу сіток полягає в покритті розрахункової області (x,t) сіткою з МхN точок, що визначить кроки по часу і простору τ і Δ відповідно. Тим самим, визначаються вузли, в яких буде здійснюватися пошук рішення. Потім треба замінити диференціальне рівняння апроксимуючим його рівняння в кінцевих різницях, для кожного (i,n)-го вузла сітки. У випадку рівняння теплопровідності досить просто замінити першу похідну по часу і другу по простору їх різницевими аналогами (такий метод дискретизації називають методом Ейлера, а отриману систему різницевих рівнянь – різницевої схемою).



Виписана схема є явною, оскільки в ній невідоме значення шуканої функції (n+1 кроці по часу) можна виразити через вже обчислені раніше значення функції на попередньому n-м кроці («шарі») по часу. Більш того, т. к. і на кожному шарі рівняння пов'язані між собою рекуррентно, тобто значення в кожному невідомому сайті на n-му шарі можна знайти, знаючи тільки обчислені раніше значення з попереднього шару. Це дозволяє організувати процес т. зв. «біжить рахунки», розв'язуючи відповідні рівняння на кожному шарі (виходячи з відомого рішення на попередньому шарі).

Наступний розрахунок реалізує описаний алгоритм явної різницевої схеми біжить рахунки в Mathcad Express:



Власне, обчислені значення матриці U і наведені на графіку, який був показаний в самому початку статті (до ката).

Стійкість
В укладення зупинимося на властивості стійкості різницевих схем. Параметр задачі, який характеризує відношення кроків різницевої схеми по простору і часу, називається коефіцієнтом (або «числом») Куранта:



Істотно, що попередні розрахунки за явною схемою проводилися при співвідношенні Куранта менше 1. А от якщо спробувати збільшити крок за часом в явній схемі, так, щоб співвідношення Куранта стало більше 1, то ми зіткнемося з катастрофою: рішення «розвалюється», і замість очікуваного (і прорахованого при менших співвідношеннях Куранта випадках) рішення виходить быстроосцилирующее, абсолютно неправдоподібне рішення.



Цей приклад демонструє важливе властивість різницевих схем, зване стійкістю. Виявляється, багато схеми (безумовно, або при певному співвідношенні параметрів, як наша схема) їм не володіють, і застосовувати такі схеми для вирішення завдань не можна.

Більш детальну інформацію про різницевих схемах та прикладах їх розрахунку в Mathcad читач може знайти в моїй книзі з обчислювальної фізики (різницевими схемами присвячені глави 5 і 6, книга у вільному доступі). Зокрема, пропоную спробувати вирішити різницевим методом «зворотне рівняння теплопровідності», яке виходить заміною змінної t на -t і яке (по ідеї) має описувати реконструкцію вихідного профілю температури стрижня за даними спостереження через деякий час.

І повторюся, що файл з розрахунками, наведеними в цій статті, ви знайдете тут, а дистрибутив Mathcad Express можна скачати за адресою.

Джерело: Хабрахабр

0 коментарів

Тільки зареєстровані та авторизовані користувачі можуть залишати коментарі.