增量谐波平衡法 IHB FDB
大谢本程序的编写者和分享者!!!!!!与程序相关的有用的参考文献：《应用于范德波方程的增量谐波平衡法》唐南中山大学研究生专刊自然科学版
clear all
clc
tic
%%%定义各参数
syms t
w0=3;
epsR=0.001;
m=;
epsilon=0.24;
r=0.4;
delta=0.56;
k=;
Cs=;
Cs1=diff(Cs,t,1);
S=;
A1=';
A2=';
A0=;
T1=;
T2=;
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=;
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=;
%求未知变量
AA=inv(Kmcr)*R; % drtA1(1)
ww=AA(1); % drtW(1)
%%赋予新变量新值
A01=A0+; % 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=;
%
AA=inv(Kmcr)*R;
ww=AA(1);
% A00=;
A01=A0+;
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
转自：http://blog.sina.com.cn/s/blog_786ce14d0102v1pa.html
请问一下，A1='; A2='; 是表示初值么？ 可以改变么？ lihaitao123 发表于 2017-5-20 21:32
请问一下，A1='; A2='; 是表示初值么？ 可以改变么？
是初值，可以修改 这个程序太早了，用的inline和quadv函数，这些函数再新版matlab里会提示：未来版本中会移除
页:
[1]