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

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

日志

duffing方程的时间历程,相图,频谱和庞加莱映射

已有 687 次阅读2009-9-13 17:51 |

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'); 

 

全部作者的其他最新日志

评论 (0 个评论)

facelist doodle 涂鸦板

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

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

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.

返回顶部