2.Вирішення звичайних диференційних рівнянь в системі MathCad

Є два типи задач, які можливо вирішувати за допомогою Mathcad:

• задача Коші — для яких визначені початкові умови на шукані функції, тобто задані значення цих функцій в початковій точці інтервалу інтеграції рівняння;

• крайові задачі — для яких задані певні співвідношення відразу на обох межах інтервалу.

ЗДР першого порядку

Диференціальне рівняння першого порядку може за визначенням містити, окрім самої шуканої функції y(t), лише її першу похідну y'(t). У переважній більшості випадків диференціальне рівняння можна записати в стандартній формі (формі Коші):

у(t)=f(y(t),t),                      (6)

і лише з такою формою вміє працювати обчислювальний процесорMathcad. Правильна з математичної точки зору постановка відповідноїзадачіКоші дляЗДРпершого порядку повинна, окрім самого рівняння, містити одну початкову умову — значення функціїy(t0) в деякій точціt0. Потрібно явно визначити функцію(t) на інтервалі відt0 доtx.

Для чисельної інтеграції однієї ЗДРу користувачаMathcad(починаючи з версіїMathcad2000Pro) є вибір — або використовувати обчислювальний блокGiven/odesoive, або вбудовані функції Rkfixed, Rkadapt, Bulstoer, як в колишніх версіяхMathcad. Першийшлях кращийпереважно з міркувань наочності представленнязадачіі результатів, а другийдає користувачеві більше важеліввпливуна параметри чисельного методу.

Всі методи засновані на апроксимації диференціальних рівнянь різницевими аналогами. Залежно від конкретної форми апроксимації, виходять алгоритми різної точності і швидкодії. У Mathcad використаний найбільш популярний алгоритм Рунге-Кутта четвертого порядку. Він забезпечує малу похибку для широкого класу систем ЗДР за винятком жорстких систем.

Обчислювальний блок Given/Odesolve

Обчислювальний блок для вирішення одного ЗДР, що реалізовує чисельний метод Рунге-Кутта, складається з трьох частин:

• Given — ключове слово;

• ЗДР

• і початкову умову, записану за допомогою логічних операторів, причому початкова умова має бути у формі в(t1) = b; odesolve(t, t1) — вбудована функція для вирішення ОДУ відносно змінної t на інтервалі (t0,t1).

Припустимо, і навіть часто бажано, задання функції Odesolve (t, t1, step) з трьома параметрами, де step— внутрішній параметр чисельного методу, що визначає кількість кроків, в яких метод Рунге-Кутта розраховуватиме рішення диференціального рівняння. Чим більше step, тим з кращою точністю буде отриманий результат, але тим більше часу буде витрачено на його пошук. Підбором цього параметру можна у декілька разів прискорити розрахунки без істотного погіршення їх точності.

Приклад рішення задачі Коші для ЗДР першого порядку у' =у-у² за допомогою обчислювального блоку наведений в лістингу 1.

Лістінг 1. Вирішення задачі Коші для ЗДР первого порядку

Mathcad вимагає, щоб кінцева точка інтеграції ЗДР лежала правіше початкової: t0<t1 (у лістингу 11.1 t0=0,t1=10), інакше буде видано повідомлення про помилку. Результатом роботи блоку Given/odesoive є функція y(t), визначена на проміжку (t0,t1). Слід скористатися звичайними засобами Mathcad, щоб побудувати її графік або отримати значення функції в якій-небудь точці вказаного інтервалу, наприклад: y(3) =0.691.

Користувач має можливість вибирати між двома модифікаціями чисельного методу Рунге-Кутта. Для зміни методу необхідно натисненням правої кнопки миші на області функції odesolve викликати контекстне меню і вибрати в ньому один з двох пунктів: Fixed (Фіксований крок) або Adaptive (Адаптивний). За умовчанням застосовується перший з них, тобто метод Рунге - Кутта з фіксованим кроком.

Графік рішення даного рівняння показаний на мал. 1. Зверніть увагу, що він відповідає отриманню рішення в матричному вигляді (лістинг 1), тому по осях відкладені відповідні стовпці, виділені з матриці оператором <>.

Рис. 1. Рішення рівняння y' =y-y² (лістінг 1)

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

Системи ЗДР першого порядку

Для MathCad система диференціальних рівнянь має бути представлена в стандартній формі.

Задання системи еквівалентно наступному векторному представленню, де Y і В '' — відповідні невідомі векторні функції змінної t розміру N*1, а р — векторна функція того ж розміру і (N+i) кількості змінних (N компонент вектора і, можливо, t). Саме векторне представлення використовується для введення системи ЗДР в середовищі Mathcad.

Задача сформульована для систем ЗДР першого порядку. Якщо в систему входять рівнян-ня вищих порядків, її можна звести до системи більшого числа рівнянь першого порядку.

Вбудовані функції для рішення систем ЗДР

У Mathcad є три вбудовані функції, які дозволяють вирішувати задачу Коші для систем рівнянь різними чисельними методами.

• rkfixed(y0, t0, t1, M, D) — метод Рунге-Кутта з фіксованим кроком,

• Rkadapt(y0, t0, t1, M, D) — метод Рунге-Кутта зі змінним кроком;

• Buistoer(y0, t0, t1, M, D) — метод Булирша-Штера;

  • у0 — вектор початкових значень в точці to розміру NXI;

  • t0 — початкова точка розрахунку,

  • t1 — кінцева точка розрахунку,

  • M — число кроків, на яких чисельний метод знаходить рішення;

  • D — векторна функція розміру N*1 двох аргументів — скалярного t та векторного у . При цьому у — шукана векторна функція аргумента t того ж розміру N*1.

Кожна з наведених функцій видає рішення у вигляді матриці розміру (M+1)х(N+1). У її лівому стовпці знаходяться значення аргументу t, що ділять інтервал на рівномірні кроки, а в останніх N стовпцях — значення шуканих функцій y0(t),y1(t) ..,yn-1(t), розраховані для цих значень аргументу.Оскільки всього точок (окрім початкової)M, то рядківуматриці рішення буде всього M+1.

Лістінг 2. Рішення системи двох ЗДР за допомогою функції rkfixed.

По-перше, функція D, що входить до числа параметрів вбудованих функцій для рішення ЗДР, має бути функцією обов'язково двох аргументів. По-друге, другий її аргумент має бути вектором того ж розміру, що і сама функція D. По-третє, такий самий розмір має бути і у вектора початкових значень уо (він визначений у першому рядку лістингу правіше).

У другому рядку лістингу визначено число кроків, на яких розраховується рішення, а його останній рядок надає матричній змінній результат дії функції rkfixed. Рішення системи ЗДР буде здійснено на проміжку (0, 50).

Рис.2. Рішення системи з двох ЗДР (лістінг 2)за допомогою функції rkfixed.

Фазовий портрет типу змальованого на рис. 2 (останній графік), має одну стаціонарну точку (аттрактор), на яку "накручується" рішення. В теорії динамічних систем аттрактор такого типу називається фокусом.

Приклад вирішення ЗДР у MathCad за допомогою функції rkfixed

(це вектор з одного елемета)

Рис.3. Рішення одного ЗДР і побудова графіка за допомогою функції rkfixed.

Проглянути всі компоненти матриці і, які не вміщуються на екрані, можна за допомогою вертикальної смуги прокрутки. Щоб побудувати графік рішення, треба відкласти відповідні компоненти матриці рішення по координатних осях: значення аргументу — уздовж осіх, а і— уздовж осіy. Для виділення з матриці стовпчика слід скористатися на панелі матриць або функцієюsubmatrix. Як відомо, вирішення звичайних диференціальних рівнянь часто зручніше змальовувати не у такому вигляді, а у фазовому просторі, по кожній з осей якого відкладаються значення кожній із знайдених функцій. При цьому аргумент входить в них лише параметрично. Для його побудови потрібно було лише поміняти мітки осей на і, відповідно.