諧波平衡 matlab,增量諧波平衡法程式

根據網上那個範德波極限環程式編的，為什麼我的幅頻圖是很奇怪的折線

clear all

close all

clc

tic

%=====輸入基本參數(已知條件)===================================

%t=1;

syms t

syms w0

m0=0.2712e6;k0=0.3707e4;c0=0.4439e4;v21=1.6039;

mu1=0.05;mu2=0.005;

m1=m0*mu1;m2=m0*mu2;

xi1=0.1195;xi2=0.2177;%阻尼比

lambda1=0.9266;lambda2=0.8117;%頻率比

w00=sqrt(k0/m0);

w1=lambda1*w00;w2=lambda2*w00;

k1=w1^2*m1;

k2=w1^2*m2;

c1=2*xi1*m1*w1;

c2=2*xi2*m2*w2;

%=======================品質矩陣m=======================================

m=diag([m0/v21 m1 m2]);

%=======================阻尼矩陣c=======================================

c=[c0/v21+c1*v21+c2*v21 -c1*v21 -c2*v21;

-c1*v21 c1*v21 0;

-c2*v21 0 c2*v21];

%========================線性剛度矩陣k====================================

k=[k0/v21+k1*v21 -k1*v21 0;

-k1*v21 k1*v21 0;

0 0 0];

%========================外激勵f====================================

f=[0.5836;0;0]*1e6;

%===============================================================

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

for i=1:3

for j=1:6

S(i,j+6*(i-1))=Cs(j);

end

end

S1=diff(S,t,1);%對S矩陣進行求1階導數

S2=diff(S,t,2);%對S矩陣進行求2階導數

%==================假定初始值============================================

A1=[1 1 1 1 1 1]';

A2=[1 1 1 1 1 1]';

A3=[1 1 1 1 1 1]';

A0=[A1;A2;A3];

%========================三次非線性剛度矩陣kn====================================

nonlinear_k=[k2*v21 0 -k2*v21;0 0 0;-k2*v21 0 k2*v21];

x0=S*A0;

q00=(x0(3)-x0(1))^2;

kn=nonlinear_k*q00;

%======================M=========================================

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

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

%======================C=========================================

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

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

%=====================K==========================================

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

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

%=====================nonlinear_K==========================================

fk_nonlinear=inline(S'*kn*S);

KN=quadv(fk_nonlinear,0,2*pi);

%=============================F==================================

fF=inline(S'*f*cos(t));

F=quadv(fF,0,2*pi);

w0=0.4;

%=============w0=1以及A0=========================================

Kmc=w0^2*M+w0*C+K+3*KN;

R=F-(w0^2*M+w0*C+K+KN)*A0;

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

%===============================================================

deta_A=inv(Kmc)*(R+Rmc*0);     %drtA1(1)

A01=A0+deta_A;

n=1;

tol=1;

epsR=0.001;

while tol>epsR

A0=A01;

%======================M=========================================

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

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

%======================C=========================================

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

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

%=====================K==========================================

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

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

%=====================nonlinear_K==========================================

fk_nonlinear=inline(S'*kn*S);

KN=quadv(fk_nonlinear,0,2*pi);

%=============================F==================================

fF=inline(S'*f*cos(t));

F=quadv(fF,0,2*pi);

%===============================================================

%%%%%帶入推導出的公式

Kmc=w0^2*M+w0*C+K+3*KN;

R=F-(w0^2*M+w0*C+K+KN)*A0;

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

%%%%%

tol=norm(R);

if(n>1000)

disp('疊代步數太多，可能不收斂')

return;

end

%===============================================================

deta_A=inv(Kmc)*(R+Rmc*0);     %drtA1(1)

A01=A0+deta_A;

n=n+1;

end

X0=S*A0;

dX0=S1*A0;

figure(1)

plot(A01);hold on;grid on;

tt=0:.1:100;

x01=subs(X0(1),tt);

x02=subs(X0(2),tt);

x03=subs(X0(3),tt);

dx01=subs(dX0(1),tt);

dx02=subs(dX0(2),tt);

dx03=subs(dX0(3),tt);

figure(2)

plot(x01,dx01,'b','linewidth',2)

hold on

plot(x02,dx02,'b','linewidth',2)

hold on

plot(x03,dx03,'b','linewidth',2)

xlabel('x0')