Рішення систем диференціальних рівнянь в MATLAB
Для вирішення систем диференціальних рівнянь існує кілька вбудованих процедур. Розглянемо застосування процедури ode45. Як один з можливих форматів виклику, можна запропонувати такий: [t, r] = ode45 (@DiffEquationFunction, [Tstart, Tfinish], StartVector). Відзначимо наступне, що процедура ode45 може вирішити систему рівнянь такого вигляду:, де - є векторна функція.
Розглянемо приклад, який ілюструє створення вихідної функції DiffEquationFunction для виклику її процедурою ode45. Нехай деяка точка маси рухається в гравітаційному полі нерухомих точок з масами і (див. Рис.). Рівняння сил в гравітаційному полі точок і буде наступним:.
Як бачимо, дане диференціальне рівняння має другий порядок. Але його можна звести до системи диференціальних рівнянь першого порядку:
Будемо вважати, що дана задача є "плоскою" і введемо такі позначення:,, і. Виходячи з таких позначень, систему диференціальних рівнянь руху точки в гравітаційному полі можна представити таким чином:.
Без шкоди для змісту рішення, значення гравітаційної постійної приймемо рівної 1,,,,
У такому вигляді систему рівнянь можна вже записати як файл-функцію, що ми і зробили, назвавши її threepoint (t, x).
M1 = 50; M2 = 0; C1x = 5; C1y = 0; C2x = 0; C2y = 10;
І вирішимо систему диференціальних рівнянь, викликавши процедуру ode45 з файлу-функції dynpoint.m.
x1 = 5; y1 = 0; x2 = 0; y2 = 100;
При таких початкових параметрах () наша точка рухається в поле одного об'єкта.
Те, що ми отримали цілком можна назвати прийнятним результатом. Судити про коректність результату з точки зору фізики надамо фізикам. Нашою метою було коректна підготовка системи диференціальних рівнянь для її вирішення процедурою ode45. З іншого боку, застосовуючи інші вбудовані процедури можна отримати отлічающеся в корені результати.
Спробуємо трохи поекспериментувати і введемо значення. Це повинно внести обурення в орбіту рухається точки.
Хрестиком і зірочкою на графіку відзначені відповідно і. Крім процедури ode45 існує ще ряд вбудованих процедур вирішення систем диференціальних рівнянь. Їх опис можна знайти в книгах, зазначених у списку літератури. Оскільки, рішення систем диференційних рівнянь є дуже важливим моментом, то ми не обмежимося одним прикладом.
Розглянемо наступну задачу, знову ж з фізики. Розглянемо траєкторію руху кулі під дією сили тяжіння. При відсутності опору повітря, як відомо з курсу середньої школи, це буде парабола. Ми ж розглянемо випадок, коли сила опору повітря пропорційна квадрату швидкості і протилежна напрямку руху. Як каже шкільний курс фізики, рівняння балансу сил буде наступним:
З огляду на те, що прискорення - похідна швидкості за часом, розпишемо це рівняння в векторному вигляді:
Де - щільність повітря, - маса кулі, - площа поперечного перерізу. Розпишемо це рівняння покоординатно:
У такому вигляді система диференціальних рівнянь готова для того, щоб спробувати вирішити її за допомогою процедури ode45. Нехай маса кулі - 10 грам, діаметр - 1 см, щільність повітря - 1 кілограм на кубічний метр. З такими даними ми складемо файл-функцію airpoint.m.
g = 10; ro = 1; s = 0.0001; m = 0.01; k = ro * s / 2 / m;
У такій системі рівнянь аргументом є швидкість і час. Якщо ми хочемо знайти траєкторію руху, то ми повинні взяти до уваги:
Отримані масиви точок,, можна в подальшому обробити. Провести інтерполяцію, апроксимацію і т.д. Але ми інтегрування замінимо підсумовуванням:
В цьому випадку початковий момент часу дорівнює 0, куля знаходиться на початку координат. Створимо файл-функцію dynairpoint, в якій викличемо процедуру ode45. У початковий момент часу швидкість кулі по горизонталі 800, по вертикалі - 100.