增量谐波平衡法（IHB）

46 views (last 30 days)
rong
rong on 14 Dec 2024
Commented: Lee on 4 Dec 2025 at 5:57
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.
  8 Comments
Show 6 older comments
Walter Roberson
Walter Roberson on 3 Dec 2025 at 18:42
Moved: Walter Roberson on 3 Dec 2025 at 18:44
Approximate translation:
Hello, I apologize for bothering you, but I'd like to ask if your paper has been published? I'm a beginner in IHBM and am learning the IHBM code. I'd like to combine your code with the paper and reproduce the result. Could you please send me the paper's .doi file? Thank you.
Lee
Lee 24 minutes ago
sry Mr Roberson, I didn't mean to bother you, its my fault that i accidently comment at the wrong place. I'm trying to cntect with the original poster, now im working on my paper using IHBM, and im a new beginner with matlab. I think you are a helpful and kind person. I've seen you help many people solve code problems, and your programming skills are extremely high.

Sign in to comment.

Sign in to answer this question.

Answers (0)

Categories

Find more on Programming in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!