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

 Réponse acceptée

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 commentaires

@VBBV, there's a typo in the comment and I think the x and y label commands should be included in the loop.
yes, xlabel and ylabel commands need to be present inisde the loop
Thanks for the assistance. the result is not producing a phase space plot.

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

En savoir plus sur Graphics Performance dans Centre d'aide et File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by