maple非线性振动

1、已知：l ，m ，c ，k ，F=F 0sin(ωt)

求：系统的振动微分方程，ω0，ω=ω0时质点的振幅 解：●建模：质点作圆周运动，杆作定轴转动。 ●Maple > restart:

> J[0]:=m*l^2: > Fe:=-k*3*l*theta(t): >Fd:=-c*2*l*

diff(theta(t),t): > F:=F0*sin(omega*t):

>eq:=J[0]*diff(theta(t ),t$2)=Fd*2*l+Fe*3*l+F *3*l:

>eq:=subs(diff(theta(t),t$2)=DDtheta,

diff(theta(t),t)=Dtheta, theta(t)=theta,eq):

>eq:=m*l^2*DDtheta+4*c*l^2*Dtheta+9*k*l^2*theta=

3*F0*sin(omega*t)*l: > eq:=expand(eq/(m*l^2));

> b:=h/sqrt((omega0^2-omega^2)^2+4*delta^2*omega^2): > b1:=subs(omega=omega0,b);

> omega0:=sqrt(9*k/m);

> delta:=(2*c)/m: > h:=(3*f0)/(m*l):

> B:=simplify(l*b1,symbolic);

答：系统振动微分方程为θ+4c θm

+9k θm

=3F 0sin ⁡(ωt)ml

，

ω0=3 k

m

，B=0

m

4c k

。 2、已知：m=0.5kg ，h=0.1m ，k=0.8kN/m ，β=30° 求：系统的固有频率和振幅，物块的运动方程。

:= eq = + + DDtheta 4c Dtheta 9k θ3F0()

sin ωt :=

b114

h 4δ2

ω0

2

:= ω03

k m

:= B 14f0m c k

●Maple 程序 > restart:

> delta[0]:=m*g*sin(beta)/k: >eq:=m*diff(x(t),t$2)=

m*g*sin(beta)- k*(delta[0]+x): > eq:=lhs(eq)-rhs(eq)=0: >eq:=subs(diff(x(t),t$2)=DDx,

eq):

> eq:=simplify(eq);

> X:=A*sin(omega[0]*t+theta): > omega[0]:=sqrt(k/m): > x[0]:=-delta[0]: > v[0]:=sqrt(2*g*h):

> A:=sqrt(x[0]^2+(v[0]/omega[0])^2): > theta:=arctan(omega[0]*x[0]/v[0]): > m:=0.5: > h:=0.1: > k:=0.8e3: > beta:=Pi/6: > g:=9.8:

>

omega[0]:=evalf(omega[0],4);

> A:=evalf(A,4);

> theta:=evalf(theta,4);

> X:=eval(X);

答：ω0=40rad/s ，A=35.1mm ，

物块的运动方程为x=35.1sin(40t-0.087)mm 。

3、吸引子的仿真。以杜芬方程为例，杜芬方程表示如下

x ′′+c x ′+ax+=Acos ωt

●Maple 程序 > restart:

> with(plots):

> de1:=diff(x(t),t)=y(t):

>de2:=diff(y(t),t)=-a*x(t)-b*x(t)^3-c*y(t)+A*cos(Omeg

:= eq = + m DDx x k 0 := A .03513 := θ-.08724 := X .03513()sin - 40.00t .08724b x 3

:= ω040.00

a*t):

> a:=-1:b:=1:c:=0.15:A:=0.3:Omega:=1: >

duffing:=dsolve({de1,de2,y(0)=-0.5,x(0)=-1},{x(t),y(t )},type=numeric,method=lsode): >

duffplot:=odeplot(duffing,[x(t),y(t)],0..200,numpoint s=4000): > duffplot;

答：杜芬方程相图如图所示。 4、已知：l=l 0−vt 求：摆的运动方程

解：●建模：小球作平面运动 自由度f=1， 取广义坐标φ

●Maple 程序 > restart: > x[rho]:=l:

> x[phi]:=l*phi:

> x[rho]:=subs(l=l(t),x[rho]):

> x[phi]:=subs(phi=phi(t),x[phi]): > v[rho]:=diff(x[rho],t): > v[phi]:=diff(x[phi],t):

> V:=vector([v[rho],v[phi]]): > v[A]:=sqrt(v[rho]^2+v[phi]^2): > T:=1/2*m*v[A]^2;

> T:=subs(diff(phi(t),t)=Dphi,phi(t)=phi,T):

:= T 12m ⎛⎝ ⎫⎭

⎪⎪⎪ + ⎛⎝ ⎫⎭⎪⎪∂∂t ()l t 2l 2⎛⎝ ⎫⎭⎪⎪∂∂t ()φt 2