How to fix an issue with plotting the Duffing Frequency Response Plot?

Dear all,

I am trying to plot the Frequency Response Plot of the Duffing Oscillator (I have scripted my own code but it did not go well). I got a code from sb else which looks more sophisticated but this one also runs into issues.

It does not plot the function entirely,

Any hints?

I attach the code down here as well as the correct plot it is supposed to look like:

Thanks in davance!

But I guess the incomplete plot like the follwoing:

clear all;
% Given values
eps = 1/30; % Nonlinear Parameter
alpha =0.937; % Forcing amplitude
T = 1000; % linear damping
ib=2.8; % Primary resonance
m=0.02;
w0=30;
a=linspace(0.01,8.5,100); % Range of a
%Finding the values of excitation frequency Omega/w0= F and G.
% And corresponding eigen values
%
for ii=1:1:length(a)
F(ii)=1-(eps^2)/2*((a(ii)^2)/12-(T*ib*m/(2*a(ii)))^2)+eps*T*ib*m/(2*a(ii))*sqrt((eps*T*ib*m/(2*a(ii)))^2+4-(eps*a(ii))^2/6);
lamF1=w0*sqrt(-(F(ii)-1+(1+2*eps)/24*(eps*a(ii))^2)*(F(ii)-1+((eps*a(ii))^2)/24))-(eps*alpha);
lamF2=-w0*sqrt(-(F(ii)-1+(1+2*eps)/24*(eps*a(ii))^2)*(F(ii)-1+((eps*a(ii))^2)/24))-(eps*alpha);
G(ii)=1-eps^2/2*((a(ii)^2)/12-(T*ib*m/(2*a(ii)))^2)-eps*T*ib*m/(2*a(ii))*sqrt((eps*T*ib*m/(2*a(ii)))^2+4-((eps*a(ii))^2)/6);
lamG1=w0*sqrt(-(G(ii)-1+(1+2*eps)*((eps*a(ii))^2)/24)*(G(ii)-1+((eps*a(ii))^2)/24))-(eps*alpha);
lamG2=-w0*sqrt(-(G(ii)-1+(1+2*eps)*((eps*a(ii))^2)/24)*(G(ii)-1+((eps*a(ii))^2)/24))-(eps*alpha);
if lamF1==conj(lamG1)% % To terminate the loop
break
else
%if gamma>0 % For spring hardening effect
plot(G(ii),a(ii),'b.');
hold on
FG=(F(ii)-1+(1+2*eps)/24*(eps*a(ii))^2)*(F(ii)-1+((eps*a(ii))^2)/24)*w0^2+(eps*a(ii))^2;
if FG<0
plot(F(ii),a(ii),'r.');
else
plot(F(ii),a(ii),'b.');
end
%else
plot(F(ii),a(ii),'b.');
hold on
GF=(G(ii)-1+(1+2*eps)/24*(eps*a(ii))^2)*(G(ii)-1+((eps*a(ii))^2)/24)*w0^2+(eps*a(ii))^2;
if GF<0
plot(G(ii),a(ii),'r.');
   else
plot(G(ii),a(ii),'b.');
end
 plot(F(ii),a(ii),'b.');
%end
end
end

xlim([-2.5,2.5]); % Range of x-axis
ylim([0,70]); %Range of y-axis
xlabel('\Omega/\omega_{0}'); % Label of x-axis
ylabel('a'); % Label of y-axis
title('Frequency response of a nonlinear system with \beta=0.05, q=0.2\gamma=0.5, \epsilon=1');

Answer 1

[Akshat Dalal]

Hello Alireza,

I understand that you have want to plot the Duffing Frequency Response Plot. You could refer to the following paper by Dr. Ashok Kumar Pandey from IIT Hyderabad to get a better understanding - https://people.iith.ac.in/ashok/NLO/Duffing_Oscillator.pdf