Drawing a bode plot of the Duffing equation, which is a special application of the integro-differential equation model.

I want to see the frequency responses of the Y1 signal in the modeling of this nonlinear duffing equation I prepared. I applied fft with y(:,1) but I'm not sure what to do to draw the bode diagram. Here I am having difficulty in creating the frequency vector corresponding to the response data.

% Parameters
c1 = 0.3; 
k1 = 3; 
k2 = 0.5; 
w = 1.7; % rad/s
A = 0.5;

% Initial conditions
y0 = [0; 0]; % y1(0) = 0, y2(0) = 0

%Time span
tspan = [0, 50]; 

% ODE Solver
[t, y] = ode45(@duffing, tspan, y0, [], c1, k1, k2, w, A); 

Y1 = fft(y(:,1)); % y1 signal Fourier transform

% Duffing equation function
function dydt = duffing(t,y,c1,k1,k2,w,A)
dydt = zeros(2,1); 
dydt(1) = y(2); 
dydt(2) = -c1*y(2) - k1*y(1) - k2*y(1)^3 + A*sin(w.*t); 
end

Answer 1

[Star Strider]

The problem with using a two-element vector for ‘tspan’ is that the MATLAB ODE integrators do not produce regularly-spaced time values with it, and the fft function requires that they must be regularly-spaced (constant sampling intervals and constant sampling frequency). To get them to be regularly-spaced, it is necessary to define them as such.

Try this —

% Parameters
c1 = 0.3; 
k1 = 3; 
k2 = 0.5; 
w = 1.7; % rad/s
A = 0.5;

% Initial conditions
y0 = [0; 0]; % y1(0) = 0, y2(0) = 0

%Time span
% tspan = [0, 50]; 
Fs = 1000;                                          % Sampling Frequency (Hz)
L = 60;                                             % Signal Length (s)
tspan = linspace(0, L*Fs, L*Fs+1)/Fs                % Time Vector

tspan = 1×60001

0 0.0010 0.0020 0.0030 0.0040 0.0050 0.0060 0.0070 0.0080 0.0090 0.0100 0.0110 0.0120 0.0130 0.0140 0.0150 0.0160 0.0170 0.0180 0.0190 0.0200 0.0210 0.0220 0.0230 0.0240 0.0250 0.0260 0.0270 0.0280 0.0290

% ODE Solver
[t, y] = ode45(@duffing, tspan, y0, [], c1, k1, k2, w, A); 

% Y1 = fft(y(:,1)); % y1 signal Fourier transform

[Y1,Fv] = FFT1(y(:,1),t);

[Y1max,idx] = max(mag2db(abs(Y1)));
Maximum_Magnitude_dB = Y1max

Maximum_Magnitude_dB = -8.1218

Peak_Frequency_Hz = Fv(idx)

Peak_Frequency_Hz = 0.2747

figure
subplot(2,1,1)
semilogx(Fv, mag2db(abs(Y1)))
ylabel('Magnitude (dB)')
grid
axis('padded')
subplot(2,1,2)
semilogx(Fv, rad2deg(unwrap(angle(Y1))))
grid
xlabel('Frequency (Hz)')
ylabel('Phase  (°)')
axis('padded')

% Duffing equation function
function dydt = duffing(t,y,c1,k1,k2,w,A)
dydt = zeros(2,1); 
dydt(1) = y(2); 
dydt(2) = -c1*y(2) - k1*y(1) - k2*y(1)^3 + A*sin(w.*t); 
end

function [FTs1,Fv] = FFT1(s,t)
    s = s(:);
    t = t(:);
    L = numel(t);
    Fs = 1/mean(diff(t));
    Fn = Fs/2;
    NFFT = 2^nextpow2(L);
    FTs = fft((s - mean(s)).*hann(L), NFFT)/sum(hann(L));
    Fv = linspace(0, 1, NFFT/2+1)*Fn;
    Iv = 1:numel(Fv);
    FTs1 = FTs(Iv);
end

.