fm=inline(S'*m*S2);
There is no documented meaning for inline() of a symbolic expression. You should be using matlabFunction() instead of inline()
增量谐波平衡法（IHB）
28 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
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.
5 comentarios
Walter Roberson
el 14 de Dic. de 2024
inline is only documented for character vectors, not for symbolic expressions.
Walter Roberson
el 15 de Dic. de 2024
Editada: Walter Roberson
el 15 de Dic. de 2024
On iteration 37, Kmc and R both get larger, and after a few iterations that triggers "matrix is close to singular". A small number of iterations later, Kmc contains NaN, at which point everything blows up.
Note: I have recoded your inline() to matlabFunction. I have recoded your quadv() to integral('arrayvalued', true). I have recoded your inv(Kmc)*R to Kmc\R
format long g
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 = matlabFunction(S'*m*S2,'vars',t);
M = integral(fm,0,2*pi,'arrayvalued',true);%质量矩阵
fk = matlabFunction(S'*k*S,'vars',t);
K = integral(fk,0,2*pi,'arrayvalued',true);%线性刚度矩阵
fc = matlabFunction(S'*c*S1,'vars',t);
C = integral(fc,0,2*pi,'arrayvalued',true);%阻尼矩阵
fkn3 = matlabFunction(S'*kn3*S,'vars',t);
KN3 = integral(fkn3,0,2*pi,'arrayvalued',true);%非线性矩阵
fkn5 = matlabFunction(S'*kn5*S,'vars',t);
KN5 = integral(fkn5,0,2*pi,'arrayvalued',true);%非线性矩阵
fkn7 = matlabFunction(S'*kn7*S,'vars',t);
KN7 = integral(fkn7,0,2*pi,'arrayvalued',true);%非线性矩阵
fkn9 = matlabFunction(S'*kn9*S,'vars',t);
KN9 = integral(fkn9,0,2*pi,'arrayvalued',true);%非线性矩阵
fkn11 = matlabFunction(S'*kn11*S,'vars',t);
KN11 = integral(fkn11,0,2*pi,'arrayvalued',true);%非线性矩阵
ff = matlabFunction(S'*f*cos(t),'vars',t);
F = integral(ff,0,2*pi,'arrayvalued',true);%激励矩阵
%=========带入公式，公式推导可见陈的书
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 = 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 = matlabFunction(S'*kn3*S,'vars',t);
%disp('KN3');
KN3 = integral(fkn3,0,2*pi,'arrayvalued',true);%非线性矩阵
fkn5 = matlabFunction(S'*kn5*S,'vars',t);
%disp('KN5');
KN5 = integral(fkn5,0,2*pi,'arrayvalued',true);%非线性矩阵
fkn7 = matlabFunction(S'*kn7*S,'vars',t);
%disp('KN7');
KN7 = integral(fkn7,0,2*pi,'arrayvalued',true);%非线性矩阵
fkn9 = matlabFunction(S'*kn9*S,'vars',t);
%disp('KN9');
KN9 = integral(fkn9,0,2*pi,'arrayvalued',true);%非线性矩阵
fkn11 = matlabFunction(S'*kn11*S,'vars',t);
%disp('KN11');
KN11 = integral(fkn11,0,2*pi,'arrayvalued',true);%非线性矩阵
%====再一次计算
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 = 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 = matlabFunction(S'*kn3*S,'vars',t);
%disp('KN3b');
KN3 = integral(fkn3,0,2*pi,'arrayvalued',true);%非线性矩阵
fkn5 = matlabFunction(S'*kn5*S,'vars',t);
%disp('KN5b');
KN5 = integral(fkn5,0,2*pi,'arrayvalued',true);%非线性矩阵
fkn7 = matlabFunction(S'*kn7*S,'vars',t);
%disp('KN7b');
KN7 = integral(fkn7,0,2*pi,'arrayvalued',true);%非线性矩阵
fkn9 = matlabFunction(S'*kn9*S,'vars',t);
%disp('KN9b');
KN9 = integral(fkn9,0,2*pi,'arrayvalued',true);%非线性矩阵
fkn11 = matlabFunction(S'*kn11*S,'vars',t);
%disp('KN11b');
KN11 = integral(fkn11,0,2*pi,'arrayvalued',true);%非线性矩阵
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 = 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 = matlabFunction(S'*kn3*S,'vars',t);
%disp('KN3c');
KN3 = integral(fkn3,0,2*pi,'arrayvalued',true);%非线性矩阵
fkn5 = matlabFunction(S'*kn5*S,'vars',t);
%disp('KN5c');
KN5 = integral(fkn5,0,2*pi,'arrayvalued',true);%非线性矩阵
fkn7 = matlabFunction(S'*kn7*S,'vars',t);
%disp('KN7c');
KN7 = integral(fkn7,0,2*pi,'arrayvalued',true);%非线性矩阵
fkn9 = matlabFunction(S'*kn9*S,'vars',t);
%disp('KN9c');
KN9 = integral(fkn9,0,2*pi,'arrayvalued',true);%非线性矩阵
fkn11 = matlabFunction(S'*kn11*S,'vars',t);
%disp('KN11c');
KN11 = integral(fkn11,0,2*pi,'arrayvalued',true);%非线性矩阵
%%%%%%%%%%%%%%%%
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);
if i >= 34
[i, n, rcond(Kmc)]
Kmc
R
end
Delta_X = Kmc \ R;
if any(isnan(Delta_X(:)))
error('Delta_X contains nan on i = %d', i);
end
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
Respuestas (0)
Ver también
Categorías
Más información sobre Loops and Conditional Statements en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
También puede seleccionar uno de estos países/idiomas:
América
- América Latina (Español)
- Canada (English)
- United States (English)
Europa
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom(English)
Asia-Pacífico
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)