Определение условия нестабильности  системы  с 2 степенями свободы

>

restart;

    Процедура дифференцирования по t любого порядка

>	Dif:=proc(x,i) if i=0 then x else diff(x,t$i) fi end: ------------------------------------------------------------------------------------------

    Процедура дифференцирования кинетической энергии по переменной z

>	DL:=proc(T,z) subs(xx=z,diff(subs(z=xx,T),xx));end:

    -------------------------------------------------------------------------------------------     
  with(LinearAlgebra):   x1:=X(t):x2:=Y(t):

      with(PDEtools,declare):with(plots):declare(x1,x2):

   Процедура получения i-го уравнения  Лагранжа 2-го рода
   (кинетическая энергия T, обобщенные координаты x1,x2, обобщенные силы Q)  

>	Lagr2:=proc(T,Q,i) simplify(diff(DL(T,diff(x||i,t)),t)-DL(T,x||i))-Q: end: -------------------------------------------------------------------------------------------

Warning, the name changecoords has been redefined

Процедура нахождения условия нестабильности с подстановкой ускорений (xtt,ytt) из уравнений движения

>	NC:=proc(M) local A,i,j,Dtr;

>	A:=Matrix(2):

>	for i to 2 do

>	for j to 2 do

>	A[i,j]:=subs(xx=Dif(x||j,M),diff(subs(Dif(x||j,M)=xx,eq[i]),xx));

>

od:

>

od:

     Dtr:=simplify(Determinant(A));

      collect(numer(simplify(subs({diff(x1,t$2)=xtt,diff(x2,t$2)=ytt},Dtr))),X(t));

>	end proc: -------------------------------------------------------------------------------------------

    Константы задачи    

>	mu1:=10: mu2:=10: m1:=3:

>	m2,m3,m4,R:=1$4:

>	I1:=m1*R^2/2:

>	a:=R^2*(m2+3/8*m3):

>	B:=m4+3/8*m3:

>	C:=3/8*R*m3:

>

F:=B/2:

>	G:=C*sin(x2):

>	H:=(I1+a*sin(x2)^2)/2:

>	с,M0:=100$2:

     Кинетическая энергия

>	T:=diff(x1,t)^2*F+diff(x1,t)*diff(x2,t)*G+diff(x2,t)^2*H:

    Обобщенные силы

>	Q1:=-с*x1-mu1*diff(x1,t): Q2:=M0-mu2*x2*diff(x2,t):

    Исследуемая  система уравнения в форме Лагранжа

>	for i to 2 do eq[i]:=Lagr2(T,Q||i,i) od:

   Решение системы дифференциальных уравнений

>	sol:=dsolve({eq[1],eq[2],X(0)=0,D(X)(0)=0, Y(0)=0,D(Y)(0)=1},{X(t),Y(t)},type=numeric,output=listprocedure):

   Построение графиков

>	opt:=0..5,numpoints=500,thickness=2,title=cat("mu2=",convert(mu2,string)):

>	s1:=odeplot(sol,[t,diff(x1,t)],opt,color=black,legend=["Скорость тележки v1"]):

>	s2:=odeplot(sol,[t,diff(x2,t)],opt,color=green,legend=["Угловая скорость маховика v2"]):

>	s3:=odeplot(sol,[t,X(t)],opt,color=blue,legend=["Координата x"]):

>	s4:=odeplot(sol,[t,Y(t)],opt,color=red,legend=["Угол Fi"]):

>	display(s1,s2,s3,s4);

   Выражение ускорений из уравнений движения

>	EQ:=subs(diff(x1,t$2)=xtt,diff(x2,t$2)=ytt,{eq[1],eq[2]}):

>	assign(solve(EQ,{xtt,ytt}));

>	odeplot(sol,[t,NC(0)],opt,color=red,legend=["1/2"]);

>	odeplot(sol,[t,NC(1)],opt,color=magenta,legend=["0/2"]);