How do you create phase planes in MATLAB?

I am new to MATLAB and creating phase planes seems to be harder than the simple commands in the tutorial. For a homework problem, I am to plot (where R is tho and B is beta)

d2x/dR2 = -x - x^3 - B dx/dR

for 2 cases (B=.1 and B=10)

What does the code look like for a phase plane in MATLAB? Thank you!

Answer 1

[Star Strider]

See if this does what you want:

ODE = @(t,x,B) [x(2); -x(1) - x(1).^3 - B*x(2)];

X0 = [1; 1]*0.1;
tspan = linspace(0, 50, 250);

B = 0.1;
[T,X1] = ode45(@(t,x) ODE(t,x,B), tspan, X0);

B = 1;
[T,X2] = ode45(@(t,x) ODE(t,x,B), tspan, X0);

figure(1)
plot(T,X1)
hold on
plot(T,X2)
hold off
grid
title('Time Series')
xlabel('\itt\rm')
legend('\itx(t)\rm', '\itdx(t)/dt\rm','\itx(t)\rm', '\itdx(t)/dt\rm')

figure(2)
plot(X1(:,1), X1(:,2))
hold on
plot(X2(:,1), X2(:,2))
hold off
grid
axis equal
title('Phase Plane')
xlabel('\itx(t)\rm')
ylabel('\itdx(t)/dt\rm')
legend('B = 0.1', 'B =    1', 'Location','SW')

[Amanda]

kind of. I get both phase planes but they are not arrows. As well that my b=.01 and b=10 should be two separate phase planes. More specifically: "Use Matlab to plot the phase plane of (2) for β = .1 and for β = 10. (So you are making two phase planes here, one for each value of β.)"

Thank you for the work though! I got the following from changing your code to be 2 Beta values but it's not working correctly, what do I do?::

ODE = @(t,x,B) [x(2); -x(1) - x(1).^3 - B*x(2)];
X0 = [1; 1]*0.1;
tspan = linspace(0, 50, 250);
B = 0.1;
[T,X1] = ode45(@(t,x) ODE(t,x,B), tspan, X0);
B = 1;
[T,X2] = ode45(@(t,x) ODE(t,x,B), tspan, X0);
figure(1)
plot(T,X1)
hold on
grid
title('Beta = 0.1')
xlabel('\itt\rm')
legend('\itx(t)\rm', '\itdx(t)/dt\rm','\itx(t)\rm', '\itdx(t)/dt\rm')
figure(2)
plot(T,X2)
hold on
grid
axis equal
title('Beta = 10')
xlabel('\itx(t)\rm')
ylabel('\itdx(t)/dt\rm')

[Star Strider]

I didn’t realise you wanted arrows. You also changed my code so that it doesn’t plot the phase planes!

This version plots them. Note that the phase planes are in figure(2), not figure(1). I plotted figure(1) to be sure the integration looked correct, and kept it in so you could verify that as well. You can eliminate figure(1) if you wish, since it is not otherwise important to the code.

The correct code:

ODE = @(t,x,B) [x(2); -x(1) - x(1).^3 - B*x(2)];
B = 0.1;
X0 = [1; 1]*0.1;
tspan = linspace(0, 50, 250);
[T,X1] = ode45(@(t,x) ODE(t,x,B), tspan, X0);
for k1 = 1:length(T)
    DX1(k1,:) = ODE(T(k1),X1(k1,:),B)';
end

B = 10;
[T,X2] = ode45(@(t,x) ODE(t,x,B), tspan, X0);
for k1 = 1:length(T)
    DX2(k1,:) = ODE(T(k1),X2(k1,:),B)';
end

figure(1)
plot(T,X1)
hold on
plot(T,X2)
hold off
grid
title('Time Series')
xlabel('\itt\rm')
legend('\itx(t)\rm', '\itdx(t)/dt\rm','\itx(t)\rm', '\itdx(t)/dt\rm')

figure(2)
quiver(X1(:,1), X1(:,2), DX1(:,1), DX1(:,2))
hold on
q2 = quiver(X2(:,1), X2(:,2), DX2(:,1), DX2(:,2));
hold off
grid
set(q2, 'AutoScaleFactor',10)
axis equal
title('Phase Plane')
xlabel('\itx(t)\rm')
ylabel('\itdx(t)/dt\rm')
legend('B = 0.1', 'B =  10', 'Location','SW')

The quiver plot requires (x,y) pairs for the origins of the arrows, and (dx,dy) pairs for their directions. I added the for loops to calculate those from the original ‘ODE’ function.

I originally had the second Beta value set to 1 because 10 isn’t very interesting (and I can never rule out typos in posts). I set it back to 10 in this code. I also used the set call to change the arrow scaling factor ('AutoScaleFactor') to 10 (empirically) in the plot for the second series (B=10) so you can see them. You can change that if you want to.

Also, you didn’t specify a time span (‘tspan’) or initial conditions (‘X0’) in your Question, so I chose mine arbitrarily. Change those as well if you want to.