Кирсанов М.Н. Решебник.Теоретическая механика.  М.:ФИЗМАТЛИТ, 2002.
 Иллюстрации и программы к книге

14.1. Задача 4

>    restart:

>    read "ris.m":

     Массы (кг)   

>    ma:=3:mb:=6:

     Жесткость пружин (Н/м)

>     c:=7:

     Кинетические энергии цилиндров  A, B

>    Ta:=3/4*ma*va^2: Tb:=3/4*mb*vb^2:

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

>    T:=Ta+Tb:

     Инерционные коэффициенты (Тарг С.М.[19, c.394])

>    a11:=coeff(diff(T,va),va);

>    a12:=coeff(diff(T,va),vb);

>    a22:=coeff(diff(T,vb),vb);

     Потенциальная энергия

>    Pt:=c/2*xa^2+c/2*(2*xb-xa)^2+c/2*(2*xb)^2+c/2*xb^2:

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

>    Q1:=-diff(Pt,xa):Q2:=-diff(Pt,xb):

       Коэффициенты жесткости

>     c11:=-coeff(Q1,xa);
 c12:=-coeff(Q2,xa);
 c22:=-coeff(Q2,xb);

        Уравнения частот

>     ur:=(c11-a11*omega^2)*(c22-a22*omega^2)-c12^2=0:

       Решение уравнения частот

>     sol:=fsolve(ur,omega=0..infinity);

a11 := 9/2

a12 := 0

a22 := 9

c11 := 14

c12 := -14

c22 := 63

sol := 1.455853300, 2.826942214

     Дифференциальные уравнения колебаний

>    eq1:=a11*diff(xa(t),t$2)=subs(xa=xa(t),xb=xb(t),Q1);

>    eq2:=a22*diff(xb(t),t$2)=subs(xa=xa(t),xb=xb(t),Q2);

eq1 := 9/2*diff(xa(t),`$`(t,2)) = -14*xa(t)+14*xb(t)

eq2 := 9*diff(xb(t),`$`(t,2)) = -63*xb(t)+14*xa(t)

    Численное решение системы  с начальными данными

>    r:=dsolve({eq1,eq2,xa(0)=0,xb(0)=0,D(xa)(0)=10,D(xb)(0)=3},
          {xa(t),xb(t)},type=numeric,
                      output=listprocedure):

>    with(plots):

>    odeplot(r,[[t,xa(t)],[t,xb(t)]],0..4.3,numpoints=150,
                                          labels=[t,x],
                                          color=black);

>    Xa:= subs(r,xa(t)): Xb:= subs(r,xb(t)):

                                 Изображение механизма в движении      
     

Warning, the name changecoords has been redefined

[Maple Plot]

               Подписи точек на рисунке
nam:=[A,B]:
RA:=40: RB:=20: h:=100:
                      Координаты

>      y[1]:=RA:y[4]:=RA:y[2]:=RB:y[3]:=2*RA:

>   
            Количество кадров K

>    K:=12: with(plots):with(plottools):

           Создаем все кадры

>    for i from 0 to K do

>    t:=i/K*4.3:

>    x[1]:=Xa(t):

>    x[2]:=Xb(t)+h:x[4]:=x[2]:
P[i]:=display(Cir(1,RA),Cir(2,RB),
pruzh(-60,x[1],y[1],4,24),
pruzh(x[1],x[4],y[1],4,24),
pruzh(x[4],160,y[4],4,24),
pruzh(x[2],160,y[2],4,24),
Cir(2,2),Cir(1,2),
cir4(1,RA,-Xa(t)/RA),cir4(2,RB,-Xb(t)/RB),
seq(TEXT([x[j]+2,y[j]+8],nam[j]),j=1..2)):
od:

                                 Изображение механизма в движении

>    PP:=display(seq(P[i],i=0..K),insequence=true,
                            thickness=2,
                            scaling=constrained,
                            axes=none):

Warning, the name arrow has been redefined

>    display(PP,Поверхность(-60,0,220,6),Стенка(-60,0,80,6),
Стенка(160,0,48,-6));

[Maple Plot]

>   

>