|
function dy=duffing(t,x)
omega=1;%定义参数
f1=x(2);
f2=-omega^2*x(1)-x(1)^3;
dy=[f1;f2];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
step=0.01
y0=[0.1;0];
tspan=[0:step:500];
[t,y]=ode45('duffing',tspan,y0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(2,2,1)
plot(t,y(:,1));%绘制y的时间历程
xlabel('t')%横轴为t
ylabel('y')%纵轴为y
grid;%显示网格线
axis([460 500 -Inf Inf])%显示范围设置
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(2,2,2)
plot(y(end-1000:end,1),y(end-1000:end,2));%绘制y的时间历程
xlabel('y')%横轴为y
ylabel('dy/dt')%纵轴为dy/dt
grid;%显示网格线
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(2,2,3)
yy=fft(y(end-1000:end,1));
N=length(yy);
power=abs(yy);
freq=(1:N-1)*1/step/N;
plot(freq(1:N/2),power(1:N/2));
xlabel('f(y)')
ylabel('y')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(2,2,4)
count=find(t>100);%截取稳态过程
y=y(count,:);
n=length(y(:,1));%计算点的总数
for i=2:n-1
if y(i-1,1)+eps<y(i,1) && y(i,1)>y(i+1,1)+eps % 简单的取出局部最大值
plot(y(i,1),y(i,2),'.');hold on;
end
end
xlabel('y');ylabel('dy/dt');
GMT+8, 2024-11-2 19:27 , Processed in 0.043398 second(s), 15 queries , Gzip On.
Powered by Discuz! X3.4
Copyright © 2001-2021, Tencent Cloud.