TA的每日心情 | 郁闷 15-11-2017 17:51 |
---|
签到天数: 4 天 [LV.2]偶尔看看I
|
发表于 30-10-2013 20:41:40
|
显示全部楼层
大林木 发表于 30-10-2013 20:40
我遇到个NATLAB问题 想在线请教您一下
主函数
global DF LEN_SEGM FFT_SEGMENT flag
flag=0;
DF=1000/1024;
LEN_SEGM=0;
x0=[0,0,2,3]';
[t,x1]=odeRK4(@shake,[0:0.001:10],x0);
FFT_SEGMENT=zeros(1,1024);
flag=1;
DF=1000/1024;
LEN_SEGM=0;
x0=[0,0,2,3]';
[t,x2]=odeRK4(@shake,[0:0.001:10],x0);
n=length(t);
figure(1)
plot(t,x1(:,1),t,x2(:,1))
legend('原系统','吸振系统')
被调用的函数1
function dx=shake(t,x)
global FFT_SEGMENT LEN_SEGM flag
M1=500;
M2=45;
K2=22e3;
K1=192e3;
c1=5e3;
c2=0.5e3;
M=[M1 0;0 M2];
K=[K1+K2 -K2;-K2 K2];
C=[c1+c2 -c2;-c2 c2];
F=[50*sin((50)*pi*t);0];
n=2;
if ((flag==1) & (LEN_SEGM==1024))
K2=ModifyK2();
end
K=[K1+K2 -K2;-K2 K2];
dx(1:n,1)=x(n+1:2*n,1);
dx(n+1:2*n,1)=M\(-C*x(n+1:2*n,1)-K*x(1:n,1)+F);
被调用函数2
function [T Y]=odeRK4(f,tspan,y0)
global FFT_SEGMENT LEN_SEGM flag
nRow=length(tspan);
nCol=length(y0);
Y=zeros(nRow,nCol);
h=tspan(2)-tspan(1);
for i=1:nRow-1
Y(i,:)=y0;
t=tspan(i);
K1=f(t,y0);
K2=f(t+h/2,y0+K1*h/2);
K3=f(t+h/2,y0+K2*h/2);
K4=f(t+h,y0+h*K3);
y0=y0+(K1+2*K2+2*K3+K4)*h/6;
if flag==1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if LEN_SEGM<1024
LEN_SEGM=LEN_SEGM+1;
FFT_SEGMENT(LEN_SEGM)=y0(1);
else
FFT_SEGMENT=[FFT_SEGMENT(2:1024) y0(1)];
% K2=ModifyK2();
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
end
Y(nRow,:)=y0;
T=tspan;
被调用函数3
function K2=ModifyK2()
global FFT_SEGMENT DF
BLOKE=zeros(1,100);
M2=45;
X=fft(FFT_SEGMENT);
X=abs(X(1:512));
[Xmax idx]=max(X);
f=(idx-1)*DF*2*pi;
K2=f.^2*M; |
|