Моделювання динамічних систем — це моделювання поведінки динамічної системи в будь-який довільний змінний момент часу.^[1][2] Модель, як правило, описується системою звичайних диференціальних рівнянь, аргументом яких є час. Така система відображає реальний об'єкт лише з деяким наближенням, яке може бути задовільним або незадовільним для певного дослідження.^[3][4][5]Тривіальним прикладом об'єкта, описаного за допомогою диференціальних рівнянь, може бути басейн, який заповнюється водою з труб.^[6]
Також при моделюванні динамічних систем використовують різницеві рівняння й системи рівнянь у випадку, коли зміна процесу відбувається стрибкоподібно, або дискретно. Такі динамічні процеси зустрічаються в онкології, динаміці популяцій, економіці, банківській справі.^[7]
Процес створення математичної моделі динамічної системи містить три основні частини:

• Емпірична;
• Теоретична;
• Математична.

В емпіричній частині зібрані дані, які були отримані зі спостережень та експериментів з ціллю дослідження об'єкта. Емпіричні закономірності та явища об'єднуються у теоретичній частині за допомогою розвитку основних концепцій. У математичній частині конструюються моделі для перевірки основних математичних концепцій. На цьому етапі відбувається процес обробки експериментальних даних, планування експериментів та спостережень.
Важлива перевага методів моделювання динамічних систем полягає в тому, що вони дозволяють різко скоротити обсяг і масштаби натурних експериментів.^[7]

В задачах моделювання важливим пунктом є аналіз властивостей і специфіки чисельної реалізації. Часто ми отримуємо представлення моделі, виходячи з її фізичних властивостей. Таке подання не завжди є зручним для чисельних експериментів. В моделюванні динамічних систем використовуються методи еквівалетного перетворення, а саме: перетворення диференціального рівняння ${\displaystyle n}$-го порядку до системи диференціальних рівнянь 1-го порядку, перетворення диференціальних моделей в інтегральні моделі, перетворення інтегральної моделі Вольтерра другого роду з ядром, що розділяється в диференціальну модель.
Методи еквівалентного перетворення диференціальних моделей в інтегральні моделі включають у себе:

• метод перетворення з розщепленням;
• метод послідовного інтегрування;
• метод старшої похідної.^[8]

Нехай подано звичайне диференціальне рівняння, що описує модель динамічного типу таким чином:
${\displaystyle D[y]=y^{(n)}(x)+\sum _{i=1}^{n}a_{i}y^{(n-1)}(x)=f(x),y^{(i)}(0)=C_{i},i={\overline {0,n-1,}}}$
або в операторній формі
${\displaystyle D[y]=f.}$
Для отримання ряду еквівалентних залежностей, що містять інтегральний оператор, застосовують прийом, що базується на різноманітних розщепленнях вихідного диференціального оператора. Дійсно, розщіплюючи оператор ${\displaystyle D}$ з використанням суми двох операторів ${\displaystyle D=D_{1}+D_{2},}$ отримуємо диференціальне рівняння
${\displaystyle D[y]=\psi ,}$
де ${\displaystyle \psi (x)=f(x)-D_{2}(y).}$ Отже, отримаємо розв'язок ${\displaystyle y=D_{1}^{-1}[\psi ].}$
Розглянемо даний метод детальніше на прикладі рівняння, яке запишемо у вигляді:
${\displaystyle y^{(n)}(x)+\sum _{i=1}^{n}a_{i}y^{(n-1)}(x)=f(x)-\sum _{i=m+1}^{n}a_{i}y^{(n-1)}(x).}$
Після заміни змінних ${\displaystyle u(x)=y^{(n-m)}(x),{\dot {u}}(x)=y^{(n-m+1)}(x),...,u^{(m)}(x)=y^{(n)}(x)}$ отримаємо рівняння ${\displaystyle m}$-го порядку:
${\displaystyle u^{m}(x)+\sum _{i=1}^{m}a_{i}u^{(n-1)}(x)=\psi (x)}$,
де ${\displaystyle \psi (x)=f(x)-\sum _{i=m+1}^{n}a_{i}y^{(n-i)}(x).}$ Коли ми перейдемо до еквівалентної системи диференціальних рівнянь і використаємо фундаментальний розв'язок стосовно канонічної системи диференціальних рівнянь, отримуємо рівняння з ядром експоненціального виду:
${\displaystyle u(x)=e^{At}u_{0}+\int _{0}^{x}e^{As}\Phi (a,u,s)\,ds,}$
де ${\displaystyle u(x)=[{\dot {u}}(x),{\ddot {u}}(x),\ldots ,u^{m}(x)],u_{0}(x)=[{\dot {u}}(0),{\ddot {u}}(0),\ldots ,u^{m}(0)],\Phi =(a,u,s)=[0,0,\ldots ,\psi (x)].}$
Матриця ${\displaystyle A}$ порядку ${\displaystyle m}$ матиме вигляд:

${\displaystyle A=\left({\begin{array}{cccc}0&1&\ldots &0\\0&0&\ldots &0\\\ldots &\ldots &\ldots &\ldots \\-a_{m}&-a_{m-1}&\ldots &-a_{1}\end{array}}\right).}$

Врахуємо залежність ${\displaystyle y(x)=\int _{0}^{x}\ldots \int _{0}^{x}u(s)\,ds={\frac {1}{(n-m-1)!}}\int _{a}^{x}(x-s)^{n-m-1}u(s)\,ds}$. Тоді, легко бачити, що еквівалентне перетворення відбувається шляхом зміни значення ${\displaystyle m\in {\overline {1,n.}}}$^[9]

Нехай у методі перетворення з розщепленням зробимо заміну ${\displaystyle m=n}$. Тоді розщеплення оператора ${\displaystyle D}$ зводиться до розв'язання вихідного рівняння відносно старшої похідної. При цьому, розв'язок рівняння полягає у послідовному ${\displaystyle n}$-кратному інтегруванні, в результаті якого отримаємо інтегральне рівняння виду:
${\displaystyle y(x)+\int _{0}^{x}K(x-s)y(s)\,ds=F(x),}$
де
${\displaystyle K(x-s)=\sum _{i=1}^{n}q_{i}{\frac {(x-s)^{i-1}}{(i-1)!}},}$
${\displaystyle F(x)=\int _{0}^{t}{\frac {(x-s)^{n-1}}{(n-1)!}}f(s)\,ds+\sum _{i=0}^{n-1}C_{i}{\frac {x^{i}}{i!}}+C_{0}\sum _{i=1}^{n-1}q_{i}{\frac {x^{i}}{i!}}+\ldots +C_{1}\sum _{i=1}^{n-2}q_{i}{\frac {x^{i+1}}{(i+1)!}}+\ldots +C_{n-2}{\frac {x^{n-1}}{(n-1)!}}}$.^[10]

Лінійні динамічні системи зазвичай описують системами лінійних звичайних диференціальних рівнянь, рівнянь у частинних похідних і лінійними різницевими та інтегральними рівняннями. Для лінійних моделей точні розв'язки можна знайти в аналітичній формі. Більше того, в деяких випадках, нелінійні процеси апроксимують лінійними.
Нелінійні динамічні системи називаються хаотичними, якщо їхня поведінка є випадковою, попри те, що вона визначається детерміністичними законами.^[11]

Нехай ${\displaystyle \delta x(0)}$  — нескінченно мала відстань між двома точками у фазовому просторі, які належать різним фазовим траєкторіям у момент часу ${\displaystyle t=0}$; ${\displaystyle \delta x(t)}$ — відстань між цими точками у момент часу ${\displaystyle t}$. Тоді, запишемо наступне:
${\displaystyle |\delta x(t)|\approx |\delta x(0)|\exp(\lambda t)}$,
де параметр ${\displaystyle \lambda }$ називається показником Ляпунова. Якщо ${\displaystyle \lambda >0}$, тоді дві фазові траєкторії, що виходять з малого околу певної точки простору, з часом розходяться з експоненціальною швидкістю. Із співвідношення ${\displaystyle |\delta x(t)|\approx |\delta x(0)|\exp(\lambda t)}$ отримаємо формулу для розрахунку показника Ляпунова:
${\displaystyle \lambda =\lim _{{N\to \infty },{\delta x(0)\to 0}}{\frac {1}{t}}\ln {\Big |}{\frac {\delta x(t)}{\delta x(0)}}{\Big |}}$.
Цей показник є функцією початкової координати у загальному випадку.
Розгляд хаотичних динамічних систем зручно почати з простих прикладів одновимірних дискретних відображень, які мають вигляд:
${\displaystyle x_{n+1}=f(x_{n})}$.
Остаточний вигляд формули для розрахунку показника Ляпунова у випадку одновимірного відображення, поданого вище, є таким: ${\displaystyle \lambda =\lim _{N\to \infty }{\frac {1}{N}}\sum _{i=0}^{N-1}\ln |{\dot {f}}(x_{i})|}$.^[12]

Розглянемо сукупність n матеріальних точок. Як відомо, положення точки у просторі визначається її радіус-вектором ${\displaystyle {\overrightarrow {r}}=(x,y,z)}$. Щоб визначити положення системи n матеріальних точок у просторі, треба знати n радіус-векторів або 3n координат. Кількість незалежних величин, що визначають положення системи у просторі, називається кількістю степенів вільності системи. У загальному випадку це можуть бути й недекартові координати (полярні, сферичні тощо).
Наведемо приклад механічної системи з одним степенем вільності. Одним із загальних принципів, що дозволяє побудувати рівняння руху, тобто створити математичну модель функціонування динамічної системи, є принцип найменшої дії (Гамільтона). Згідно із ним кожна механічна система характеризується деякою визначеною функцією ${\displaystyle L(q,{\dot {q}},t)}$, де q — узагальнені координати, ${\displaystyle {\dot {q}}}$ — узагальнена швидкість, t — момент часу.^[13]

Наведемо приклад популяції, яка є ізольованою. Популяцією є сукупність індивідів, що можуть давати потомство й піддаються впливу однакових внутрішніх і зовнішніх факторів. Припустимо, що ареал їх проживання обмежений. Основним припущенням, що використовується при побудові математичних моделей динаміки зміни чисельності популяцій, є балансове співвідношення між різними групами у структурі популяцій під впливом факторів різної природи.
Одним з найпростіших прикладів таких моделей є робота Томаса Мальтуса «Досвід закону про народонаселення» (1797). Автор вплинув на формування концепції Чарльза Дарвіна про розуміння природного відбору як рушійної сили еволюції. У цій роботі модель мала вигляд звичайного скалярного лінійного диференціального рівняння зі сталим коефіцієнтом ${\displaystyle {\dot {x}}(t)=kx(t),x(0)=x_{0}>0}$, де ${\displaystyle x(t)}$ — чисельність популяції в момент t, k  — інтенсивність народжуваності (смертності). Розв'язок рівняння має вигляд:${\displaystyle x(t)=x_{0}e^{kt},t\geq 0}$.^[14]