非线性振动 matlab,数学与非线性科学 - 声振论坛 - 振动,动力学,声学,信号处理,故障诊断 - Powered by Discuz!...

17.png (2.32 KB, 下载次数: 4)

2018-2-28 16:07 上传　　Matlab代码：

clear all

clc

tic

%定义各参数

syms t

w0=3;

epsR=0.001;

m=[1 0；0 1];

epsilon=0.24;

r=0.4;

delta=0.56;

k=[1+epsilon*r -epsilon*r;-epsilon*r 1+epsilon*delta];

Cs=[cos(t) cos(3*t) sin(t) sin(3*t)];

Cs1=diff(Cs，t，1);

S=[cos(t) cos(3*t) sin(t) sin(3*t) 0 0 0 0;0 0 0 0 cos(t) cos(3*t) sin(t) sin(3*t)];

A1=[1 1 1 1]';

A2=[1 1 1 1]';

A0=[A1；A2];

T1=[eye(4，4) zeros(4，4)];

T2=[zeros(4，4) eye(4，4)];

S2=diff(S，t，2);

fm=inline(S'*m*S2);

M=quadv(fm，0，2*pi);

fk=inline(S'*k*S);

K=quadv(fk，0，2*pi);

S1=diff(S，t，1);

c=[1 0；0 1];

fc=inline(S'*c*S1);

C=quadv(fc，0，2*pi);

c3=diag(S*A0).^2;

fc3=inline(S'*c3*S1);

C3=quadv(fc3，0，2*pi);

k2=diag(S*A0).*diag(S1*A0);

fk2=inline(S'*k2*S);

K2=quadv(fk2，0，2*pi);

%代入推导出的公式

Kmc=w0^2*M+epsilon*w0*(C3-C)+K+2*epsilon*w0*K2;

R=-(w0^2*M+epsilon*w0*(C3-C)+K)*A0;

Rmc=-(2*w0*M+epsilon*(C3-C))*A0;

%AA首元素已知a1=0.0，求ww

a1=0.0;

%变换矩阵，使ww变量代替a1

Kmc11=-Rmc(:，1);

Kmcr=[Kmc11 Kmc(:，2:size(Kmc，2))];

%求未知变量

AA=inv(Kmcr)*R;

%drtA1(1)

ww=AA(1);

%drtW(1)

%赋予新变量新值

A01=A0+[a1; AA(2:length(A0)，1)];

%A(1)+drtA(1)

% Aw0=AA+A00;

%A1(0)+drtA1(1)=A1(1)

w01=w0+ww;

%W+drtW(1)

n=1;

tol=1;

while tol>epsR

A0=A01;

w0=w01;

c3=diag(S*A0).^2;

fc3=inline(S'*c3*S1);

C3=quadv(fc3，0，2*pi);

k2=diag(S*A0).*diag(S1*A0);

fk2=inline(S'*k2*S);

K2=quadv(fk2，0，2*pi);

%带入推导出的公式

Kmc=w0^2*M+epsilon*w0*(C3-C)+K+2*epsilon*w0*K2;

R=-(w0^2*M+epsilon*w0*(C3-C)+K)*A0;

Rmc=-(2*w0*M+epsilon*(C3-C))*A0;

%%%%%

tol=norm(R);

if(n>1000)

disp('迭代步数太多，可能不收敛')

return;

end

Kmc11=-Rmc(:，1);

Kmcr=[Kmc11 Kmc(:，2:size(Kmc，2))];

AA=inv(Kmcr)*R;

ww=AA(1);

%A00=[w0；A0(2:6，1)];

A01=A0+[a1；AA(2:length(A0)，1)];

w01=w0+ww;

n=n+1;

end

X0=S*A0;

dX0=S1*A0;

%绘范德波图

tt=0:.1:10;

xo1=subs(X0(1)，tt);

xo2=subs(X0(2)，tt);

dxo1=subs(dX0(1)，tt);

dxo2=subs(dX0(2)，tt);

figure(1)

plot(xo1，dxo1，'b'，'linewidth'，2)

hold on

plot(xo2，dxo2，'b'，'linewidth'，2)

axis([-3 3 -3 3])

title('范德波极限环')

xlabel('x0')

ylabel('dx0')

toc

运行结果：范德波极限环

18.png (38.84 KB, 下载次数: 0)

2018-2-28 16:08 上传　　——以上代码由声振之家会员zhangwenjing分享，代码未经验证。