% frequency response curve (Primary Resonance) of yabuno model controlled by ppf
% plotting a1 and a2 versus segma1 with different values of segma 2
clear all
% constants of the system
alpha=.894 ; k1=.5 ; k2=.547 ; k3=.447 ; f=.05 ; mu1=.02 ; mu2=.001 ; c1=.3 ; c2=.3 ; w=1 ; w2=1 ; s2=0.0;
i=1; j=1; k=1;
for s1=-.3:.00005:.3
Z=roots([36*alpha^2*w2^8*s1^8/c2^8+36*alpha^2*w2^8*s2^8/c2^8+36*alpha^2*w2^8*mu2^8/c2^8+2520*alpha^2*w2^8*s2^4*s1^4/c2^8-2016*alpha^2*w2^8*s2^3*s1^5/c2^8-288*alpha^2*w2^8*s2^7*s1/c2^8+1008*alpha^2*w2^8*s2^6*s1^2/c2^8-2016*alpha^2*w2^8*s2^5*s1^3/c2^8+1008*alpha^2*w2^8*s2^2*s1^6/c2^8-288*alpha^2*w2^8*s2*s1^7/c2^8-864*alpha^2*w2^8*mu2^4*s2^3*s1/c2^8+1296*alpha^2*w2^8*mu2^4*s2^2*s1^2/c2^8-864*alpha^2*w2^8*mu2^4*s2*s1^3/c2^8+216*alpha^2*w2^8*mu2^4*s1^4/c2^8+144*alpha^2*w2^8*mu2^6*s2^2/c2^8-288*alpha^2*w2^8*mu2^6*s2*s1/c2^8+144*alpha^2*w2^8*mu2^6*s1^2/c2^8+216*alpha^2*w2^8*mu2^4*s2^4/c2^8+144*alpha^2*w2^8*mu2^2*s2^6/c2^8-864*alpha^2*w2^8*mu2^2*s2^5*s1/c2^8+2160*alpha^2*w2^8*mu2^2*s2^4*s1^2/c2^8-2880*alpha^2*w2^8*mu2^2*s2^3*s1^3/c2^8+2160*alpha^2*w2^8*mu2^2*s2^2*s1^4/c2^8-864*alpha^2*w2^8*mu2^2*s2*s1^5/c2^8+144*alpha^2*w2^8*mu2^2*s1^6/c2^8 ...
-48*s1^7*w2^6*alpha/c2^6-144*s1*w2^6*mu2^4*alpha*s2^2/c2^6+288*s1^2*w2^6*mu2^4*alpha*s2/c2^6-864*s1^3*w2^6*mu2^2*alpha*s2^2/c2^6+576*s1^4*w2^6*mu2^2*alpha*s2/c2^6-144*s1^5*w2^6*mu2^2*alpha/c2^6-48*s1*w2^6*mu2^6*alpha/c2^6-144*s1^3*w2^6*mu2^4*alpha/c2^6-144*s1*w2^6*mu2^2*alpha*s2^4/c2^6+576*s1^2*w2^6*mu2^2*alpha*s2^3/c2^6-12*alpha*w2^5*s2^5*c1/c2^5+60*alpha*w2^5*s2^4*c1*s1/c2^5-120*alpha*w2^5*s2^3*s1^2*c1/c2^5-48*s1*w2^6*s2^6*alpha/c2^6+288*s1^2*w2^6*s2^5*alpha/c2^6-720*s1^3*w2^6*s2^4*alpha/c2^6+960*s1^4*w2^6*s2^3*alpha/c2^6-720*s1^5*w2^6*s2^2*alpha/c2^6+120*alpha*w2^5*s2^2*s1^3*c1/c2^5-60*alpha*w2^5*s2*s1^4*c1/c2^5+12*alpha*w2^5*s1^5*c1/c2^5-12*alpha*w2^5*mu2^4*c1*s2/c2^5+288*s1^6*w2^6*s2*alpha/c2^6+12*alpha*w2^5*mu2^4*c1*s1/c2^5-24*alpha*w2^5*mu2^2*s2^3*c1/c2^5+72*alpha*w2^5*mu2^2*s2^2*c1*s1/c2^5-72*alpha*w2^5*mu2^2*s2*s1^2*c1/c2^5+24*alpha*w2^5*mu2^2*s1^3*c1/c2^5 ...
16*mu1^2*w2^4*mu2^4/c2^4+c1^2*mu2^2*w2^2/c2^2+16*mu1^2*w2^4*s1^4/c2^4+16*mu1^2*w2^4*s2^4/c2^4+16*s1^2*w2^4*mu2^4/c2^4+16*s1^2*w2^4*s2^4/c2^4+c1^2*w2^2*s1^2/c2^2+c1^2*w2^2*s2^2/c2^2+32*s1^4*w2^4*mu2^2/c2^4-64*s1^3*w2^4*s2^3/c2^4-64*s1^5*w2^4*s2/c2^4+96*s1^4*w2^4*s2^2/c2^4-8*s1^4*w2^3*c1/c2^3+16*s1^6*w2^4/c2^4+32*mu1^2*w2^4*mu2^2*s1^2/c2^4+8*mu1*w2^3*mu2^3*c1/c2^3+96*mu1^2*w2^4*s2^2*s1^2/c2^4-64*mu1^2*w2^4*s2^3*s1/c2^4-64*mu1^2*w2^4*s2*s1^3/c2^4+8*mu1*w2^3*s1^2*c1*mu2/c2^3+32*s1^2*w2^4*mu2^2*s2^2/c2^4-64*s1^3*w2^4*mu2^2*s2/c2^4+8*s1*w2^3*mu2^2*c1*s2/c2^3-8*s1^2*w2^3*mu2^2*c1/c2^3+8*s1*w2^3*s2^3*c1/c2^3-24*s1^2*w2^3*s2^2*c1/c2^3-2*c1^2*w2^2*s2*s1/c2^2+24*s1^3*w2^3*s2*c1/c2^3-64*mu1^2*w2^4*mu2^2*s2*s1/c2^4+8*mu1*w2^3*s2^2*c1*mu2/c2^3+32*mu1^2*w2^4*mu2^2*s2^2/c2^4-16*mu1*w2^3*s2*s1*c1*mu2/c2^3 ...
-k1^2*f^2*w2^2*s1^2/c2^2+2*k1^2*f^2*w2^2*s2*s1/c2^2-k1^2*f^2*w2^2*mu2^2/c2^2-k1^2*f^2*w2^2*s2^2/c2^2]);
if imag(Z(1))==0 && real(Z(1))>0
s_1(i)=s1;
z_1(i)=Z(1)^.5;
y_1(i)=2*Z(1)^.5/c2*sqrt(mu2^2+(s2-s1)^2);
i=i+1;
end
if imag(Z(2))==0 && real(Z(2))>0
s_2(j)=s1;
z_2(j)=Z(2)^.5;
y_2(j)=2*Z(2)^.5/c2*sqrt(mu2^2+(s2-s1)^2);
j=j+1;
end
if imag(Z(3))==0 && real(Z(3))>0
s_3(k)=s1;
z_3(k)=Z(3)^.5;
y_3(k)=2*Z(3)^.5/c2*sqrt(mu2^2+(s2-s1)^2);
k=k+1;
end
end
% combining vectors
s1=[ s_1 s_2 s_3 ];
y1=[ y_1 y_2 y_3 ];
y2=[ z_1 z_2 z_3 ];
% checking stability
j=1; k=1;
for i=1:length(s1)
% defining sines and cosines
sp1=2/k1/f*(mu1*y1(i)+c1*mu2/c2*y2(i)^2/y1(i));
cp1=2/k1/f*(-s1(i)*y1(i)+3/8*alpha*y1(i)^3-c1/c2*(s2-s1(i))*y2(i)^2/y1(i));
sp2=-2*mu2/c2*y2(i)/y1(i);
cp2=2/c2*(s2-s1(i))*y2(i)/y1(i);
% defining Jacobian matrix
J=[-mu1 k1/2*f*cp1 c1/2*sp2 c1/2*cp2*y2(i);-3/4*alpha*y1(i)-c1/2*y2(i)/y1(i)^2*cp2-k1/2/y1(i)^2*f*cp1 -k1/2/y1(i)*f*sp1...
c1/2/y1(i)*cp2 -c1/2*y2(i)/y1(i)*sp2 ; -c2/2*sp2 0 -mu2 -c2/2*y1(i)*cp2 ; -3/4*alpha*y1(i)-c1/2*y2(i)/y1(i)^2*cp2-k1*f/2/y1(i)^2*cp1-c2/2/y2(i)*cp2...
-k1*f/2/y1(i)*sp1 c1/2/y1(i)*cp2+c2/2*y1(i)/y2(i)^2*cp2 -c1/2*y2(i)/y1(i)*sp2+c2/2*y1(i)/y2(i)*sp2];
% defining stability criterion
if real(eig(J))<0
s1s(j)=s1(i);
y1s(j)=y1(i);
y2s(j)=y2(i);
j=j+1;
else
s1us(k)=s1(i);
y1us(k)=y1(i);
y2us(k)=y2(i);
k=k+1;
end
end
% plotting a1 with segma1
plot(s1s,y1s,'k.',s1us,y1us,'k--')
axis([-.3 .3 0 1])
% plotting a2 with segma1
figure
plot(s1s,y2s,'k.',s1us,y2us,'k--')
axis([-.3 .3 0 1])
