求助：增量谐波平衡法编程 修改

zhangwenjing

怎么迭代啊，颤振的计算。求帮助, 请各位同学不吝赐教

clear all; clc
syms t
w0=1; ww=0.1; Q=10; epsR=0.01;
m=[1 0.25;0.25 0.5]; c=[0.1 0;0 0.1]; k=[0.2 0.1*Q;0 0.5-0.04*Q];
Cs=[1 cos(t) cos(2*t) sin(t) sin(2*t) sin(3*t)];
S=[Cs zeros(1,6);zeros(1,6) Cs];
A1=[0.1 0.1 0.1 0.1 0.1 0.1]'; A2=[0.1 0.1 0.1 0.1 0.1 0.1]';
T1=[zeros(6,6) eye(6,6)]; T2=[eye(6,6) zeros(6,6)];
S2=diff(S,t,2); fm=inline(S'*m*S2); M=quadv(fm,0,2*pi);
S1=diff(S,t,1); fc=inline(S'*c*S1); C=quadv(fc,0,2*pi);
fk=inline(S'*k*S); K=quadv(fk,0,2*pi);
A0=[A1;A2]; %X0=S*A0;
k3=[10*Cs*A1 0;0 20*Cs*A2];
fk3=inline(S'*k3*S);
K3=quadv(fk3,0,2*pi);
Kmc=w0^2*M+w0*C+K+K3;
R=-(w0^2*M+w0*C+K+K3)*A0;
Rmc=(2*w0*M+C)*A0;
AA=inv(Kmc)*(R-Rmc*ww);
A01=AA+A0; A10=T1*A01; A20=T2*A01;
n=1; tol=1;
while tol>epsR
    A0=A01; A1=A10; A2=A20;
    k3=[10*Cs*A1 0;0 20*Cs*A2];
    fk3=inline(S'*k3*S); K3=quadv(fk3,0,2*pi);
    Kmc=w0^2*M+w0*C+K+K3;
    R=-(w0^2*M+w0*C+K+K3)*A0;
    Rmc=(2*w0*M+C)*A0;
   AA=inv(Kmc)*(R-Rmc*ww);
   tol=norm(R)
   A01=AA+A0; A10=T1*A01; A20=T2*A01;
   n=n+1; oo(n)=tol;
   if(n>500), disp('迭代步数太多，可能不收敛'); return; end
end
fA=norm(AA), plot(1:n,oo)

ChaChing

请先看下:@)
建议提问的网友分清 编程问题 和 专业问题
http://forum.vibunion.com/forum/ ... p;extra=&page=1

zhangwenjing

附未完成之程序
clear all
clc
syms t
w0=0.2;
Q=10;
epsR=0.01;
m=[1 0.25;0.25 0.5];
c=[0.1 0;0 0.1];
k=[0.2 0.1*Q;0 0.5-0.04*Q];
Cs=[1 cos(t) cos(2*t) sin(t) sin(2*t) sin(3*t)];
S=[Cs zeros(1,6);zeros(1,6) Cs];
A1=[0.1 0.1 0.1 0.1 0.1 0.1]';
A2=[0.1 0.1 0.1 0.1 0.1 0.1]';
T1=[zeros(6,6) eye(6,6)];
T2=[eye(6,6) zeros(6,6)];
%A0=[A1;A2];
%X0=S*A0;
S2=diff(S,t,2);
fm=inline(S'*m*S2);
M=quadv(fm,0,2*pi);
S1=diff(S,t,1);
fc=inline(S'*c*S1);
C=quadv(fc,0,2*pi);
fk=inline(S'*k*S);
K=quadv(fk,0,2*pi);
A0=[A1;A2];
%X0=S*A0;
k3=[10*(Cs*A1).^2 0;0 20*(Cs*A2).^2];
fk3=inline(S'*k3*S);
K3=quadv(fk3,0,2*pi);
%%%%%
Kmc=w0^2*M+w0*C+K+3*K3;
R=-(w0^2*M+w0*C+K+K3)*A0;
Rmc=(2*w0*M+C)*A0;
%AA=inv(Kmc)*(R-Rmc*ww);
%AA首元素已知a1=0.0,求ww
a1=0.0;
Kmc11=Rmc(:,1);
Kmc=[Kmc11 Kmc(:,2:size(Kmc,2))];
AA=inv(Kmc)*R   ;  %drtA1(1)
ww=AA(1)     ;     %drtW(1)
%%%A00=[w0;A0(2:12,1)] ;            %A1(0)
A01=A0+[a1;AA(2:12,1)]  ;      %A(1)+drtA(1)
%%%Aw0=AA+A00;              %A1(0)+drtA1(1)=A1(1)
A10=T2*A01;                  
A20=T1*A01;
w01=w0+ww;         %W+drtW(1)
n=1;
tol=1;
   
while tol>epsR
    A0=A01;
    %A0=[a1 A00(2:12,1)]
    A1=A10;
    A2=A20;
    w0=w01;

    k3=[10*(Cs*A1).^2 0;0 20*(Cs*A2).^2];
fk3=inline(S'*k3*S);
K3=quadv(fk3,0,2*pi);
    Kmc=w0^2*M+w0*C+K+3*K3;
R=-(w0^2*M+w0*C+K+K3)*A0;
Rmc=(2*w0*M+C)*A0;
tol=norm(R)
Kmc11=Rmc(:,1);
Kmc=[Kmc11 Kmc(:,2:size(Kmc,2))];
AA=inv(Kmc)*R;
ww=AA(1);
%%%A00=[w0;A0(2:12,1)];
A01=A0+[a1;AA(2:12,1)];
A10=T2*A01;
A20=T1*A01;
w01=w0+ww;
n=n+1;
oo(n)=tol;
if(n>1000)
    disp('迭代步数太多，可能不收敛')
   return;
end
end
遇到的问题是：迭代不收敛，请各位熟悉IHB的同学帮忙看看

xautmcf

我也在做有关IHB的计算。看了你的程序很有启发，但是你在论坛里留的有关K3矩阵的计算方法好像不清楚！！程序中两处也不一样？？能说的详细一下吗？？我的邮箱chaofeng_ma@163.com

xautmcf

把你的程序运行了一遍，是不收敛，建议在计算delta_a时，把Rmc项也考虑上，我计算了一次，仅迭代了30次，得到的R的范数基本上趋于一个小于1的数！！

皮卡丘

xautmcf 发表于 2010-11-18 10:24

                               
登录/注册后可看大图

把你的程序运行了一遍，是不收敛，建议在计算delta_a时，把Rmc项也考虑上，我计算了一次，仅迭代了30次，得 ...

请问有无解决  可不可以分享一下

Angelaxiaosu

不知道是否解决了，这个问题

twtp100

现在这个问题解决了吗？

twtp100

zhangwenjing 发表于 2010-8-29 13:51
附未完成之程序
clear all
clc

a1=0.0;
Kmc11=Rmc(:,1);
Kmc=[Kmc11 Kmc(:,2:size(Kmc,2))];
AA=inv(Kmc)*R   ;  %drtA1(1)
这一步是主要是干什么？
ww=AA(1)     ;     %drtW(1)
%%%A00=[w0;A0(2:12,1)] ;            %A1(0)
A01=A0+[a1;AA(2:12,1)]  ;      %A(1)+drtA(1)
%%%Aw0=AA+A00;              %A1(0)+drtA1(1)=A1(1)
A10=T2*A01;                  
A20=T1*A01;
w01=w0+ww;         %W+drtW(1)

thdl520

不知道各位大神是否已经有成功的案例没？程序问题好头疼啊！

薛永宁

还有继续做的吗？求赐教啊

alittlehug

这个只能得到稳定解吧

zhouyuan1996

alittlehug 发表于 2019-3-4 17:40
这个只能得到稳定解吧

可以互相学习下吗

黄小文

你好，请问这个问题解决了么？

博博士

增量谐波平衡法matlab程序，单自由度，多自由度均可以，弧长增量法求解多值区间，不稳定解，获得复杂幅频曲线，并判断解的稳定性和分岔。可以参考：https://item.taobao.com/item.htm ... ;abbucket=19#detail

补充内容 (2019-12-14 20:44):
参考TB店铺（增量谐波平衡法）：https://item.taobao.com/item.htm ... amp;id=610346631802