I got a blank figure while trying to plot the phase space

3 views (last 30 days)
I got a blank figure while trying to plot the phase space of this code.
This is the function file
function dxdt = IP(t,x,B,C,D,r,z,A,K)
dxdt = zeros(4,1);
dxdt(1) = x(3);
dxdt(2) = x(4);
dxdt(3) = (z*D*exp(-z*x(1))*(2*A*exp(2*z*r)*exp(-z*x(1)))+(B*exp(z*r))*(1+exp(-(z*x(1))))+2*C)/((1-exp(-(z*x(1))))^3)+K*(x(2)-x(1));
dxdt(4) = (z*D*exp(-z*x(2))*(2*A*exp(2*z*r)*exp(-z*x(2)))+(B*exp(z*r))*(1+exp(-(z*x(2))))+2*C)/((1-exp(-(z*x(2))))^3)+K*(x(1)-x(2));
end
this is the script:
close all
D=0.02;
r=0.04;
A=2.5;
B=0.001;
z=3.0;
C=3.25;
q1=0;
q2=0;
p1=0.632;
E=0.2;
p2=sqrt(2*E- p1^2);
dt=1/5;
tspan=[10:dt:100];
K=[0 0.050 0.5 1.5];
m=length(K);
for i=1:m;
[t x]=ode45(@(t,x)IP(t,x,B,C,D,r,z,A,K(i)),[0 100],[q1 q2 p1 p2]);
figure(1)
subplot(2,2,(i))
plot(x(:,1),x(:,3),'-b',x(:,2),x(:,4),'-r ')
xlabel('q');
ylabel('p');
end

Accepted Answer

VBBV
VBBV on 18 May 2023
Choose better initial conditions e.g. for q1 and q2 they are both zero
close all
D=0.02;
r=0.04;
A=2.5;
B=0.01;
z=3.0;
C=3.25;
%>>>>>>>>> choose hetter intial coditions
q1=1;
q2=0.5;
%-->>>>>>>>>
p1=0.632;
E=0.2;
p2=sqrt(2*E- p1^2);
dt=1/5;
tspan=[10:dt:100];
K=[0 0.050 0.5 1.5];
m=length(K);
figure(1)
hold on
for i=1:m;
[t,x]=ode45(@(t,x)IP(t,x,B,C,D,r,z,A,K(i)),[0 100],[q1 q2 p1 p2]);
subplot(2,2,(i))
plot(x(:,1),x(:,3),'-b',x(:,2),x(:,4),'-r ')
end
xlabel('q');
ylabel('p');
function dxdt = IP(t,x,B,C,D,r,z,A,K)
dxdt = zeros(4,1);
dxdt(1) = x(3);
dxdt(2) = x(4);
dxdt(3) = (z*D*exp(-z*x(1))*(2*A*exp(2*z*r)*exp(-z*x(1)))+(B*exp(z*r))*(1+exp(-(z*x(1))))+2*C)/((1-exp(-(z*x(1))))^3)+K*(x(2)-x(1));
dxdt(4) = (z*D*exp(-z*x(2))*(2*A*exp(2*z*r)*exp(-z*x(2)))+(B*exp(z*r))*(1+exp(-(z*x(2))))+2*C)/((1-exp(-(z*x(2))))^3)+K*(x(1)-x(2));
end
  3 Comments
VBBV
VBBV on 18 May 2023
yes, xlabel and ylabel commands need to be present inisde the loop
OGUNGBEMI Ezekiel
OGUNGBEMI Ezekiel on 18 May 2023
Thanks for the assistance. the result is not producing a phase space plot.

Sign in to comment.

More Answers (0)

Categories

Find more on Mathematics 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!