增量谐波平衡法（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.

[Walter Roberson]

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

There is no documented meaning for inline() of a symbolic expression. You should be using matlabFunction() instead of inline()

[rong]

In MATLAB, inline is a command used to create inline functions. Inline functions refer to simple functions defined directly in the command line without the need to write a complete function file.

[PCG]

inline function is not recommanded in the Documentation. It suggests to use Anonymous Functions for clarifty.

[Walter Roberson]

inline is only documented for character vectors, not for symbolic expressions.

[Walter Roberson]

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

增量步=1 迭代次数=6 本增量步弧长=0.02 已计算到频率=0.005
增量步=2 迭代次数=2 本增量步弧长=0.02 已计算到频率=0.0233
增量步=3 迭代次数=2 本增量步弧长=0.02 已计算到频率=0.0416
增量步=4 迭代次数=2 本增量步弧长=0.02 已计算到频率=0.0599

% 以下是结合了弧长法的增量谐波平衡法
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

增量步=5 迭代次数=1 本增量步弧长=0.02 已计算到频率比=0.052463
增量步=6 迭代次数=1 本增量步弧长=0.02 已计算到频率比=0.055809
增量步=7 迭代次数=1 本增量步弧长=0.02 已计算到频率比=0.060259
增量步=8 迭代次数=1 本增量步弧长=0.02 已计算到频率比=0.064067
增量步=9 迭代次数=1 本增量步弧长=0.02 已计算到频率比=0.06829
增量步=10 迭代次数=1 本增量步弧长=0.02 已计算到频率比=0.072234
增量步=11 迭代次数=1 本增量步弧长=0.02 已计算到频率比=0.076008
增量步=12 迭代次数=1 本增量步弧长=0.02 已计算到频率比=0.080285
增量步=13 迭代次数=1 本增量步弧长=0.02 已计算到频率比=0.084719
增量步=14 迭代次数=1 本增量步弧长=0.02 已计算到频率比=0.088785
增量步=15 迭代次数=1 本增量步弧长=0.02 已计算到频率比=0.0924
增量步=16 迭代次数=1 本增量步弧长=0.02 已计算到频率比=0.096593
增量步=17 迭代次数=1 本增量步弧长=0.02 已计算到频率比=0.10149
增量步=18 迭代次数=1 本增量步弧长=0.02 已计算到频率比=0.10576
增量步=19 迭代次数=1 本增量步弧长=0.02 已计算到频率比=0.1142
增量步=20 迭代次数=1 本增量步弧长=0.02 已计算到频率比=0.11878
增量步=21 迭代次数=1 本增量步弧长=0.02 已计算到频率比=0.12383
增量步=22 迭代次数=1 本增量步弧长=0.02 已计算到频率比=0.12804
增量步=23 迭代次数=1 本增量步弧长=0.02 已计算到频率比=0.1354
增量步=24 迭代次数=1 本增量步弧长=0.02 已计算到频率比=0.14082
增量步=25 迭代次数=1 本增量步弧长=0.02 已计算到频率比=0.14616
增量步=26 迭代次数=1 本增量步弧长=0.02 已计算到频率比=0.15217
增量步=27 迭代次数=1 本增量步弧长=0.02 已计算到频率比=0.15803
增量步=28 迭代次数=1 本增量步弧长=0.02 已计算到频率比=0.16355
增量步=29 迭代次数=1 本增量步弧长=0.02 已计算到频率比=0.16913
增量步=30 迭代次数=1 本增量步弧长=0.02 已计算到频率比=0.17496
增量步=31 迭代次数=1 本增量步弧长=0.02 已计算到频率比=0.18096
增量步=32 迭代次数=1 本增量步弧长=0.02 已计算到频率比=0.18682
增量步=33 迭代次数=1 本增量步弧长=0.02 已计算到频率比=0.19293

ans = 1×3

34 1 0.000832472831677305

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Kmc = 8×8

3.07442496358979 0.0625176408410071 -3.53788592271287e-16 9.14550841427191e-17 0.089942896512981 -2.60208521396521e-18 8.23993651088983e-18 2.60208521396521e-18 -0.062517640841007 3.07442496358979 1.09771475532078e-16 2.53593591977693e-16 0.0444535719428989 -0.0572423098724829 -2.49366499671666e-18 -5.63785129692462e-18 -4.05789216512659e-16 8.78076990560145e-17 2.70119496397232 0.125035281682014 0.000315208063140634 -2.49366499671666e-18 -0.0572423098724828 4.22838847269347e-18 8.64033842209936e-17 1.8544682515783e-16 -0.125035281682014 2.70119496397232 0.000389817640602169 -5.63785129692462e-18 4.22838847269347e-18 -0.0572423098724828 -0.0572423098724828 -2.60208521396521e-18 8.23993651088983e-18 2.60208521396521e-18 -0.00501613929517051 0.00240913236343865 1.00330631245916e-05 9.68067493569551e-06 -2.60208521396521e-18 -0.0572423098724829 -2.49366499671666e-18 -5.63785129692462e-18 0.00510180010813734 0.0040911150198678 2.17121421697791e-06 2.06186839309913e-06 8.23993651088983e-18 -2.49366499671666e-18 -0.0572423098724828 4.22838847269347e-18 -3.87249992969476e-05 2.17121421697792e-06 0.00663885403834666 0.00249629880138431 2.60208521396521e-18 -5.63785129692462e-18 4.22838847269347e-18 -0.0572423098724828 4.44674274215452e-05 2.06186839309902e-06 -0.00250511246589626 0.00664658111122844

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

R = 8×1

1.0e+00 * -9.6934948429217e-08 -6.02702957085536e-08 3.99559977882295e-07 4.30458722538031e-08 6.59038878482197e-07 -2.31314160354037e-06 1.71628362248177e-05 -6.08609016267479e-07

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

增量步=34 迭代次数=1 本增量步弧长=0.02 已计算到频率比=0.19869

ans = 1×3

35 1 0.000851476172632113

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Kmc = 8×8

3.06530896050632 0.0647676084126644 -3.49437352334983e-16 9.24669484818719e-17 0.0930915549528194 -2.60208521396521e-18 8.23993651088983e-18 2.60208521396521e-18 -0.0647676084126644 3.06530896050632 1.11584049114765e-16 2.58921804709046e-16 0.0456350268022624 -0.0572423098724829 -2.49366499671666e-18 -5.63785129692462e-18 -4.04884658502844e-16 8.8092149691482e-17 2.66473095163847 0.129535216825329 -0.000610436646085905 -2.49366499671666e-18 -0.0572423098724828 4.22838847269347e-18 8.6695563122308e-17 1.86547483898593e-16 -0.129535216825329 2.66473095163847 0.000220282228900084 -5.63785129692462e-18 4.22838847269347e-18 -0.0572423098724828 -0.0572423098724828 -2.60208521396521e-18 8.23993651088983e-18 2.60208521396521e-18 -0.00510636642178352 0.00247835550761723 -1.55440607772182e-05 7.87878341508535e-06 -2.60208521396521e-18 -0.0572423098724829 -2.49366499671666e-18 -5.63785129692462e-18 0.00499208848217003 0.00409957304811971 -3.47827887767147e-06 1.75294957507187e-06 8.23993651088983e-18 -2.49366499671666e-18 -0.0572423098724828 4.22838847269347e-18 -3.57598305587269e-05 -3.47827887767147e-06 0.00663287270508995 0.00258622028787487 2.60208521396521e-18 -5.63785129692462e-18 4.22838847269347e-18 -0.0572423098724828 -6.30328201507478e-05 1.75294957507175e-06 -0.00259518838513829 0.0066403988874158

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

R = 8×1

1.0e+00 * -1.5018037231862e-06 -8.21856292510126e-07 -8.06555521311345e-07 1.83955543158138e-06 3.14048135617079e-06 -5.9290581340229e-06 -5.96159554563798e-06 2.05747691296829e-05

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

增量步=35 迭代次数=1 本增量步弧长=0.02 已计算到频率比=0.20519

ans = 1×3

36 1 0.00122057999020293

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Kmc = 8×8

3.05850015413942 0.0663984040479582 -3.4619232271063e-16 9.32315025736719e-17 -0.0572423098724828 0.0953029588478128 8.23993651088983e-18 2.60208521396521e-18 -0.0663984040479582 3.05850015413942 1.12940532670303e-16 2.6288652111866e-16 -2.60208521396521e-18 0.0466379724923455 -2.49366499671666e-18 -5.63785129692462e-18 -4.04206223686534e-16 8.83090006327851e-17 2.63749572617087 0.132796808095916 8.23993651088983e-18 4.69945692750575e-05 -0.0572423098724828 4.22838847269347e-18 8.69151234282402e-17 1.87370946341022e-16 -0.132796808095917 2.63749572617087 2.60208521396521e-18 -0.00080147251661237 4.22838847269347e-18 -0.0572423098724828 -0.0572423098724828 -2.60208521396521e-18 8.23993651088983e-18 2.60208521396521e-18 0.0091637250731164 -0.00520474775030154 -1.21235862584795e-06 -2.10019520815407e-05 -2.60208521396521e-18 -0.0572423098724829 -2.49366499671666e-18 -5.63785129692462e-18 -0.000120784959741116 0.0049035132482177 -2.47780346867045e-07 -4.72065702420196e-06 8.23993651088983e-18 -2.49366499671666e-18 -0.0572423098724828 4.22838847269347e-18 -1.21235862584816e-06 8.64133067331129e-05 0.0066284818363454 0.00265139214062799 2.60208521396521e-18 -5.63785129692462e-18 4.22838847269347e-18 -0.0572423098724828 -2.10019520815405e-05 -9.83256832952439e-06 -0.00266048018320867 0.00663584707071047

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

R = 8×1

1.0e+00 * -1.75554898396157e-09 -1.71602290305338e-09 -1.83320543896327e-07 -1.94021380675135e-07 1.20134007978009e-05 1.31577945503982e-06 -2.03299044197231e-05 -1.6131447158779e-05

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

增量步=36 迭代次数=1 本增量步弧长=0.02 已计算到频率比=0.21161

ans = 1×3

37 1 4.53638558199135e-06

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Kmc = 8×8

-51.5425215506293 1.31139255553861 2.45001996593025e-14 8.3151243536515e-15 -0.0572423098724828 -2.60208521396521e-18 3.12562135262721 2.60208521396521e-18 -1.31139255553862 -51.5425215506294 1.16237784849135e-14 2.84957570503432e-14 -2.60208521396521e-18 -0.0572423098724829 5.72307290192964 -5.63785129692462e-18 5.7059913609495e-15 2.87267425453717e-15 -215.766591092904 2.62278511107723 8.23993651088983e-18 -2.49366499671666e-18 0.356137843178012 4.22838847269347e-18 2.16408630547226e-15 7.11754907678935e-15 -2.62278511107723 -215.766591092904 2.60208521396521e-18 -5.63785129692462e-18 0.379086158863104 -0.0572423098724828 -0.0572423098724828 -2.60208521396521e-18 8.23993651088983e-18 2.60208521396521e-18 0.0409440325015182 0.0489235940343938 -0.0518265029820299 0.000691140347894789 -2.60208521396521e-18 -0.0572423098724829 -2.49366499671666e-18 -5.63785129692462e-18 -0.00353210818715083 0.0685935265978579 -0.018599858720303 0.00131622217025467 8.23993651088983e-18 -2.49366499671666e-18 -0.0572423098724828 4.22838847269347e-18 0.000716822032502853 0.00125088291253864 -0.00205216097088551 0.0529532372386194 2.60208521396521e-18 -5.63785129692462e-18 4.22838847269347e-18 -0.0572423098724828 0.000691140347894792 0.00131622217025467 0.000393825722812202 0.0528747013160192

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

R = 8×1

-5.99165512015018 -11.7956704529291 -0.728678699745335 -0.793626958588706 0.22217022090205 0.241344812046751 0.0114782792903054 0.00461964065644804

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

ans = 1×3

37 2 2.4087535230428e-06

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Kmc = 8×8

-2.01891059507366 0.404870733876842 1.98643632809383e-15 8.18444278069471e-16 -0.0572423098724828 -2.60208521396521e-18 1.37090155766482 2.60208521396521e-18 -0.404870733876842 -2.01891059507366 1.17143366998967e-15 2.95531162033226e-15 -2.60208521396521e-18 -0.0572423098724829 2.56432561617624 -5.63785129692462e-18 1.51383553841923e-16 3.27555105179746e-16 -17.6721472706815 0.809741467753684 8.23993651088983e-18 -2.49366499671666e-18 0.152945189763104 4.22838847269347e-18 2.74117270485322e-16 8.25667556816235e-16 -0.809741467753684 -17.6721472706815 2.60208521396521e-18 -5.63785129692462e-18 0.164447106075729 -0.0572423098724828 -0.0572423098724828 -2.60208521396521e-18 8.23993651088983e-18 2.60208521396521e-18 0.0822766060771253 0.0524696212610081 -0.034871455075073 0.00119461076654952 -2.60208521396521e-18 -0.0572423098724829 -2.49366499671666e-18 -5.63785129692462e-18 0.0362747919059344 0.132078621306665 -0.0114852745660001 0.00235596707491647 8.23993651088983e-18 -2.49366499671666e-18 -0.0572423098724828 4.22838847269347e-18 0.00137505875710394 0.00227753341655217 -0.00129463623056605 0.0180141062433657 2.60208521396521e-18 -5.63785129692462e-18 4.22838847269347e-18 -0.0572423098724828 0.00119461076654952 0.00235596707491647 0.000254475148581731 0.105919077745201

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

R = 8×1

-0.320474096745202 -1.01485634565629 -0.082128062777041 -0.0978973569078619 0.0661070434519298 0.227569661121666 0.00572755553165938 0.00717253003544974

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

ans = 1×3

37 3 2.97537640308107e-15

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Kmc = 8×8

1.0e+00 * -113.595354640814 -1.91551499091699 5.22762044404616e-14 1.84793177474687e-14 -0.0572423098724828 -2.60208521396521e-18 -247.158281718397 2.60208521396521e-18 1.91551499091698 -113.595354640814 2.49538887723875e-14 5.91848812218929e-14 -2.60208521396521e-18 -0.0572423098724829 -301.80406631147 -5.63785129692462e-18 1.2912812868005e-14 6.44714128713069e-15 -463.977923453645 -3.83102998183397 8.23993651088983e-18 -2.49366499671666e-18 -10.2900128536383 4.22838847269347e-18 4.64889556179024e-15 1.51216964455504e-14 3.83102998183396 -463.977923453645 2.60208521396521e-18 -5.63785129692462e-18 -10.0362440847866 -0.0572423098724828 -0.0572423098724828 -2.60208521396521e-18 8.23993651088983e-18 2.60208521396521e-18 12245931221.6238 12118763810.0923 -0.0394583106012512 934201079.523372 -2.60208521396521e-18 -0.0572423098724829 -2.49366499671666e-18 -5.63785129692462e-18 12118763810.1689 17514578029.6439 0.0255843689250553 768359441.111766 8.23993651088983e-18 -2.49366499671666e-18 -0.0572423098724828 4.22838847269347e-18 297222091.130269 227183068.691815 0.000367331707233437 -2940262641.79326 2.60208521396521e-18 -5.63785129692462e-18 4.22838847269347e-18 -0.0572423098724828 934201079.523372 768359441.111766 0.000269039530171873 21330710092.8082

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

R = 8×1

1.0e+00 * -741.10443731277 -888.950722802049 -31.2835602444747 -30.2565794947039 11537519600.9503 14428812042.6315 265374792.984389 730073586.153924

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

ans = 1×3

37 4 1.83275623637216e-106

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Kmc = 8×8

1.0e+00 * -1.80038750887787e+21 7520579656.59886 815192.59622369 278358.447489432 -0.0572423098724828 -2.60208521396521e-18 7.67098711261595e+21 2.60208521396521e-18 -7520808805.24938 -1.80038750887787e+21 381748.727986927 918582.876721774 -2.60208521396521e-18 -0.0572423098724829 9.69601650350976e+21 -5.63785129692462e-18 203798.14905531 95437.1819961678 -7.20155003551147e+21 15041316883.9617 8.23993651088983e-18 -2.49366499671666e-18 3.16835116408602e+20 4.22838847269347e-18 69589.6118724891 229645.71917967 -15041460039.7347 -7.20155003551147e+21 2.60208521396521e-18 -5.63785129692462e-18 3.14167604511445e+20 -0.0572423098724828 -0.0572423098724828 -2.60208521396521e-18 8.23993651088983e-18 2.60208521396521e-18 1.42702948587817e+109 1.4306050087834e+109 5679793.92844287 1.06678705379814e+108 -2.60208521396521e-18 -0.0572423098724829 -2.49366499671666e-18 -5.63785129692462e-18 1.4306050087834e+109 2.10122333143589e+109 8858389.71981873 8.84433720186986e+107 8.23993651088983e-18 -2.49366499671666e-18 -0.0572423098724828 4.22838847269347e-18 3.28352809884526e+107 2.62784941928752e+107 105791.381667565 -3.74424579532039e+108 2.60208521396521e-18 -5.63785129692462e-18 4.22838847269347e-18 -0.0572423098724828 1.06678705379814e+108 8.84433720186986e+107 -183952.951635323 2.51524725913353e+109

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

R = 8×1

1.0e+00 * -9.18183146011467e+31 -1.16057018559818e+32 -3.79237586611222e+30 -3.76044693142364e+30 1.25556908531097e+109 1.61588228414682e+109 2.8001145285291e+107 7.63834470822325e+107

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.832756e-106.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

ans = 1×3

37 5 NaN

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Kmc = 8×8

1.0e+00 * -1.42401001894909e+203 -9.06221501108151e+186 6.44773649378467e+187 2.20166611982891e+187 -0.0572423098724828 -2.60208521396521e-18 -1.21346792840842e+204 2.60208521396521e-18 -9.06221501108151e+186 -1.42401001894909e+203 3.01942782147965e+187 7.26549819543541e+187 -2.60208521396521e-18 -0.0572423098724829 -1.53380586976563e+204 -5.63785129692462e-18 1.61193412344617e+187 7.54856955369913e+186 -5.69604007579636e+203 -5.66142716527435e+186 8.23993651088983e-18 -2.49366499671666e-18 -5.01199189504365e+202 4.22838847269347e-18 5.50416529957228e+186 1.81637454885885e+187 -5.66142716527435e+186 -5.69604007579636e+203 2.60208521396521e-18 -5.63785129692462e-18 -4.96979471637668e+202 -0.0572423098724828 -0.0572423098724828 -2.60208521396521e-18 8.23993651088983e-18 2.60208521396521e-18 NaN NaN -4.61350289930706e+97 NaN -2.60208521396521e-18 -0.0572423098724829 -2.49366499671666e-18 -5.63785129692462e-18 NaN NaN -7.02710768236552e+97 NaN 8.23993651088983e-18 -2.49366499671666e-18 -0.0572423098724828 4.22838847269347e-18 NaN NaN -9.4085742042866e+95 NaN 2.60208521396521e-18 -5.63785129692462e-18 4.22838847269347e-18 -0.0572423098724828 NaN NaN 1.23772811895241e+96 NaN

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

R = 8×1

1.0e+00 * -1.29175423182433e+305 -1.63275862236047e+305 -5.3353381566363e+303 -5.29041864715606e+303 NaN NaN NaN NaN

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN.

Delta_X contains nan on i = 37

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