增量谐波平衡法（IHB）

Hello everyone, I am a beginner in Matlab programming and have been trying to solve the differential equation system for incremental harmonic balance (IHB) given below. I have a MATLAB program, but there seems to be a problem and the analysis results are not satisfactory. We would greatly appreciate any assistance or related work. Thank you in advance.

matlab code
clear all
close all
clc
tic
syms t
%========输入基本参数（已知条件）======
%========duffing 方程得参数===========

a=0.95; %负刚度弹簧在静平衡位置的长度无量纲
lambda=0.085;
lambda_s=-4*lambda*(a-1)/(4*lambda*(a-1)+a);
miu=0.00001; %质量比
xi=0.05;%正刚度阻尼比
xi_b1=0.001;%负刚度装置阻尼比
z=0.06;%幅值
p0=-4*lambda*(1/(a^2)^(1/2)-1);
p1=(2*lambda)/(a^2)^(3/2);
p2=-(3*lambda)/(2*(a^2)^(5/2));
p3=(5*lambda)/(4*(a^2)^(7/2));
p4=-(35*lambda)/(32*(a^2)^(9/2));
p5=(63*lambda)/(64*(a^2)^(11/2));

m=[1,0;0,miu];% m惯性项系数
k=[1+lambda_s,-lambda_s;-lambda_s,lambda_s+p0];% k一次项系数
f=[z;0];% f激励幅值
c=[2*xi,0;0,2*xi_b1];% c阻尼系数

%=====控制参数
omg0=0.005;domg=0.0183;%频率比初始值与增量
%%%%%%能否收敛，delta_s和Nd的取值至关重要
%Nd一般要大于2，易收敛时Nd不宜太大，Nd越小，取值点越密集。非线性较强，Nd取值应稍大一些
delta_s=0.02;%弧长增量值
Nd=1;
Num_Pre_step=4; %预测解需要预测Num_Pre_step步
Num_Incremental_step=160; %程序总共计算Num_Incremental_step步

%========设置谐波矩阵=================

Cs=[ cos(t) sin(t) cos(2*t) sin(2*t)];
S=blkdiag(Cs,Cs);
S1=diff(S,t,1);
S2=diff(S,t,2);
%========设置A矩阵初值================
A1=[ 0.1 0.1 0.1 0.1];%上部结构位移响应的谐波系数
A2=[ 0.1 0.1 0.1 0.1];%调谐装置位移响应的谐波系数
A=[A1,A2]’;%傅里叶系数矩阵
length_A=length(A);
%========质量矩阵m====================

%========刚度矩阵k====================

%========阻尼矩阵k====================

%========外激励矩阵f==================

%========非线性刚度矩阵===============

A0=[zeros(1,length(A1)),A(1:length(A2),:)’]’;
q=(S*A0)’;
kn3=[0;p1]*q.^2;
kn5=[0;p2]*q.^4;
kn7=[0;p3]*q.^6;
kn9=[0;p4]*q.^8;
kn11=[0;p5]*q.^10;

%================
%========积分过程==
fm=inline(S’*m*S2);
M=quadv(fm,0,2*pi);%质量矩阵
fk=inline(S’*k*S);
K=quadv(fk,0,2*pi);%线性刚度矩阵
fc=inline(S’*c*S1);
C=quadv(fc,0,2*pi);%阻尼矩阵
fkn3=inline(S’*kn3*S);
KN3=quadv(fkn3,0,2*pi);%非线性矩阵
fkn5=inline(S’*kn5*S);
KN5=quadv(fkn5,0,2*pi);%非线性矩阵
fkn7=inline(S’*kn7*S);
KN7=quadv(fkn7,0,2*pi);%非线性矩阵
fkn9=inline(S’*kn9*S);
KN9=quadv(fkn9,0,2*pi);%非线性矩阵
fkn11=inline(S’*kn11*S);
KN11=quadv(fkn11,0,2*pi);%非线性矩阵
ff=inline(S’*f*cos(t));
F=quadv(ff,0,2*pi);%激励矩阵

%=========带入公式，公式推导可见陈的书
Kmc=omg0^2*M+omg0*C+K+3*KN3+5*KN5+7*KN7+9*KN9+11*KN11;
R=-F-(omg0^2*M+omg0*C+K+KN3+KN5+KN7+KN9+KN11)*A;
Rmc=-(2*omg0*M+C)*A;
Delta_A=inv(Kmc)*R;

%=======开始迭代过程
epsR=1e-4;
i=1;
X=zeros(length_A+1,4); %建立0矩阵便于保存四个解用于预测
s=zeros(1,3);
Harmonic_A=[]; %用于保存每一个解的谐波项系数
Result_A1=[ ];
for i=1:Num_Pre_step %该部分没有应用弧长法，预先求得一部分解便于弧长法的预测
n=1;tol=1;
while tol>epsR
A=A+Delta_A;
%====回代非线性刚度矩阵，计算下一次循环过程的非线性刚度矩阵
A0=[zeros(1,length(A1)),A(1:length(A2),:)’]’;
q=(S*A0)’;
kn3=[0;p1]*q.^2;
kn5=[0;p2]*q.^4;
kn7=[0;p3]*q.^6;
kn9=[0;p4]*q.^8;
kn11=[0;p5]*q.^10;
fkn3=inline(S’*kn3*S);
KN3=quadv(fkn3,0,2*pi);%非线性矩阵
fkn5=inline(S’*kn5*S);
KN5=quadv(fkn5,0,2*pi);%非线性矩阵
fkn7=inline(S’*kn7*S);
KN7=quadv(fkn7,0,2*pi);%非线性矩阵
fkn9=inline(S’*kn9*S);
KN9=quadv(fkn9,0,2*pi);%非线性矩阵
fkn11=inline(S’*kn11*S);
KN11=quadv(fkn11,0,2*pi);%非线性矩阵

%====再一次计算
Kmc=omg0^2*M+omg0*C+K+3*KN3+5*KN5+7*KN7+9*KN9+11*KN11;
R=-F-(omg0^2*M+omg0*C+K+KN3+KN5+KN7+KN9+KN11)*A;
Rmc=-(2*omg0*M+C)*A;
tol=norm(R);
Delta_A=inv(Kmc)*R;

if(n>60)
disp(‘迭代步数太多，可能不收敛’)
return
end
n=n+1;
end
I3=n-1;
disp([‘增量步=’ num2str(i),’ 迭代次数=’ num2str(I3),’ 本增量步弧长=’ num2str(delta_s),’ 已计算到频率=’ num2str(omg0)])
Harmonic_A=[Harmonic_A;A’];
%%%%%%%%%%%%%保存最新的四组解,便于弧长法预测
for p=1:3
X(:,p)=X(:,p+1);
end
X(1:length_A,4)=A;
X(length_A+1,4)=omg0;
p=0;
%%%%%%%%%%%%%
Frequency(i)=omg0;
% Amplitude11(i)=sqrt(A(2)^2+A(4)^2);
% Amplitude12(i)=sqrt(A(3)^2+A(5)^2);
Amplitude(i)=sqrt(A(2)^2+A(1)^2+A(3)^2+A(4)^2);
Result_A1(:,i)=A;
omg0=omg0+domg;
i=i+1;
end

% 以下是结合了弧长法的增量谐波平衡法
Result_A2=[];
for i=Num_Pre_step+1:Num_Incremental_step %%%%取Num_Incremental_step个增量步
n=1;tol=1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%
for kk=1:3
s(kk)=norm(X(:,kk+1)-X(:,kk));
end
tt(1)=0;tt(2)=s(1);tt(3)=tt(2)+s(2);tt(4)=tt(3)+s(3);tt(5)=tt(4)+delta_s;
PreValue_X=zeros(length_A+1,1);
for ii=1:4
aa=1;
for jj=1:4
if jj~=ii
aa=aa*((tt(5)-tt(jj))/(tt(ii)-tt(jj)));
end
end
PreValue_X=PreValue_X+aa*X(:,ii);
end
A=PreValue_X(1:length_A);
omg0=PreValue_X(length_A+1);
%%%%%%%%%%%%%%%%%%%%%%%%%% 以上这部分为弧长法，根据已经计算得到的四组解预测下一个解的过程，可见陈的书P177
%%%%%计算非线性刚度矩阵
A0=[zeros(1,length(A1)),A(1:length(A2),:)’]’;
q=(S*A0)’;
kn3=[0;p1]*q.^2;
kn5=[0;p2]*q.^4;
kn7=[0;p3]*q.^6;
kn9=[0;p4]*q.^8;
kn11=[0;p5]*q.^10;
fkn3=inline(S’*kn3*S);
KN3=quadv(fkn3,0,2*pi);%非线性矩阵
fkn5=inline(S’*kn5*S);
KN5=quadv(fkn5,0,2*pi);%非线性矩阵
fkn7=inline(S’*kn7*S);
KN7=quadv(fkn7,0,2*pi);%非线性矩阵
fkn9=inline(S’*kn9*S);
KN9=quadv(fkn9,0,2*pi);%非线性矩阵
fkn11=inline(S’*kn11*S);
KN11=quadv(fkn11,0,2*pi);%非线性矩阵

Kmc=omg0^2*M+omg0*C+K+3*KN3+5*KN5+7*KN7+9*KN9+11*KN11;
R=-F-(omg0^2*M+omg0*C+K+KN3+KN5+KN7+KN9+KN11)*A;
Rmc=-(2*omg0*M+C)*A;

%%%%%%%%%%%%判断并寻找控制变量
DELTA_X=PreValue_X-X(:,4); %DELTA_X是预测解与上一增量步的差值，只用来寻找最大值元素，与Delta_X的意义不同。Delta_X是每一个迭代步产生的插值
[~,flag]=max(abs(DELTA_X(1:length_A)));%找到绝对值最大的元素及其下标索引Note_flag,注意要把omg0排除在外，找delta_A中的最大值元素
Note_flag=flag;

%%%%%%%%%%%%%%%%%%%%处理求得的解，得到我们所需要的Delta_A和domg
Kmc(:,Note_flag)=-Rmc(:,1);
Kmc=Kmc(1:length_A,1:length_A);
Delta_X=inv(Kmc)*R;
Delta_X(length_A+1)=0;
Delta_X(length_A+1)=Delta_X(Note_flag);
Delta_X(Note_flag)=0;
Delta_A=Delta_X(1:length_A);
domg=Delta_X(length_A+1);

%A00=[w0；A0(2:6，1)];
A=A+Delta_A;
omg0=omg0+domg;

%%%%%%%%%%%%%%%%%%%%%%%%%下面是每个增量步内的循环迭代
while tol>epsR

%====回代非线性刚度矩阵，计算下一次循环过程的非线性刚度矩阵

A0=[zeros(1,length(A1)),A(1:length(A2),:)’]’;
q=(S*A0)’;
kn3=[0;p1]*q.^2;
kn5=[0;p2]*q.^4;
kn7=[0;p3]*q.^6;
kn9=[0;p4]*q.^8;
kn11=[0;p5]*q.^10;
fkn3=inline(S’*kn3*S);
KN3=quadv(fkn3,0,2*pi);%非线性矩阵
fkn5=inline(S’*kn5*S);
KN5=quadv(fkn5,0,2*pi);%非线性矩阵
fkn7=inline(S’*kn7*S);
KN7=quadv(fkn7,0,2*pi);%非线性矩阵
fkn9=inline(S’*kn9*S);
KN9=quadv(fkn9,0,2*pi);%非线性矩阵
fkn11=inline(S’*kn11*S);
KN11=quadv(fkn11,0,2*pi);%非线性矩阵

%%%%%%%%%%%%%%%%
Kmc=omg0^2*M+omg0*C+K+3*KN3+5*KN5+7*KN7+9*KN9+11*KN11;
R=-F-(omg0^2*M+omg0*C+K+KN3+KN5+KN7+KN9+KN11)*A;
Rmc=-(2*omg0*M+C)*A;
tol=norm(R);

Kmc(:,Note_flag)=-Rmc(:,1);
Kmc=Kmc(1:length_A,1:length_A);
Delta_X=inv(Kmc)*R;
Delta_X(length_A+1)=0;
Delta_X(length_A+1)=Delta_X(Note_flag);
Delta_X(Note_flag)=0;
Delta_A=Delta_X(1:length_A);
domg=Delta_X(length_A+1);

%A00=[w0；A0(2:6，1)];
A=A+Delta_A;
omg0=omg0+domg;

if(n>60)
disp(‘迭代步数太多，可能不收敛’)
return
end
n=n+1;
end
I3=n-1;
disp([‘增量步=’ num2str(i),’ 迭代次数=’ num2str(I3),’ 本增量步弧长=’ num2str(delta_s),’ 已计算到频率比=’ num2str(omg0)])
Harmonic_A=[Harmonic_A;A’];
delta_s=delta_s*Nd/I3;
% [i I3 delta_s omg0]
%%%%%%%%%%%%%保存最新的四组解,用于预测
for p=1:3
X(:,p)=X(:,p+1);
end
X(1:length_A,4)=A;
X(length_A+1,4)=omg0;
p=0;
%%%%%%%%%%%%%
Frequency(i)=omg0;
% Amplitude11(i)=sqrt(A(2)^2+A(4)^2);
% Amplitude12(i)=sqrt(A(3)^2+A(5)^2);
Amplitude(i)=sqrt(A(2)^2+A(1)^2+A(3)^2+A(4)^2);
Result_A2(:,i)=A;
i=i+1;

end
Result_A= [Result_A1,zeros(length(A),Num_Incremental_step-Num_Pre_step)]+ Result_A2;

% figure(1)
% plot(Frequency,Amplitude11,’o-‘,’linewidth’,2,’color’,’r’);
% xlabel(‘Frequency’);
% ylabel(‘Amplitude’);
% figure(2)
% plot(Frequency,Amplitude12,’o-‘,’linewidth’,2,’color’,’b’);
% xlabel(‘Frequency’);
% ylabel(‘Amplitude’);
figure(3)
plot(Frequency,Amplitude,’o-‘,’linewidth’,2,’color’,’b’);
xlabel(‘Frequency’);
ylabel(‘Amplitude’);
grid on
toc
Warning: The maximum function count has been exceeded; May have singularity.Hello everyone, I am a beginner in Matlab programming and have been trying to solve the differential equation system for incremental harmonic balance (IHB) given below. I have a MATLAB program, but there seems to be a problem and the analysis results are not satisfactory. We would greatly appreciate any assistance or related work. Thank you in advance.

matlab code
clear all
close all
clc
tic
syms t
%========输入基本参数（已知条件）======
%========duffing 方程得参数===========

m=[1,0;0,miu];% m惯性项系数
k=[1+lambda_s,-lambda_s;-lambda_s,lambda_s+p0];% k一次项系数
f=[z;0];% f激励幅值
c=[2*xi,0;0,2*xi_b1];% c阻尼系数

%========设置谐波矩阵=================

%========刚度矩阵k====================

%========阻尼矩阵k====================

%========外激励矩阵f==================

%========非线性刚度矩阵===============

A0=[zeros(1,length(A1)),A(1:length(A2),:)’]’;
q=(S*A0)’;
kn3=[0;p1]*q.^2;
kn5=[0;p2]*q.^4;
kn7=[0;p3]*q.^6;
kn9=[0;p4]*q.^8;
kn11=[0;p5]*q.^10;

%=========带入公式，公式推导可见陈的书
Kmc=omg0^2*M+omg0*C+K+3*KN3+5*KN5+7*KN7+9*KN9+11*KN11;
R=-F-(omg0^2*M+omg0*C+K+KN3+KN5+KN7+KN9+KN11)*A;
Rmc=-(2*omg0*M+C)*A;
Delta_A=inv(Kmc)*R;

%====再一次计算
Kmc=omg0^2*M+omg0*C+K+3*KN3+5*KN5+7*KN7+9*KN9+11*KN11;
R=-F-(omg0^2*M+omg0*C+K+KN3+KN5+KN7+KN9+KN11)*A;
Rmc=-(2*omg0*M+C)*A;
tol=norm(R);
Delta_A=inv(Kmc)*R;

Kmc=omg0^2*M+omg0*C+K+3*KN3+5*KN5+7*KN7+9*KN9+11*KN11;
R=-F-(omg0^2*M+omg0*C+K+KN3+KN5+KN7+KN9+KN11)*A;
Rmc=-(2*omg0*M+C)*A;

%A00=[w0；A0(2:6，1)];
A=A+Delta_A;
omg0=omg0+domg;

%%%%%%%%%%%%%%%%%%%%%%%%%下面是每个增量步内的循环迭代
while tol>epsR

%====回代非线性刚度矩阵，计算下一次循环过程的非线性刚度矩阵

%%%%%%%%%%%%%%%%
Kmc=omg0^2*M+omg0*C+K+3*KN3+5*KN5+7*KN7+9*KN9+11*KN11;
R=-F-(omg0^2*M+omg0*C+K+KN3+KN5+KN7+KN9+KN11)*A;
Rmc=-(2*omg0*M+C)*A;
tol=norm(R);

%A00=[w0；A0(2:6，1)];
A=A+Delta_A;
omg0=omg0+domg;

end
Result_A= [Result_A1,zeros(length(A),Num_Incremental_step-Num_Pre_step)]+ Result_A2;

% figure(1)
% plot(Frequency,Amplitude11,’o-‘,’linewidth’,2,’color’,’r’);
% xlabel(‘Frequency’);
% ylabel(‘Amplitude’);
% figure(2)
% plot(Frequency,Amplitude12,’o-‘,’linewidth’,2,’color’,’b’);
% xlabel(‘Frequency’);
% ylabel(‘Amplitude’);
figure(3)
plot(Frequency,Amplitude,’o-‘,’linewidth’,2,’color’,’b’);
xlabel(‘Frequency’);
ylabel(‘Amplitude’);
grid on
toc
Warning: The maximum function count has been exceeded; May have singularity. Hello everyone, I am a beginner in Matlab programming and have been trying to solve the differential equation system for incremental harmonic balance (IHB) given below. I have a MATLAB program, but there seems to be a problem and the analysis results are not satisfactory. We would greatly appreciate any assistance or related work. Thank you in advance.

matlab code
clear all
close all
clc
tic
syms t
%========输入基本参数（已知条件）======
%========duffing 方程得参数===========

m=[1,0;0,miu];% m惯性项系数
k=[1+lambda_s,-lambda_s;-lambda_s,lambda_s+p0];% k一次项系数
f=[z;0];% f激励幅值
c=[2*xi,0;0,2*xi_b1];% c阻尼系数

%========设置谐波矩阵=================

%========刚度矩阵k====================

%========阻尼矩阵k====================

%========外激励矩阵f==================

%========非线性刚度矩阵===============

A0=[zeros(1,length(A1)),A(1:length(A2),:)’]’;
q=(S*A0)’;
kn3=[0;p1]*q.^2;
kn5=[0;p2]*q.^4;
kn7=[0;p3]*q.^6;
kn9=[0;p4]*q.^8;
kn11=[0;p5]*q.^10;

%=========带入公式，公式推导可见陈的书
Kmc=omg0^2*M+omg0*C+K+3*KN3+5*KN5+7*KN7+9*KN9+11*KN11;
R=-F-(omg0^2*M+omg0*C+K+KN3+KN5+KN7+KN9+KN11)*A;
Rmc=-(2*omg0*M+C)*A;
Delta_A=inv(Kmc)*R;

%====再一次计算
Kmc=omg0^2*M+omg0*C+K+3*KN3+5*KN5+7*KN7+9*KN9+11*KN11;
R=-F-(omg0^2*M+omg0*C+K+KN3+KN5+KN7+KN9+KN11)*A;
Rmc=-(2*omg0*M+C)*A;
tol=norm(R);
Delta_A=inv(Kmc)*R;

Kmc=omg0^2*M+omg0*C+K+3*KN3+5*KN5+7*KN7+9*KN9+11*KN11;
R=-F-(omg0^2*M+omg0*C+K+KN3+KN5+KN7+KN9+KN11)*A;
Rmc=-(2*omg0*M+C)*A;

%A00=[w0；A0(2:6，1)];
A=A+Delta_A;
omg0=omg0+domg;

%%%%%%%%%%%%%%%%%%%%%%%%%下面是每个增量步内的循环迭代
while tol>epsR

%====回代非线性刚度矩阵，计算下一次循环过程的非线性刚度矩阵

%%%%%%%%%%%%%%%%
Kmc=omg0^2*M+omg0*C+K+3*KN3+5*KN5+7*KN7+9*KN9+11*KN11;
R=-F-(omg0^2*M+omg0*C+K+KN3+KN5+KN7+KN9+KN11)*A;
Rmc=-(2*omg0*M+C)*A;
tol=norm(R);

%A00=[w0；A0(2:6，1)];
A=A+Delta_A;
omg0=omg0+domg;

end
Result_A= [Result_A1,zeros(length(A),Num_Incremental_step-Num_Pre_step)]+ Result_A2;

Cart

Cart