Математическая модель в виде дифференциальных уравнений

Дифференциальные уравнения и математическое моделирование

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ И МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ

научный руководитель канд. физ.-мат. наук, проф.

Лесосибирский педагогический институт – филиал ФГАОУ «Сибирский федеральный университет»

Для изучения объектов или процессов, протекающих в окружающем нас мире, широко используются методы математического моделирования. Математические модели являются мощным средством познания окружающего мира. При этом следует заметить, что построенная математическая модель не может отразить все многообразные и сложные черты изучаемого явления. При моделировании что-то является главным, а что-то – второстепенным, чем можно пренебречь.

Изучение большого круга задач естествознания, техники и механики, биологии, медицины и других отраслей научных знаний показывает, что решение многих из них сводится к математическому моделированию процессов в виде формулы, т. е. в виде функциональной зависимости.

Так, например, некоторые процессы в радиотехнике, кинетика химических реакций, динамика биологических популяций, движение космических объектов, модели экономического развития исследуются с помощью уравнений, в которых кроме независимых переменных и неизвестных функций этих переменных, содержатся производные неизвестных функций (или их дифференциалы). Такие уравнения называются дифференциальными.

Вот почему возможности применения дифференциальных уравнений для решения задач по дисциплинам естественно – научного цикла довольно широки. Обыкновенные дифференциальные уравнения моделируют явления и процессы, которые описываются одной функцией или вектор-функцией одного переменного.

В математическое исследование любой задачи реального мира можно выделить три основных этапа:

1. построение математической модели явления;

2. изучение этой математической модели и получение решения соответствующей математической задачи;

3. приложение полученных результатов к практическому вопросу, из разрешения которого возникла данная математическая модель, и отыскание других вопросов, к которым она применима.

В таблице представлены основные области наук, в которых какое-либо явление или процесс можно записать в виде дифференциального уравнения.

Характеристика составления математической модели

Пример математической модели

1. Установить величины, изменяющиеся в данном явлении, и выявить физические законы, связывающие их.

2. Выбрать независимую переменную и функцию этой искомой переменной.

3. Исходя из условий задачи, определить начальные или краевые условия.

4. Выразить все фигурирующие в условии задачи величины через независимую переменную, искомую функцию и производные этой функции.

5. Исходя из условий задачи и физического закона, которому подчиняется данное явление, составить дифференциальное уравнение.

6. Найти общее решение или общий интеграл дифференциального уравнения.

7. По начальным или краевым условиям найти частное решение.

8. Исследовать полученное решение.

А) Первый закон Ньютона:

Б) Уравнение показательного роста и показательного убывания:

где – некоторая константа.

В) Уравнение гармонического колебания:

где – положительная постоянная.

где – радиус Земли, – расстояние между центрами метеороида и Земли, – ускорение свободного падения.

1. Сделать чертёж и ввести обозначения;

2. Отделить условия. Имеющие место в произвольной точке искомой линии, от условий, выполняющихся лишь в отдельных точках;

3. Выразить все упомянутые в задаче величины через координаты произвольной точки и через значение производной в этой точке, учитывая геометрический смысл производной;

4. По условию задачи составить дифференциальное уравнение;

5. Найти общее решение этого уравнения и получить из него с помощью начальных условий уравнение искомой линии.

Формула зеркала, собирающего все параллельные лучи в одну точку:

При создании математической модели используют физические закономерности, выявленные при экспериментальном изучении объекта моделирования. Так, например, математическая модель кровообращения основано на законах гидродинамики.

Модель хищник – жертва:

где – положительные константы, – число хищников, – число жертв.

Сущность химических реакций сводится к разрыву связей в исходных веществах и возникновению новых связей в продуктах реакции. При этом общее число атомов каждого элемента до и после реакции остаётся постоянным.

Закон действующих масс:

где – концентрации веществ

– коэффициент пропорциональности.

Базовая математическая модель в области финансов формулируется в терминах стохастических процессов, приводящих, таким образом, к стохастическим дифференциальным уравнениям. Время и недостоверность являются главными элементами моделирования финансового поведения экономических агентов.

Модель фондового (биржевого) ценообразования:

где и – постоянные коэффициенты, связанные с характеристиками модели.

Компартментальное моделирование распространено в медицине и биологии. Согласно определению американского фармаколога и биохимика Шеппарда компартмент — это некоторое количество вещества, выделяемое в биологической системе и обладающее свойством единства, поэтому в процессах транспорта и химических преобразований его можно рассматривать как целое. Например, в качестве особых компартментов рассматривают весь кислород в легких, всю углекислоту в венозной крови, количество введенного препарата в межклеточной жидкости, запас гликогена в печени и т. п. Модели, в которых исследуемая система представляется в виде совокупности компартментов, потоков вещества между ними, а также источников и стоков всех веществ, называются компартментальными.

Модель роста опухоли:

где – концентрация опухолевых клеток, – внеклеточная матрица (например, IV тип коллагена).

В заключение отметим, что математическая модель является основой математически оформленной теории того или иного явления, а аппарат дифференциальных уравнений нашел большое применение в математическом моделировании.
Результативность математического моделирования подтверждена всей человеческой практикой, это сильное средство научного исследования, которое используют в каждой конкретной области науки.

Численное решение математических моделей объектов заданных системами дифференциальных уравнений

Введение:

При математическом моделировании ряда технических устройств используются системы дифференциальных нелинейных уравнений. Такие модели используются не только в технике, они находят применение в экономике, химии, биологии, медицине, управлении.

Исследование функционирования таких устройств требуют решения указанных систем уравнений. Поскольку основная часть таких уравнений являются нелинейными и нестационарными, часто невозможно получить их аналитическое решение.

Возникает необходимость использовать численные методы, наиболее известным из которых является метод Рунге — Кутты [1]. Что касается Python, то в публикациях по численным методам, например [2,3], данных по применение Рунге — Кутты крайне мало, а по его модификации — методу Рунге-Кутта-Фельберга вообще нет.

В настоящее время, благодаря простому интерфейсу, наибольшее распространение в Python имеет функцию odeint из модуля scipy.integrate. Вторая функция ode из этого модуля реализует несколько методов, в том числе и упомянутый пятиранговый метод Рунге-Кутта-Фельберга, но, вследствие универсальности, имеет ограниченное быстродействие.

Целью настоящей публикации является сравнительный анализ перечисленных средств численного решения систем дифференциальных уравнений с модифицированным автором под Python методом Рунге-Кутта-Фельберга. В публикации так же приведены решения по краевым задачам для систем дифференциальных уравнений (СДУ).

Краткие теоретические и фактические данные по рассматриваемым методам и программным средствам для численного решения СДУ

Для одного дифференциального уравнения n – го порядка, задача Коши состоит в нахождении функции, удовлетворяющей равенству:

и начальным условиям

Перед решением эта задача должна быть переписана в виде следующей СДУ

(1)

с начальными условиями

Модуль имеет две функции ode() и odeint(), предназначенные для решения систем обыкновенных дифференциальных уравнений (ОДУ) первого порядка с начальными условиями в одной точке (задача Коши). Функция ode() более универсальная, а функция odeint() (ODE integrator) имеет более простой интерфейс и хорошо решает большинство задач.

Функция odeint() имеет три обязательных аргумента и много опций. Она имеет следующий формат odeint(func, y0, t[,args=(), . ]) Аргумент func – это имя Python функции двух переменных, первой из которых является список y=[y1,y2. yn], а второй – имя независимой переменной.

Функция func должна возвращать список из n значений функций при заданном значении независимого аргумента t. Фактически функция func(y,t) реализует вычисление правых частей системы (1).

Второй аргумент y0 функции odeint() является массивом (или списком) начальных значений при t=t0.

Третий аргумент является массивом моментов времени, в которые вы хотите получить решение задачи. При этом первый элемент этого массива рассматривается как t0.

Функция odeint() возвращает массив размера len(t) x len(y0). Функция odeint() имеет много опций, управляющих ее работой. Опции rtol (относительная погрешность) и atol (абсолютная погрешность) определяют погрешность вычислений ei для каждого значения yi по формуле

Они могут быть векторами или скалярами. По умолчанию

Вторая функция модуля scipy.integrate, которая предназначена для решения дифференциальных уравнений и систем, называется ode(). Она создает объект ОДУ (тип scipy.integrate._ode.ode). Имея ссылку на такой объект, для решения дифференциальных уравнений следует использовать его методы. Аналогично функции odeint(), функция ode(func) предполагает приведение задачи к системе дифференциальных уравнений вида (1) и использовании ее функции правых частей.

Отличие только в том, что функция правых частей func(t,y) первым аргументом принимает независимую переменную, а вторым – список значений искомых функций. Например, следующая последовательность инструкций создает объект ODE, представляющий задачу Коши.

При построении численных алгоритмов будем считать, что решение этой дифференциальной задачи существует, оно единственно и обладает необходимыми свойствами гладкости.

При численном решении задачи Коши

(2)

(3)

по известному решению в точке t =0 необходимо найти из уравнения (3) решение при других t. При численном решении задачи (2),(3) будем использовать равномерную, для простоты, сетку по переменной t с шагом т > 0.

Приближенное решение задачи (2), (3) в точке обозначим . Метод сходится в точке если при . Метод имеет р-й порядок точности, если , р > 0 при . Простейшая разностная схема для приближенного решения задачи (2),(3) есть

(4)

При имеем явный метод и в этом случае разностная схема аппроксимирует уравнение (2) с первым порядком. Симметричная схема в (4) имеет второй порядок аппроксимации. Эта схема относится к классу неявных — для определения приближенного решения на новом слое необходимо решать нелинейную задачу.

Явные схемы второго и более высокого порядка аппроксимации удобно строить, ориентируясь на метод предиктор-корректор. На этапе предиктора (предсказания) используется явная схема

(5)

а на этапе корректора (уточнения) — схема

В одношаговых методах Рунге—Кутта идеи предиктора-корректора реализуются наиболее полно. Этот метод записывается в общем виде:

(6),

Формула (6) основана на s вычислениях функции f и называется s-стадийной. Если при имеем явный метод Рунге—Кутта. Если при j>1 и то определяется неявно из уравнения:

(7)

О таком методе Рунге—Кутта говорят как о диагонально-неявном. Параметры определяют вариант метода Рунге—Кутта. Используется следующее представление метода (таблица Бутчера)

Одним из наиболее распространенных является явный метод Рунге—Кутта четвертого порядка

(8)

Метод Рунге—Кутта— Фельберга

Привожу значение расчётных коэффициентов метода

(9)

С учётом(9) общее решение имеет вид:

(10)

Это решение обеспечивает пятый порядок точности, остаётся его адаптировать к Python.

Вычислительный эксперимент по определению абсолютной погрешности численного решения нелинейного дифференциального уравнения с использованием обеих функций def odein(),def oden() модуля scipy.integrate и адаптированного к Python методов Рунге—Кутта и Рунге—Кутта— Фельберга

Адаптированные к Python методы Рунге—Кутта и Рунге—Кутта— Фельберга имеют меньшую абсолютную, чем решение с применением функции odeint, но большую, чем с использованием функции edu. Необходимо провести исследование быстродействия.

Численный эксперимент по сравнению быстродействия численного решения СДУ при использовании функции ode с атрибутом dopri5 (метод Рунге – Кутты 5 порядка) и с использованием адаптированного к Python метода Рунге—Кутта— Фельберга

Сравнительный анализ проведём на примере модельной задачи, приведенной в [2]. Чтобы не повторять источник, приведу постановку и решение модельной задачи из [2].

Решим задачу Коши, описывающую движение тела, брошенного с начальной скоростью v0 под углом α к горизонту в предположении, что сопротивление воздуха пропорционально квадрату скорости. В векторной форме уравнение движения имеет вид

где – радиус вектор движущегося тела, – вектор скорости тела, – коэффициент сопротивления, вектор силы веса тела массы m, g – ускорение свободного падения.

Особенность этой задачи состоит в том, что движение заканчивается в заранее неизвестный момент времени, когда тело падает на землю. Если обозначить , то в координатной форме мы имеем систему уравнений:

К системе следует добавить начальные условия: (h начальная высота), . Положим . Тогда соответствующая система ОДУ 1 – го порядка примет вид:

Для модельной задачи положим . Опуская довольно обширное описание программы, приведу только листинг из комментариев к которому, думаю, будет ясен принцип её работы. В программу добавлен отсчёт времени работы для сравнительного анализа.

Flight time = 1.2316 Distance = 5.9829 Height =1.8542
Flight time = 1.1016 Distance = 4.3830 Height =1.5088
Flight time = 1.0197 Distance = 3.5265 Height =1.2912
Flight time = 0.9068 Distance = 2.5842 Height =1.0240
Время на модельную задачу: 0.454787

Для реализации средствами Python численного решения СДУ без использования специальных модулей, мною была предложена и исследована следующая функция:

def increment(f, t, y, tau
k1=tau*f(t,y)
k2=tau*f(t+(1/4)*tau,y+(1/4)*k1)
k3 =tau *f(t+(3/8)*tau,y+(3/32)*k1+(9/32)*k2)
k4=tau*f(t+(12/13)*tau,y+(1932/2197)*k1-(7200/2197)*k2+(7296/2197)*k3)
k5=tau*f(t+tau,y+(439/216)*k1-8*k2+(3680/513)*k3 -(845/4104)*k4)
k6=tau*f(t+(1/2)*tau,y-(8/27)*k1+2*k2-(3544/2565)*k3 +(1859/4104)*k4-(11/40)*k5)
return (16/135)*k1+(6656/12825)*k3+(28561/56430)*k4-(9/50)*k5+(2/55)*k6

Функция increment(f, t, y, tau) обеспечивает пятый порядок численного метода решения. Остальные особенности программы можно посмотреть в следующем листинге:

Время на модельную задачу: 0.259927

Предложенная программная реализация модельной задачи без использования специальных модулей имеет почти в двое большее быстродействие, чем с функцией ode, однако нельзя забывать, что ode имеет более высокую точность численного решения и возможности выбора метода решения.

Решение краевой задачи с поточно разделёнными краевыми условиями

Приведем пример некоторой конкретной краевой задачи с поточно разделенными краевыми условиями:

(11)

Для решения задачи (11) используем следующий алгоритм:

1. Решаем первые три неоднородные уравнения системы (11) с начальными условиями

Введем обозначение для решения задачи Коши:

2. Решаем первые три однородные уравнения системы (11) с начальными условиями

Введем обозначение для решения задачи Коши:

3. Решаем первые три однородные уравнения системы (11) с начальными условиями

Введем обозначение для решения задачи Коши:

4. Общее решение краевой задачи (11) при помощи решений задач Коши записывается в виде линейной комбинации решений:

где p2, p3 — некоторые неизвестные параметры.

5. Для определения параметров p2, p3, используем краевые условия последних двух уравнений (11), то есть условия при x = b. Подставляя, получим систему линейных уравнений относительно неизвестных p2, p3:
(12)
Решая (12), получим соотношения для p2, p3.

По приведенному алгоритму с применением метода Рунге—Кутта—Фельберга получим следующую программу:

y0[0]= 0.0
y1[0]= 1.0
y2[0]= 0.7156448588231397
y3[0]= 1.324566562303714
y0[N-1]= 0.9900000000000007
y1[N-1]= 0.1747719838716767
y2[N-1]= 0.8
y3[N-1]= 0.5000000000000001
Время на модельную задачу: 0.070878

Вывод

Разработанная мною программа отличается от приведенной в [3] меньшей погрешностью, что подтверждает приведенный в начале статьи сравнительный анализ функции odeint с реализованным на Python метода Рунге—Кутта—Фельберга.

3. Н.М. Полякова, Е.В. Ширяева Python 3. Создание графического интерфейса пользователя (на примере решения методом пристрелки краевой задачи для линейных обыкновенных дифференциальных уравнений). Ростов-на-Дону 2017.

Математическая модель в виде обыкновенных дифференциальных уравнений

2.5. Математическая модель в виде обыкновенных дифференциальных уравнений

Математическая модель в виде одного или нескольких обыкновенных дифференциальных уравнений (ОДУ) широко используются при изучении переходных процессов в системах автоматического регулирования (САР), при описании баллистики летательных аппаратов, а также при описании процессов движения (потоки, частицы, механические элементы).

В простейшем случае модель может иметь вид линейного дифференциального уравнения n-го порядка:

или системы дифференциальных уравнений 1-го порядка

Часто встречаются смешанные задачи, а также нелинейные ОДУ.

Модель, заданная в виде дифференциальных уравнений, должна включать в себя необходимый набор начальных условий:

Рекомендуемые файлы

Исследование моделей, заданных в виде обыкновенных дифференциальных уравнений, осуществляется аналитическими и численными методами. Наиболее полными являются аналитические решения, обеспечивающие всесторонний анализ полученных результатов. Но такие решения получены лишь для ограниченного числа дифференциальных уравнений. Численные методы решения позволяют найти лишь конкретные значения изучаемой функции при заданной комбинации исходных данных. Для анализа модели можно использовать некоторую совокупность решений. Однако, очевидно, что результаты анализа в этом случае могут зависеть от выбора этой совокупности.

Возмущающая сила, вызывающая колебания, зависит от времени f(t). Наряду с возмущающей силой f(t) на груз действует сила инерции , сила вязкого трения , усилие пружины . Все эти силы тормозят движение груза.

Согласно принципу Даламбера сумма всех сил, действующих на груз должна равняться нулю:

. (2.18)

Начальные условия характеризуют начальное положение и начальную скорость груза:

x(0) = x₀; . (2.19)

Уравнение (2.18) совместно с начальными условиями (2.19) представляет собой математическую модель рассматриваемой механической системы.

2.6. Модели, заданные в виде уравнений в частных производных

Ряд задач, связанных с использованием физических полей, приводит к моделям в виде дифференциальных уравнений в частных производных.

Особенностью таких задач является то, что изучаемые параметры изменяются не только во времени, но и зависят от координат x, y, z рассматриваемого пространства. Такие модели называются нестационарными. Модели, в которых параметры не зависят от времени, называются стационарными.

К таким моделям сводятся описания полей температур в элементах конструкции двигателя и полей скоростей при течении жидкости (газа). Уравнениями в частных производных описываются колебания элементов конструкции и поля напряжений, возникающих при работе этих элементов.

Линейное дифференциальное уравнение в частных производных имеет вид

.

Математическая модель, описанная дифференциальными уравнениями в частных производных, должна включать в себя необходимые для решения задачи краевые условия:

1. Должна быть задана область D, ограниченная поверхностью (на плоскости – кривой) G , в которой определяется решение.

2. Должны быть заданы условия на границе G этой области.

В случае нестационарного поля эти граничные условия, так же как и сама область могут меняться во времени.

Граничные условия могут быть 1-го, 2-го и 3-го рода:

а) Граничные условия 1-го рода предусматривают задание на границе величины искомой функции:

– для стационарного поля;

– для нестационарного поля.

б) Граничные условия 2-го рода – предусматривают задание производной искомой функции:

– для стационарного поля;

– для нестационарного поля.

в) Граничные условия 3-го рода – предусматривают комбинации функции и ее производной:

– для стационарного поля;

– для нестационарного поля.

3. Для нестационарных полей должны быть заданы одно или два начальных условия, характеризующих состояние поля в начальный момент времени:

(i = 1, 2, 3).

Здесь x_i – координаты пространства.

Совокупность уравнений и краевых (и начальных) условий полностью определяет модель и позволяет провести ее исследование.

Решение часто задается в виде семейств изолиний F = const (Рис. 2.11).

В качестве примера рассмотрим хорошо изолированный металлический пруток, нагреваемый с одной стороны. С другой стороны помещен измеритель температуры (Рис. 2.12). Величина подогрева x(t) в момент времени t является входным сигналом, а измеряемая на другом конце температура y(t) – выходным сигналом.

Обозначим через x расстояние от измерителя до точки прутка. Температура в этой точке z будет описываться функцией вида

Уравнение теплопроводности для одномерного случая для определения функции z будет иметь вид:

,

где K – коэффициент теплопроводности.

Начальным условием в данном случае является начальное распределение температуры (при t = 0) по прутку: z(0, x) = j(x).

Граничные условия определяются двумя условиями:

а) Нагрев прутка на правом конце

.

б) На левом конце подвод тепла отсутствует

.

в) Показания на измерителе температур (x = 0) в момент времени t определяется следующим выражением

.

Таким образом, для вычисления температуры на расстоянии L от измерителя по формуле для y(t) необходимо проинтегрировать дифференциальное уравнение с учетом начальных и граничных условий, т.е. получить функцию z(t,x). Затем следует проградуировать измеритель температуры, т.е. определить соответствие между x(t) и y(t), задавая различные значения x(t) и вычисляя .

Контрольные вопросы к лекции 5

1. Где используются математические модели в виде обыкновенных дифференциальных уравнений?

2. Что должна включать в себя математическая модель в виде обыкновенных дифференциальных уравнений?

3. Какими методами осуществляется исследование моделей, заданных в виде обыкновенных дифференциальных уравнений?

4. Запишите математическую модель движения груза массой m, закрепленного на вертикальной стенке с помощью пружины жесткостью С и совершающего колебательное движение вдоль оси х в среде с вязкостью n.

5. Какой принцип используется при построении этой модели?

6. К какому типу относится эта модель?

7. Где используются математические модели в виде дифференциальных уравнений в частных производных?

8. Что является особенностью математических моделей в виде дифференциальных уравнений в частных производных?

9. Что должна включать в себя математическая модель в виде дифференциальных уравнений в частных производных?

10. Какого типа бывают граничные условия?

11. Приведите математическую модель распределения температурного поля в металлическом прутке, нагреваемом с одной стороны.