How can I draw a bifurcation diagram for the given code?
2 views (last 30 days)
Show older comments
clc
close all
k1 = 3;
k_1 = 1;
k2 = 2.5;
k3 = 1;
k_3 = 1;
k4 = 2;
k5 = 1;
E1 = 1;
E2 = 2;
K1 = (k_1+k2)/k1;
K2 = (k_3+k4)/k3;
a = -(k5*E1)^2;
b = -((k3*E2+(k5*E1)^2)*(k1*(K1+E1)+k3*K2));
c= -(k3*E2*(k1*(K1+E1)+k4)+k1*(k5*E1)^2*(k2*E1+k3*K2*(K1+E1)));
d = -(k1*k3*(k4*K1*E2+E1*E2*(k2+k4)+k2*k5^2*K2*E1^3));
e=-(k5*E1);
f=-(k5*E1*(k1*K1+k3*(K2+E2)));
g = -(k1*k3*k5*E1*K1*(K2+E2));
m1 = -5*(a*e)^2;
m2 = (8*a^2*e*g-4*a^2*f^2-4*b^2*e^2+8*a*e^2*c);
m3 = (2*e*g-f^2)*(3*b^2-6*a*c)+3*e^2*(2*b*d-c^2)-3*a^2*g^2;
m4 = 2*(2*e*g-f^2)*(c^2-2*b*d)+2*g^2*(2*a*c-b^2)-2*e^2*d^2;
m5 = d^2*(2*e*g-f^2)+g^2*(2*b*d-c^2);
w = 0.001:0.01:0.9;
T = (m1*w^8+m2*w^6+m3*w^4+m4*w^2+m5)*((g-e*w^2)^2+(f*w)^2)*((d-b*w)^2+(c*w-a*w^3)^2);
0 Comments
Accepted Answer
Chunru
on 27 Aug 2022
When do array operation, you should use .^ or .*:
k1 = 3;
k_1 = 1;
k2 = 2.5;
k3 = 1;
k_3 = 1;
k4 = 2;
k5 = 1;
E1 = 1;
E2 = 2;
K1 = (k_1+k2)/k1;
K2 = (k_3+k4)/k3;
a = -(k5*E1)^2;
b = -((k3*E2+(k5*E1)^2)*(k1*(K1+E1)+k3*K2));
c= -(k3*E2*(k1*(K1+E1)+k4)+k1*(k5*E1)^2*(k2*E1+k3*K2*(K1+E1)));
d = -(k1*k3*(k4*K1*E2+E1*E2*(k2+k4)+k2*k5^2*K2*E1^3));
e=-(k5*E1);
f=-(k5*E1*(k1*K1+k3*(K2+E2)));
g = -(k1*k3*k5*E1*K1*(K2+E2));
m1 = -5*(a*e)^2;
m2 = (8*a^2*e*g-4*a^2*f^2-4*b^2*e^2+8*a*e^2*c);
m3 = (2*e*g-f^2)*(3*b^2-6*a*c)+3*e^2*(2*b*d-c^2)-3*a^2*g^2;
m4 = 2*(2*e*g-f^2)*(c^2-2*b*d)+2*g^2*(2*a*c-b^2)-2*e^2*d^2;
m5 = d^2*(2*e*g-f^2)+g^2*(2*b*d-c^2);
w = 0.001:0.01:0.9;
T = (m1*w.^8+m2*w.^6+m3*w.^4+m4*w.^2+m5).*((g-e*w.^2).^2+(f*w).^2).*((d-b*w).^2+(c*w-a*w.^3).^2);
plot(w, T)
More Answers (0)
See Also
Categories
Find more on Characters and Strings in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!