Методы исследования на устойчивость. Разбор практических задач в Maiple/
В данной статье представлен разбор задач по биологии методами теории устойчивости на примере задачи о конкуренции 2 видов.
Скачать:
Предварительный просмотр:
Решим задачу о конкуренции двух видов.
Подключение необходимых библиотек.
> restart:
> with(plots):
> with(DEtools):
Описание процедуры, цель которой выяснение вопроса устойчивости положения равновесия (x0,y0)
Используется метод линеаризации.
> steady:=proc(x0,y0,a1,b1,c1,a2,b2,c2) local f1,f2,f1_taylor,f2_taylor,A,sob; f1:=unapply(subs({r_1=a1,K_1=b1,a_12=c1,x=u+x0,y=v+y0},f_1),u,v);
> f2:=unapply(subs({r_2=a2,K_2=b2,a_21=c2,x=u+x0,y=v+y0},f_2),u,v);
> readlib(mtaylor):
> f1_taylor:=mtaylor(f1(u,v), [u,v],2);f2_taylor:=mtaylor(f2(u,v), [u,v],2);
> coeff(f1_taylor,u);
>A:=linalg[matrix](2,2,[coeff(f1_taylor,u),coeff(f1_taylor,v),coeff(f2_taylor,u),coeff(f2_taylor,v)]);
> sob:=linalg[eigenvalues](A);
> if sob[1]<=0 and sob[2]<=0 then 1:else 0:fi;
> end:
1 ШАГ. Задание общего вида правых частей Д. У. системы и необходимых характеристик: удельные скорости роста (r_i, i=1,2), соотношения между емкостью среды i-того вида при отсутствии конкурента (K_i, i=1,2) и коэффициентом относительного влияния друг на друга (a_12, a_21).
> f_1:=r_1*x*(1-x/K_1-a_12*y/K_1);
> f_2:=r_2*y*(1-y/K_2-a_21*x/K_2);
> L_1:=(K_1-x)/a_12;
> L_2:=K_2-a_21*x;
2 ШАГ. Задание конкретных значений характеристик модели.
L_1(red) целиком выше L_2(blue)
> L_11:=plot(subs({K_1=4,a_12=.1},L_1),x=0..10,y=0..50):
> L_22:=plot(subs({K_2=6,a_21=4},L_2),x=0..10,y=0..50,color=blue):
> display({L_11,L_22},title=`Red-L_1, blue-L_2`);
> f_11:=unapply(subs({r_1=2,K_1=4,a_12=.1},f_1),x,y);
> f_22:=unapply(subs({r_2=1,K_2=6,a_21=4},f_2),x,y);
3 ШАГ. Нахождение положений равновесия, путем решения системы алгебраических уравнений.
> fixed:=solve({f_11(x,y)=0,f_22(x,y)=0},{x,y});
4 ШАГ . Исследование полученных на предыдущем шаге и удовлетворяющих физическому смыслу задачи положений равновесия на предмет устойчивости.
Устойчивость
> if steady(0,0,2,4,.1,1,6,4)=1 then print(`(0,0)-Steady`):else print(`(0,0)-UnSteady`):fi;
> if steady(0,6,2,4,.1,1,6,4)=1 then print(`(0,6)-Steady`):else print(`(0,6)-UnSteady`):fi;
> if steady(4,0,2,4,.1,1,6,4)=1 then print(`(4,0)-Steady`):else print(`(4,0)-UnSteady`):fi;
5 ШАГ. Графическое изображение поля направлений.
>dfieldplot([D(x)(t)=f_11(x,y,t),D(y)(t)=f_22(x,y,t)],[x(t),y(t)],t=-10..10,x=0..5,y=0..5,title=`{dx/dt=f_11(x,y),dy/dt=f_22(x,y)}`);
6 ШАГ. Задание задачи Коши и изображение фазового портрета. Если необходимо, данный шаг можно повторить несколько раз для разных начальных условий задачи Коши.
>phaseportrait([D(x)(t)=f_11(x,y,t),D(y)(t)=f_22(x,y,t)],[x(t),y(t)],0..10,[[x(0)=.3,y(0)=2]],stepsize=.05,scene=[x(t),y(t)],linecolour=sin(t*Pi/2),method=classical[foreuler],title=`{dx/dt=f_11(x,y),dy/dt=f_22(x,y)} {x(0)=0.3,y(0)=2}`);
7 ШАГ. График решения системы, описывающей модель, и некоторых характерных линий.
> F:=dsolve({D(x)(t)=f_11(x,y,t),D(y)(t)=f_22(x,y,t),x(0)=3,y(0)=3},{x(t),y(t)},type=numeric):
> g_1:=plots[odeplot](F,[t,x(t)],0..5,labels=[t,x],color=red):
> g_2:=plots[odeplot](F,[t,y(t)],0..5,labels=[t,y],color=blue):
> g_3:=plot(4,0..5,color=green):
> display({g_1,g_2,g_3},title=`First kind-red, second kind-blue`);
Аналогичный алгоритм применяется для остальных типов поведения прямых L_1 и L_2.
L_2(blue) целиком выше L_1(red)
> L_11:=plot(subs({K_1=6,a_12=4},L_1),x=0..10,y=0..10):
> L_22:=plot(subs({K_2=4,a_21=.1},L_2),x=0..50,y=0..10,color=blue):
> display({L_11,L_22},title=`Red-L_1, blue-L_2`);
> f_11:=unapply(subs({r_1=2,K_1=6,a_12=4},f_1),x,y);
> f_22:=unapply(subs({r_2=1,K_2=4,a_21=.1},f_2),x,y);
> fixed:=solve({f_11(x,y)=0,f_22(x,y)=0},{x,y});
Устойчивость
> if steady(0,0,2,6,4,1,4,.1)=1 then print(`(0,0)-Steady`):else print(`(0,0)-UnSteady`):fi;
> if steady(0,4,2,6,4,1,4,.1)=1 then print(`(0,4)-Steady`):else print(`(0,4)-UnSteady`):fi;
> if steady(6,0,2,6,4,1,4,.1)=1 then print(`(6,0)-Steady`):else print(`(6,0)-UnSteady`):fi;
>dfieldplot([D(x)(t)=f_11(x,y,t),D(y)(t)=f_22(x,y,t)],[x(t),y(t)],t=-10..10,x=0..5,y=0..5,title=`{dx/dt=f_11(x,y),dy/dt=f_22(x,y)}`);
>phaseportrait([D(x)(t)=f_11(x,y,t),D(y)(t)=f_22(x,y,t)],[x(t),y(t)],0..10,[[x(0)=.2,y(0)=.2]],stepsize=.05,scene=[x(t),y(t)],linecolour=sin(t*Pi/2),method=classical[foreuler],title=`{dx/dt=f_11(x,y),dy/dt=f_22(x,y)} {x(0)=0.2,y(0)=0.2}`);
> F:=dsolve({D(x)(t)=f_11(x,y,t),D(y)(t)=f_22(x,y,t),x(0)=4.5,y(0)=4.5},{x(t),y(t)},type=numeric);
> g_1:=plots[odeplot](F,[t,x(t)],0..6,labels=[t,x],color=red):
> g_2:=plots[odeplot](F,[t,y(t)],0..6,labels=[t,y],color=blue):
> g_3:=plot(4,0..5,color=green):
> display({g_1,g_2,g_3},title=`First kind-red, second kind-blue`);
L_1 и L_2 пересекаются в положительном квадранте и L_1 круче падает чем L_2.
> L_11:=plot(subs({K_1=4.5,a_12=.1},L_1),x=0..10,y=0..30):
> L_22:=plot(subs({K_2=5,a_21=.1},L_2),x=0..50,y=0..10,color=blue):
> display({L_11,L_22},title=`Red-L_1, blue-L_2`);
> f_11:=unapply(subs({r_1=2,K_1=4.5,a_12=.1},f_1),x,y);
> f_22:=unapply(subs({r_2=1,K_2=3,a_21=.1},f_2),x,y);
> fixed:=solve({f_11(x,y)=0,f_22(x,y)=0},{x,y});
Устойчивость
> if steady(0,0,2,4.5,.1,1,3,.1)=1 then print(`(0,0)-Steady`):else print(`(0,0)-UnSteady`):fi;
> if steady(0,3,2,4.5,.1,1,3,.1)=1 then print(`(0,3)-Steady`):else print(`(0,3)-UnSteady`):fi;
> if steady(4.5,0,2,4.5,.1,1,3,.1)=1 then print(`(4.5,0)-Steady`):else print(`(4.5,0)-UnSteady`):fi;
> if steady(4.242424243,2.575757576,2,4.5,.1,1,3,.1)=1 then print(`(4.242424243,2.575757576)-Steady`):else print(`(4.242424243,2.575757576)-UnSteady`):fi;
>dfieldplot([D(x)(t)=f_11(x,y,t),D(y)(t)=f_22(x,y,t)],[x(t),y(t)],t=0..15,x=0..10,y=0..10,title=`{dx/dt=f_11(x,y),dy/dt=f_22(x,y)}`);
> F:=dsolve({D(x)(t)=f_11(x,y,t),D(y)(t)=f_22(x,y,t),x(0)=3,y(0)=3},{x(t),y(t)},type=numeric):
> g_1:=plots[odeplot](F,[t,x(t)],0..5,labels=[t,x],color=red):
> g_2:=plots[odeplot](F,[t,y(t)],0..5,labels=[t,y],color=blue):
> display({g_1,g_2},title=`First kind-red, second kind-blue`);
L_1 и L_2 пересекаются в положительном квадранте и L_2 круче падает чем L_1.
> L_11:=plot(subs({K_1=4.5,a_12=1},L_1),x=0..6,y=0..6):
> L_22:=plot(subs({K_2=5,a_21=2},L_2),x=0..6,y=0..6,color=blue):
> display({L_11,L_22},title=`Red-L_1, blue-L_2`);
> f_11:=unapply(subs({r_1=2,K_1=4.5,a_12=1},f_1),x,y);
> f_22:=unapply(subs({r_2=1,K_2=3,a_21=2},f_2),x,y);
> fixed:=solve({f_11(x,y)=0,f_22(x,y)=0},{x,y});
Устойчивость
> if steady(0,0,2,4.5,1,1,3,2)=1 then print(`(0,0)-Steady`):else print(`(0,0)-UnSteady`):fi;
> if steady(0,3,2,4.5,1,1,3,2)=1 then print(`(0,3)-Steady`):else print(`(0,3)-UnSteady`):fi;
> if steady(4.5,0,2,4.5,1,1,3,2)=1 then print(`(4.5,0)-Steady`):else print(`(4.5,0)-UnSteady`):fi;
>dfieldplot([D(x)(t)=f_11(x,y,t),D(y)(t)=f_22(x,y,t)],[x(t),y(t)],t=0..15,x=0..5,y=0..5,title=`{dx/dt=f_11(x,y),dy/dt=f_22(x,y)}`);
> F:=dsolve({D(x)(t)=f_11(x,y,t),D(y)(t)=f_22(x,y,t),x(0)=6,y(0)=5},{x(t),y(t)},type=numeric):
> g_1:=plots[odeplot](F,[t,x(t)],0..27,labels=[t,x],color=red):
> g_2:=plots[odeplot](F,[t,y(t)],0..27,labels=[t,y],color=blue):
> display({g_1,g_2},title=`First kind-red, second kind-blue`);
