声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 3298|回复: 3

[近似分析] 增量谐波平衡法 IHB FDB

[复制链接]
发表于 2016-5-12 15:36 | 显示全部楼层 |阅读模式

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

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

x
  大谢本程序的编写者和分享者!!!!!!
  与程序相关的有用的参考文献:《应用于范德波方程的增量谐波平衡法》唐南中山大学研究生专刊自然科学版

  1.   clear all

  2.   clc

  3.   tic

  4.   %%%定义各参数

  5.   syms t

  6.   w0=3;

  7.   epsR=0.001;

  8.   m=[1 0;0 1];

  9.   epsilon=0.24;

  10.   r=0.4;

  11.   delta=0.56;

  12.   k=[1+epsilon*r -epsilon*r;-epsilon*r 1+epsilon*delta];

  13.   Cs=[cos(t) cos(3*t) sin(t) sin(3*t)];

  14.   Cs1=diff(Cs,t,1);

  15.   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)];

  16.   A1=[1 1 1 1]';

  17.   A2=[1 1 1 1]';

  18.   A0=[A1;A2];

  19.   T1=[eye(4,4) zeros(4,4)];

  20.   T2=[zeros(4,4) eye(4,4)];

  21.   S2=diff(S,t,2);

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

  23.   M=quadv(fm,0,2*pi);

  24.   fk=inline(S'*k*S);

  25.   K=quadv(fk,0,2*pi);

  26.   S1=diff(S,t,1);

  27.   c=[1 0;0 1];

  28.   fc=inline(S'*c*S1);

  29.   C=quadv(fc,0,2*pi);

  30.   c3=diag(S*A0).^2;

  31.   fc3=inline(S'*c3*S1);

  32.   C3=quadv(fc3,0,2*pi);

  33.   k2=diag(S*A0).*diag(S1*A0)

  34.   fk2=inline(S'*k2*S);

  35.   K2=quadv(fk2,0,2*pi);

  36.   %%%%%带入推导出的公式

  37.   Kmc=w0^2*M+epsilon*w0*(C3-C)+K+2*epsilon*w0*K2;

  38.   R=-(w0^2*M+epsilon*w0*(C3-C)+K)*A0;

  39.   Rmc=-(2*w0*M+epsilon*(C3-C))*A0;

  40.   % AA首元素已知a1=0.0,求ww

  41.   a1=0.0;

  42.   %%%变换矩阵,使ww变量代替a1

  43.   Kmc11=-Rmc(:,1);

  44.   Kmcr=[Kmc11 Kmc(:,2:size(Kmc,2))];

  45.   %求未知变量

  46.   AA=inv(Kmcr)*R; % drtA1(1)

  47.   ww=AA(1); % drtW(1)

  48.   %%赋予新变量新值

  49.   A01=A0+[a1;AA(2:length(A0),1)]; % A(1)+drtA(1)

  50.   % Aw0=AA+A00; % A1(0)+drtA1(1)=A1(1)

  51.   w01=w0+ww; % W+drtW(1)

  52.   n=1;

  53.   tol=1;

  54.   %

  55.   while tol>epsR

  56.   A0=A01;

  57.   %

  58.   w0=w01

  59.   %

  60.   c3=diag(S*A0).^2;

  61.   fc3=inline(S'*c3*S1);

  62.   C3=quadv(fc3,0,2*pi);

  63.   k2=diag(S*A0).*diag(S1*A0)

  64.   fk2=inline(S'*k2*S);

  65.   K2=quadv(fk2,0,2*pi);

  66.   %%%%%带入推导出的公式

  67.   Kmc=w0^2*M+epsilon*w0*(C3-C)+K+2*epsilon*w0*K2;

  68.   R=-(w0^2*M+epsilon*w0*(C3-C)+K)*A0;

  69.   Rmc=-(2*w0*M+epsilon*(C3-C))*A0;

  70.   %%%%%

  71.   tol=norm(R)

  72.   if(n>1000)

  73.   disp('迭代步数太多,可能不收敛')

  74.   return;

  75.   end

  76.   Kmc11=-Rmc(:,1);

  77.   Kmcr=[Kmc11 Kmc(:,2:size(Kmc,2))];

  78.   %

  79.   AA=inv(Kmcr)*R;

  80.   ww=AA(1);

  81.   % A00=[w0;A0(2:6,1)];

  82.   A01=A0+[a1;AA(2:length(A0),1)];

  83.   w01=w0+ww;

  84.   n=n+1;

  85.   end

  86.   X0=S*A0;

  87.   dX0=S1*A0;

  88.   %%绘范德波图

  89.   tt=0:.1:10;

  90.   xo1=subs(X0(1),tt);

  91.   xo2=subs(X0(2),tt);

  92.   dxo1=subs(dX0(1),tt);

  93.   dxo2=subs(dX0(2),tt);

  94.   figure(1)

  95.   plot(xo1,dxo1,'b','linewidth',2)

  96.   hold on

  97.   plot(xo2,dxo2,'b','linewidth',2)

  98.   axis([-3 3 -3 3])

  99.   title('范德波极限环')

  100.   xlabel('x0')

  101.   ylabel('dx0')

  102.   toc
复制代码


30.jpg

31.jpg



转自:http://blog.sina.com.cn/s/blog_786ce14d0102v1pa.html

回复
分享到:

使用道具 举报

发表于 2017-5-20 21:32 | 显示全部楼层
请问一下,A1=[1 1 1 1]';   A2=[1 1 1 1]'; 是表示初值么? 可以改变么?
发表于 2020-2-2 13:53 | 显示全部楼层
lihaitao123 发表于 2017-5-20 21:32
请问一下,A1=[1 1 1 1]';   A2=[1 1 1 1]'; 是表示初值么? 可以改变么?

是初值,可以修改
发表于 2020-2-2 16:04 | 显示全部楼层
这个程序太早了,用的inline和quadv函数,这些函数再新版matlab里会提示:未来版本中会移除
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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

GMT+8, 2024-3-29 07:01 , Processed in 0.067257 second(s), 25 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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