Академия ALT Linux:

Компьютерная математика с Maxima
[+]
Компания ALT Linux
Опубликован: 24.03.2015 | Доступ: свободный | Студентов: 553 / 138 | Длительность: 19:00:00
Темы: Математика, Программное обеспечение, Физика
Специальности: Математик, Преподаватель, Физик
Лекция 6:

Моделирование с Maxima
A
|
версия для печати
< Лекция 5 || Лекция 6: 123456 || Лекция 7 >

6.3.2 Фазовые портреты динамических систем

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

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

Как правило, решение задач Коши для ОДУ и их систем — задача хорошо разработанная и с вычислительной точки зрения довольно простая. На практике чаще встречаются другие, более сложные задачи, в частности, исследование поведения динамической системы в зависимости от начальных условий. При этом в большинстве случаев бывает необходимым изучить только асимптотическое решение ОДУ, т.е. y(t\to\infty), называемое аттрактором. Очень наглядным образом можно визуализировать такую информацию на фазовой плоскости, во многом благодаря тому, что существует всего несколько типов аттракторов, и для них можно построить четкую классификацию.

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

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

Рассмотрим несколько наиболее известных классических примеров динамических систем, имея в виду. Это модели динамики популяций (Лотка-Вольтерры), генератора автоколебаний (Ван дер Поля), турбулентной конвекции (Лоренца) и химической реакции с диффузией (Пригожина). Для изучения динамических систем разработана специальная теория, центральным моментом которой является анализ фазовых портретов, т. е. решений, получающихся при выборе всевозможных начальных условий.

6.3.3 Модель динамики популяций

Модель взаимодействия "хищник—жертва" независимо предложили в 1925–1927 гг. Лотка и Вольтерра. Два дифференциальных уравнения моделируют временную динамику численности двух биологических популяций жертв x и хищников y. Предполагается, что жертвы размножаются с постоянной скоростью, а их численность убывает вследствие поедания хищниками. Хищники же размножаются со скоростью, пропорциональной количеству пищи, и умирают естественным образом.

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

Эта модель Вольтерра-Лотка с логистической поправкой описывается системой уравнений

\frac{d}{d\,t}\,y\left( t\right) &=y\left( t\right) \,\left( a-b\,z\left( t\right) \right) -\alpha\,{y\left( t\right) }^{2} \\ \frac{d}{d\,t}\,z\left( t\right) &=\left(-c + d \, y\left( t\right) \right) \,z\left( t\right) -\alpha\,{z\left( t\right) }^{2}
с условиями заданной численности "жертв" и "хищников" в начальный момент t = 0.

Решая эту задачу при различных значениях, получаем различные фазовые портреты (обычный колебательный процесс и постепенная гибель популяций). Результаты приведены на рис.6.10, рис.6.11 и рис.6.12.

(%i1)	a:4$ b:2.5$ c:2$ d:1$ alpha=0$
	eq1:'diff(y(t),t)=(a-b*z(t))*y(t)-alpha*y(t)^2;
	eq2:'diff(z(t),t)=(-c+d*y(t))*z(t)-alpha*y(t)^2;
	atvalue(y(t),t=0,3); atvalue(z(t),t=0,1);
(\%o6)\ \frac{d}{d\,t}\,\mathrm{y}\left( t\right) =\mathrm{y}\left( t\right) \,\left( 4-2.5\,\mathrm{z}\left( t\right) \right) -\alpha\,{\mathrm{y}\left( t\right) }^{2}\\ (\%o7)\ \frac{d}{d\,t}\,\mathrm{z}\left( t\right) =\left( \mathrm{y}\left( t\right) -2\right) \,\mathrm{z}\left( t\right) -\alpha\,{\mathrm{y}\left( t\right) }^{2}\\ (\%o8)\ 3\\ (\%o9)\ 1
(%i10)	desolve([eq1,eq2],[y(t),z(t)]);
rat: replaced -2.5 by -5/2 = -2.5
rat: replaced -2.5 by -5/2 = -2.5
rat: replaced -2.5 by -5/2 = -2.5
rat: replaced 2.5 by 5/2 = 2.5
(\%o10)$\\ \noindent \begin{math} [\mathrm{y}\left( t\right) =\mathrm{ilt}\left( -\frac{5\,\mathrm{laplace}\left( \mathrm{y}\left( t\right) \,\mathrm{z}\left( t\right) ,t,g2176\right) +2\,\alpha\,\mathrm{laplace}\left( {\mathrm{y}\left( t\right) }^{2},t,g2176\right) -6}{2\,g2176-8},g2176,t\right) ,\\ \mathrm{z}\left( t\right) =\mathrm{ilt}\left( \frac{\mathrm{laplace}\left( \mathrm{y}\left( t\right) \,\mathrm{z}\left( t\right) ,t,g2176\right) -\alpha\,\mathrm{laplace}\left( {\mathrm{y}\left( t\right) }^{2},t,g2176\right) +1}{g2176+2},g2176,t\right) ] \end{math}

Очевидная проблема — неразрешимость данной системы в явном виде методом преобразования Лапласа, т.к. она нелинейна.

Используем численный метод Рунге-Кутта из пакета dynamics.

Результаты решения для значений \alpha=0 и \alpha=0.02 представлены на рис.6.11 и рис.6.12.

Фазовый портрет для системы Лотка-Вольтерра

увеличить изображение
Рис. 6.10. Фазовый портрет для системы Лотка-Вольтерра
Решения системы Лотка-Вольтерра в зависимости от времени (a = 0)

увеличить изображение
Рис. 6.11. Решения системы Лотка-Вольтерра в зависимости от времени (a = 0)
Решения системы Лотка-Вольтерра в зависимости от времени (a = 0, 1)

увеличить изображение
Рис. 6.12. Решения системы Лотка-Вольтерра в зависимости от времени (a = 0, 1)

Рассмотрим командный файл для задачи моделирования системы Лотка-Вольтерра в Maxima:

a:4; b:2.5; c:2; d:1; alpha1:0;
ode1:(a-b*x)*y-alpha1*x^2$ ode2:(-c+d*y)*x-alpha1*y^2$
alpha2:0.02;
ode3:(a-b*x)*y-alpha2*x^2$ ode4:(-c+d*y)*x-alpha2*y^2$
load("dynamics");
t1:[]$ xg1:[]$ yg1:[]$ t2:[]$ xg2:[]$ yg2:[]$
l1:rk([ode1,ode2],[y,x],[1,3],[t,0,9,0.01])$
l2:rk([ode3,ode4],[y,x],[1,3],[t,0,9,0.01])$
for i:1 thru length(l1) do(t1:append(t1,[l1[i][1]]),
	xg1:append(xg1,[l1[i][2]]),yg1:append(yg1,[l1[i][3]]));
for i:1 thru length(l2) do(t2:append(t2,[l2[i][1]]),
	xg2:append(xg2,[l2[i][2]]),yg2:append(yg2,[l2[i][3]]));
load("draw");
draw2d(terminal='eps, file_name="lotka1",grid=true,xlabel = "x",
	ylabel = "y", title="Lotka-Volterra system, phaze portrait",
	key= "alpha=0",points_joined = true, point_type = none,
	line_width = 4,color = black, points(xg1,yg1),
	points_joined = true, color = black,point_type = none,
	line_width = 1,key="alpha=0.02", points(xg2,yg2))$
draw2d(terminal='eps, file_name="lotka2",grid=true,xlabel = "t",
	ylabel = "x,y", title="Lotka-Volterra system, alpha=0",
	key= "x(t)",points_joined = true, line_width = 1,
	color = black,point_type = none, points(t1,xg1),
	points_joined = true, line_width = 4, point_type = none,
	color = black, key= "y(t)", points(t1,yg1))$
draw2d(terminal='eps, file_name="lotka3",grid=true,xlabel = "t",
	ylabel = "x,y", title="Lotka-Volterra system, alpha=0.02",
	key= "x(t)",points_joined = true, point_type = none,
	line_width = 1, color = black, points(t2,xg2),
	points_joined = true, line_width = 4, point_type = none,
	color = black, key= "y(t)", points(t2,yg2))$

Дифференциальные уравнения формируются символьными выражениями, определяющими правые части. Порядок следования выражений для расчёта правых частей ОДУ в первом списке функции rk должен соответствовать друг другу. Следует обратить внимание, что результат выполнения функции rk — двухуровневый список (каждый элемент списков l1 и l2 — также список, содержащий значение независимой переменной, и соответствующие значения искомых функций), поэтому оп преобразуется в векторы, используемые для построения графических иллюстраций.

6.3.4 Движение твердого тела

Рассмотрим пример построения трехмерного фазового портрета. Находим решение задачи Эйлера свободного движения твердого тела:

\frac{dx_1}{dt}&=x_2x_3, \\ \frac{dx_2}{dt} &=-x_1x_3,\\ \frac{dx_3}{dt}&=-0.51 x_1x_2.
(%i1)	eq1:'diff(x1(t),t)=x2(t)*x3(t);
	eq2:'diff(x2(t),t)=-x1(t)*x3(t);
	eq3:'diff(x3(t),t)=-0.51*x1(t)*x2(t);
(\%o1)\ \frac{d}{d\,t}\,x1\left( t\right) =x2\left( t\right) \,x3\left( t\right)\\ (\%o2)\ \frac{d}{d\,t}\,x2\left( t\right) =-x1\left( t\right) \,x3\left( t\right)\\ (\%o3)\ \frac{d}{d\,t}\,x3\left( t\right) =-0.51\,x1\left( t\right) \,x2\left( t\right)
(%i4)	load("dynamics")$ l: rk([y*z, -x*z,0.51*x*y],
	[x,y,z],[1,2,3],[t,0,4,0.1])$

Фазовый портрет для данной динамической системы (трехмерная кривая) представлен на рис.6.13.

Фазовый портрет трехмерной динамической системы

увеличить изображение
Рис. 6.13. Фазовый портрет трехмерной динамической системы

6.3.5 Аттрактор Лоренца

Одна из самых знаменитых динамических систем предложена в 1963 г. Лоренцем в качестве упрощенной модели конвективных турбулентных движений жидкости в нагреваемом сосуде тороидальной формы. Система состоит из трех ОДУ и имеет три параметра модели. Задаём правые части уравнений модели Лоренца:

(%i2)	eq: [s*(y-x), x*(r-z) -y, x*y - b*z];
(\%o2)\quad [s\,\left( y-x\right) ,x\,\left( r-z\right) -y,x\,y-b\,z]

Задаём временные параметры решения

(%i3)	t_range: [t,0,50,0.05];
(\%o3)\quad [t,0,50,0.05]

Задаём параметры системы

(%i4)	s: 10.0; r: 28.0; b: 2.6667;
(\%o6)\quad 10.028.02.6667

Задаём начальные значения x,y,z

(%i7)	init: [1.0,0,0];
(\%o7)\quad [1.0,0,0]

Выполняем решение

(%i8)	sol: rk(eq, [x,y,z],init,t_range)$
(%i9)	len:length(sol);
(\%o9)\quad 1001
Выделяем компоненты решения и строим графические иллюстрации:
(%i10)	t:makelist(sol[k][1],k,1,len)$
	x:makelist(sol[k][2],k,1,len)$
	y:makelist(sol[k][3],k,1,len)$
	z:makelist(sol[k][4],k,1,len)$
	plot2d([discrete,t,x])$ plot2d([discrete,t,y])$

Результаты решения (хаотические колебания x,y,z) представлен на рис.6.14 и рис.6.15 (фазовый портрет системы). На рисунках объединены в одних осях кривые x(t), y(t), z(t).

Пример формирования динамического хаоса (аттрактор Лоренца)

увеличить изображение
Рис. 6.14. Пример формирования динамического хаоса (аттрактор Лоренца)
Трехмерный фазовый портрет (аттрактор Лоренца)

увеличить изображение
Рис. 6.15. Трехмерный фазовый портрет (аттрактор Лоренца)

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

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

6.3.6 Модель автоколебательной системы: уравнение Ван дер Поля

Рассмотрим решение уравнения Ван дер Поля, описывающего электрические колебания в замкнутом контуре, состоящем из соединенных последовательно конденсатора, индуктивности, нелинейного сопротивления и элементов, обеспечивающих подкачку энергии извне. Неизвестная функция времени y(t) имеет смысл электрического тока, а в параметре \mu заложены количественные соотношения между составляющими электрической цепи, в том числе и нелинейной компонентой сопротивления:

\frac{d^2y(t)}{dx^2}-\mu (1-y(t)^2) \frac{dy(t)}{dt}+y(t) = 0
Решение уравнения Ван дер Поля

увеличить изображение
Рис. 6.16. Решение уравнения Ван дер Поля

Решением уравнения Ван дер Поля являются колебания, вид которых для \mu=1 показан на рис.6.16. Они называются автоколебаниями и принципиально отличаются от рассмотренных ранее (например, численности популяций в модели Вольтерpa) тем, что их характеристики (амплитуда, частота, спектр) не зависят от начальных условий, а определяются исключительно свойствами самой динамической системы. Через некоторое время расчетов после выхода из начальной точки решение выходит на один и тот же цикл колебаний, называемый предельным циклом. Аттрактор типа предельного цикла является замкнутой кривой на фазовой плоскости. К нему асимптотически притягиваются все окрестные траектории, выходящие из различных начальных точек, как изнутри (рис.6.17), так и снаружи предельного цикла.

Фазовый портрет уравнения Ван дер Поля

увеличить изображение
Рис. 6.17. Фазовый портрет уравнения Ван дер Поля

Использованный командный файл Maxima (для построения графической иллюстрации использован пакет draw):

load("dynamics")$ load("draw")$
mu:1$ s:rk ([v,mu*(1-y^2)*v-y],[y,v],[1,0],[t,0,40,0.2])$
time:makelist(s[k][1],k,1,length(s))$
y:makelist(s[k][2],k,1,length(s))$
v:makelist(s[k][3],k,1,length(s))$
draw2d(points_joined = true, point_type=6, key= "y-v",
	xlabel="y",ylabel="v",points(y,v), terminal=eps)$
Дальше >>
< Лекция 5 || Лекция 6: 123456 || Лекция 7 >

Вопросы и ответы

вопросов: 0