求助 matlab求解duffing方程

微风breeze

刚接触非线性振动不久，就试着用matlab求解一个Duffing方程，结果在幅频曲线上没有出现跳跃现象，是怎么回事？
matlab程序：
Duffing.m
function dy = Duffing( t,y )
dy=zeros(2,1);
dy(1)=y(2);
%f=0.35*cos(2*pi*0.1*t);
f=0;
dy(2)=f-(0.2*y(2)+y(1)+0.5*y(1)^3);
end

Duffing_simulation.m
clear;
clc;
[T,Y]=ode45('Duffing',[0,50],[0 1]);
H1=figure;
%相图
plot(Y(:,1),Y(:,2));
grid;
figure;
%动态响应图
plot(T,Y(:,1));
grid;
figure;
%快速傅立叶变换
fs=30;
N=1024;
YY=fft(Y(:,1),N);
mag=abs(YY)*2/N;
mag(1)=mag(1)/2;
f=(0:(N/2))'*fs/N;
mag=mag(1:N/2+1);
%幅频特性曲线
stem(f,mag)
grid;
%相频特性曲线
for i=1:(N/2+1)
    Theta(i)=angle(YY(i));
    Theta(i)=Theta(i)*180/pi;
end
figure;
plot(f,Theta)
grid;
是程序哪块出现错误了？还是什么别的问题？求助，谢谢

无水1324

幅频曲线的跳跃和你程序中的FFT计算的结果可能不是一回事

微风breeze

回复 2 # 无水1324 的帖子

那怎样才能得到具有跳跃现象的幅频图了？如果不用FFT的话？还有其他的方法能将时域信息转换到频域上了？

微风breeze

回复 2 # 无水1324 的帖子

想了一下，知道了，用FFT变换是没法得到具有跳跃现象的幅频图的。。。

振动联盟浪子

回复 4 # 微风breeze 的帖子

我也想用matlab解非线性方程，matlab用的不多，还不怎么会用，我先简单看了一下，你的第一个文件中，f=0;是干什么用的？

微风breeze

回复 5 # 振动联盟浪子 的帖子

就是自由振动，没有外界激励

振动联盟浪子

回复 6 # 微风breeze 的帖子

你这个程序，如果考虑不是自由振动，有激励，并且如果非线性部分比较复杂，也能解吗？ 这个解法是直接用软件里面的ode45来实现，它和用一般的数值解法如龙格库塔法区别大吗？

无水1324

ode45就是龙格库塔法的。这个程序通用性比较强，非线性比较复杂，你还可以试一下，ode15，ode23，ode87，等等，应该可以求解

yzsldj

不需要FFT，Matlab程序：

close
ezplot('4*(.0025*x^2+(x-1-.3*y^2/8)^2)*y^2-1')  
axis([0,3,0,6]);
title('4 \omega_0 [\alpha^2 \omega^2+(\omega- \omega_0-3\epsilon\omega_0A^2/8)^2]A^2=F^2');
xlabel('频率 \omega');    %加X轴说明
ylabel('振幅A');                  

yzsldj

向左拐头，Matlab程序：

close
ezplot('4*(.0025*x^2+(x-1+.6*y^2/8)^2)*y^2-1')  
axis([-3,3,0,9]);
title('4 \omega_0 [\alpha^2 \omega^2+(\omega- \omega_0-3\epsilon\omega_0A^2/8)^2]A^2=F^2');
xlabel('频率 \omega');    %加X轴说明
ylabel('振幅A');                  

微风breeze

回复 10 # yzsldj 的帖子

请问一下，你这个ezplot里面的式子表示的什么？还有这个式子是怎么得到的？

微风breeze

回复 11 # 微风breeze 的帖子

这个式子是已有的解出的幅频关系式？如果模型不止一个自由度，没法解出解析解（或者近似解），那么又该如何得到幅频特性曲线了？

yzsldj

这个式子是已经解出的幅频关系式，多自由度没有考虑过

zaitala

这个需要仔细 了解一下才行。

zzqmessi

冒昧的问一下 对于跳跃现象 用什么方法可以得到