杭州锐达数字技术有限公司
查看: 4290|回复: 12
打印 上一主题 下一主题

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

[复制链接]
跳转到指定楼层
楼主
发表于 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的同学帮忙看看
地板
发表于 2010-11-11 15:36 | 只看该作者
我也在做有关IHB的计算。看了你的程序很有启发,但是你在论坛里留的有关K3矩阵的计算方法好像不清楚!!程序中两处也不一样??能说的详细一下吗??我的邮箱chaofeng_ma@163.com
5
发表于 2010-11-18 10:24 | 只看该作者
把你的程序运行了一遍,是不收敛,建议在计算delta_a时,把Rmc项也考虑上,我计算了一次,仅迭代了30次,得到的R的范数基本上趋于一个小于1的数!!

评分

1

查看全部评分

6
发表于 2012-12-3 17:23 | 只看该作者

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

可以互相学习下吗
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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