声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 9066|回复: 16

[混合编程] 求助:增量谐波平衡法编程 修改

[复制链接]
发表于 2010-8-28 09:56 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
本帖最后由 ChaChing 于 2010-8-28 23:59 编辑

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

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)

本帖被以下淘专辑推荐:

回复
分享到:

使用道具 举报

发表于 2010-8-29 00:00 | 显示全部楼层
请先看下:@)
建议提问的网友分清 编程问题 和 专业问题
http://forum.vibunion.com/forum/ ... p;extra=&page=1

评分

1

查看全部评分

 楼主| 发表于 2010-8-29 13:51 | 显示全部楼层

求助:用IHB法编程

附未完成之程序
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的同学帮忙看看
112.jpg
发表于 2010-11-11 15:36 | 显示全部楼层
我也在做有关IHB的计算。看了你的程序很有启发,但是你在论坛里留的有关K3矩阵的计算方法好像不清楚!!程序中两处也不一样??能说的详细一下吗??我的邮箱chaofeng_ma@163.com
发表于 2010-11-18 10:24 | 显示全部楼层
把你的程序运行了一遍,是不收敛,建议在计算delta_a时,把Rmc项也考虑上,我计算了一次,仅迭代了30次,得到的R的范数基本上趋于一个小于1的数!!

评分

1

查看全部评分

发表于 2012-12-3 17:23 | 显示全部楼层

请问有无解决  可不可以分享一下
发表于 2012-12-24 15:47 | 显示全部楼层
不知道是否解决了,这个问题
发表于 2014-4-13 13:30 | 显示全部楼层
现在这个问题解决了吗?
发表于 2014-4-13 13:31 | 显示全部楼层
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)
发表于 2014-10-27 22:04 | 显示全部楼层
不知道各位大神是否已经有成功的案例没?程序问题好头疼啊!
发表于 2014-11-25 17:15 | 显示全部楼层
还有继续做的吗?求赐教啊
发表于 2019-3-4 17:40 | 显示全部楼层
这个只能得到稳定解吧
发表于 2019-9-16 11:20 | 显示全部楼层
alittlehug 发表于 2019-3-4 17:40
这个只能得到稳定解吧

可以互相学习下吗
发表于 2019-10-21 11:24 | 显示全部楼层
你好,请问这个问题解决了么?
发表于 2019-12-10 21:48 | 显示全部楼层
增量谐波平衡法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

点评

大家都是不懂,来学习交流的,你在这边卖程序,一个程序卖5000,你良心不会痛吗?不怕遭报应毕不了业吗  发表于 2020-6-13 09:29
回复 支持 0 反对 1

使用道具 举报

您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-6-14 20:25 , Processed in 0.060783 second(s), 23 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表