Задача о колебании  цепи из n звеньев

Кирсанов М.Н., МЭИ(ТУ), 2001

Курс "Моделирование мехатронных систем", семестр 6

    

>   restart: n:=5:  a:=array(1..n,1..n):

Кинетическая энергия вращения звена вокруг центра масс
Tw:=add(w[i]^2,i=1..n)/24:

Кинетическая энергия центра масс
Tv:= add( (add(w[j],j=1..i-1)+w[i]/2)^2,i=1..n)/2:
T:=Tv+Tw:
for i to n do for j to n do  a[i,j]:=coeff(diff(T,w[i]),w[j]): od:od: evalm(6*op(a)):

Потенциальная энергия
 PI:= add( (add(phi[j]^2,j=1..i-1)+(1/2)*phi[i]^2) ,i=1..n)/2:

Обобщенные силы
 Q:=vector(n,[seq( subs(phi[i]=phi[i](t),diff(PI,phi[i])),i=1..n)]):

Уравнения Лагранжа 2 рода
for i to n do eq[i]:=add( diff(phi[j](t),t$2)*a[i,j],j=1..n)=-Q[i]; od:

Численное решение уравнений при нулевых начальных скоростях

A - амплитуда начального отклонения
A:=1.2:
r:=dsolve({seq(eq[i],i=1..n),
for i to n do f[i]:= subs(r,phi[i](t)): od:

Длина звена
L:=1:

Анимация из К кадров
with(plottools):with(plots): K:=50:
for i from 1 to K do
t:=i/K*15:
x[1]:=0:y[1]:=0:
   for j from 2 to n+1 do
   y[j]:=y[j-1]-L*cos(f[j-1](t));
   x[j]:=x[j-1]+L*sin(f[j-1](t));
   od:
opor:=CURVES ([ [-1,0] , [1,0]]):
P[i]:=PLOT (opor,seq(circle([x[m],y[m]],0.02,color=red),m=1..n),seq(CURVES ([ [x[m],y[m]] , [x[m+1],y[m+1]]]),m=1..n) ) :od:

Warning, the names arrow and changecoords have been redefined

                                         Изображение колебаний цепи
display(seq(P[i],i=1..36),insequence=true,thickness=2,scaling=constrained,axes=none);

[Maple Plot]

>   

>