我要加入 登录
声振论坛 返回首页

zhangwenjing的个人空间 http://home.vibunion.com/?124723 [收藏] [复制] [分享] [RSS]

日志

IHB FDB

热度 3已有 3833 次阅读2010-9-1 13:57 |

clear all
clc
tic
%%%定义各参数
syms t
w0=3;
epsR=0.001;
m=[1 0;0 1];
epsilon=0.24;
r=0.4;
delta=0.56;
k=[1+epsilon*r -epsilon*r;-epsilon*r 1+epsilon*delta];
Cs=[cos(t) cos(3*t) sin(t) sin(3*t)];
Cs1=diff(Cs,t,1);
S=[cos(t) cos(3*t) sin(t) sin(3*t) 0 0 0 0;0 0 0 0 cos(t) cos(3*t) sin(t) sin(3*t)];
A1=[1 1 1 1]';
A2=[1 1 1 1]';
A0=[A1;A2];
T1=[eye(4,4) zeros(4,4)];
T2=[zeros(4,4) eye(4,4)];
S2=diff(S,t,2);
fm=inline(S'*m*S2);
M=quadv(fm,0,2*pi);
fk=inline(S'*k*S);
K=quadv(fk,0,2*pi);
S1=diff(S,t,1);
c=[1 0;0 1];
fc=inline(S'*c*S1);
C=quadv(fc,0,2*pi);
c3=diag(S*A0).^2;
fc3=inline(S'*c3*S1);
C3=quadv(fc3,0,2*pi);
k2=diag(S*A0).*diag(S1*A0)
fk2=inline(S'*k2*S);
K2=quadv(fk2,0,2*pi);
%%%%%带入推导出的公式
Kmc=w0^2*M+epsilon*w0*(C3-C)+K+2*epsilon*w0*K2;
R=-(w0^2*M+epsilon*w0*(C3-C)+K)*A0;
Rmc=-(2*w0*M+epsilon*(C3-C))*A0;
%AA首元素已知a1=0.0,求ww
 a1=0.0;
 %%%变换矩阵,使ww变量代替a1
 Kmc11=-Rmc(:,1);
 Kmcr=[Kmc11 Kmc(:,2:size(Kmc,2))];
%求未知变量
 AA=inv(Kmcr)*R;     %drtA1(1)
 ww=AA(1);     %drtW(1)
        
   %%赋予新变量新值        
 A01=A0+[a1;AA(2:length(A0),1)];        %A(1)+drtA(1)
% Aw0=AA+A00;              %A1(0)+drtA1(1)=A1(1)
 w01=w0+ww;         %W+drtW(1)
 n=1;
 tol=1;
%    
 while tol>epsR
     A0=A01;
%    
     w0=w01

c3=diag(S*A0).^2;
fc3=inline(S'*c3*S1);
C3=quadv(fc3,0,2*pi);
k2=diag(S*A0).*diag(S1*A0)
fk2=inline(S'*k2*S);
K2=quadv(fk2,0,2*pi);
%%%%%带入推导出的公式
Kmc=w0^2*M+epsilon*w0*(C3-C)+K+2*epsilon*w0*K2;
R=-(w0^2*M+epsilon*w0*(C3-C)+K)*A0;
Rmc=-(2*w0*M+epsilon*(C3-C))*A0;
%%%%%
 tol=norm(R)
 if(n>1000)
     disp('迭代步数太多,可能不收敛')
    return;
 end
 
 
 Kmc11=-Rmc(:,1);
 Kmcr=[Kmc11 Kmc(:,2:size(Kmc,2))];
%
 AA=inv(Kmcr)*R;
ww=AA(1);
% A00=[w0;A0(2:6,1)];
 A01=A0+[a1;AA(2:length(A0),1)];
 w01=w0+ww;
 n=n+1;
 end
 
  X0=S*A0;
  dX0=S1*A0;
 %%绘范德波图
tt=0:.1:10;
xo1=subs(X0(1),tt);
xo2=subs(X0(2),tt);
dxo1=subs(dX0(1),tt);
dxo2=subs(dX0(2),tt);
figure(1)
plot(xo1,dxo1,'b','linewidth',2)
hold on
plot(xo2,dxo2,'b','linewidth',2)
axis([-3 3 -3 3])
title('范德波极限环')
xlabel('x0')
ylabel('dx0')
toc
 

发表评论 评论 (6 个评论)

回复 zhangwenjing 2010-9-2 15:22
增量谐波平衡法,MATLAB编程
回复 大魁天下vip 2015-12-12 10:18
感谢楼主无私分享,真是及时雨啊!
回复 崇老先生 2016-3-22 15:26
谢谢楼主!!!!
回复 崇老先生 2016-3-22 15:27
zhangwenjing: 增量谐波平衡法,MATLAB编程
楼主可以加个口口探讨一哈吗?
回复 崇老先生 2016-3-22 15:29
zhangwenjing: 增量谐波平衡法,MATLAB编程
楼主可以加个口口探讨一哈吗?QQ:2079950490
微信号:cfq20131008
回复 ild_cyh 2016-10-2 19:56
zhangwenjing: 增量谐波平衡法,MATLAB编程
你好,请问算例中对k2的处理方法是建立在阻尼阵C是对角阵的基础上的吗?如果阻尼阵C不是对角阵的话,程序中的k2应该如何处理?谢谢。

——————————————————
k2=diag(S*A0).*diag(S1*A0)
fk2=inline(S'*k2*S);
K2=quadv(fk2,0,2*pi);
——————————————————

facelist doodle 涂鸦板

您需要登录后才可以评论 登录 | 我要加入

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

GMT+8, 2024-11-3 05:38 , Processed in 0.034051 second(s), 16 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

返回顶部