Программа 21. Задача 37

>    restart;

>    Graf:=proc(S,L,f,w)

>     [seq(v[S[-1],j]=v[S[1],j]-

>      add(L[i]*sin(f[i]-Pi/2*(j-1))*w[i],i=1..nops(f)),j=1..2)]

>    end:

Неподвижные точки (нулевые скорости)

>    v[C,1],v[C,2],v[A,2],v[K,2]:=0$4:

>    Gr1:=Graf([K,A,P,B,C],[r$4],[Pi/2,Phi,Phi,0],[w2,w2,w3,w3]):

>    eq1:=Gr1[1]; eq2:=Gr1[2];

eq1 := 0 = v[K,1]-r*w2-r*sin(Phi)*w2-r*sin(Phi)*w3

eq2 := 0 = r*cos(Phi)*w2+r*cos(Phi)*w3+r*w3

>    Gr2:=Graf([A,B,C],[2*r,r],[Phi,0],[w0,w3]):

>    eq3:=Gr2[1]; eq4:=Gr2[2];

eq3 := 0 = v[A,1]-2*r*sin(Phi)*w0

eq4 := 0 = 2*r*cos(Phi)*w0+r*w3

>    S:=solve({eq1,eq2,eq3,eq4},{w2,w3,v[K,1],v[A,1]}); assign(S):

S := {w3 = -2*cos(Phi)*w0, w2 = 2*cos(Phi)*w0+2*w0, v[K,1] = 2*r*cos(Phi)*w0+2*r*w0+2*r*sin(Phi)*w0, v[A,1] = 2*r*sin(Phi)*w0}

>    J2:=m[2]*r^2/2:# Момент инерции

>    T1:=m[1]*v[K,1]^2/2:# Кинетическая энергия

>    T2:=m[2]*v[A,1]^2/2+J2*w2^2/2:

>    T:=factor(combine(T1+T2,trig));

T := 1/2*r^2*w0^2*(8*m[1]+8*m[1]*cos(Phi)+4*m[1]*sin(2*Phi)+8*m[1]*sin(Phi)+5*m[2]-m[2]*cos(2*Phi)+4*m[2]*cos(Phi))

>    z1:=diff(T,w0):

>    z2:=diff(T,Phi):

>    with(PDEtools): declare(pfi(t)):

>    Phi:=phi(t);

>    w0:=diff(Phi,t);

pfi(t)*`will now be displayed as`*pfi

Phi := phi(t)

w0 := phi[t]

>    Q:=-expand(M*w2/w0);# Обобщенная сила

>    Уравн:=collect(diff(z1,t)-z2,w0)=Q;

Q := -2*M*cos(phi(t))-2*M

`Уравн` := 1/2*r^2*phi[t]^2*(-8*m[1]*sin(phi(t))+8*m[1]*cos(2*phi(t))+8*m[1]*cos(phi(t))+2*m[2]*sin(2*phi(t))-4*m[2]*sin(phi(t)))+r^2*phi[t,t]*(8*m[1]+8*m[1]*cos(phi(t))+4*m[1]*sin(2*phi(t))+8*m[1]*sin...
`Уравн` := 1/2*r^2*phi[t]^2*(-8*m[1]*sin(phi(t))+8*m[1]*cos(2*phi(t))+8*m[1]*cos(phi(t))+2*m[2]*sin(2*phi(t))-4*m[2]*sin(phi(t)))+r^2*phi[t,t]*(8*m[1]+8*m[1]*cos(phi(t))+4*m[1]*sin(2*phi(t))+8*m[1]*sin...

>    r:=1: m[1]:=1: m[2]:=2: M:=1:

>    НачУсл:=phi(0)=Pi/2,D(phi)(0)=0:

>    Sol:=dsolve({Уравн,НачУсл},phi(t),type=numeric,output=operator):

>    assign(Sol):

>    with(plots):with(plottools):

>    График1:=odeplot(Sol,[t,phi(t)],0..9,thickness=2):

>    График2:=odeplot(Sol,[t,w0],0..9):

>    display(График1,График2);

[Maple Plot]