How to write a phase plane in matlab?
11 views (last 30 days)
Show older comments
I'm trying to write a phase plane for ODE.
The equation goes like: da/dt = -2*k1*a^2 + k2*b; db/dt = 2*k1*a^2 - k2*b, where k1 & k2 are const.
I want it to have arrows and lines and stuff like this:![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/287173/image.jpeg)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/287173/image.jpeg)
I'm still very new to matlab so it's quite painful to read long long unfamiliar codes. If you can put detailed annotations with codes I'd appreciate it very much!
Thank you!
3 Comments
Accepted Answer
Raunak Gupta
on 30 Apr 2020
Edited: Raunak Gupta
on 30 Apr 2020
Hi,
I have compiled some code which essentially plot the same figure you have. For plotting the straight lines, you must choose the starting points in either a or b. I am choosing those in variable y10 and y20 as seen from the graph. Finally, for getting the orange line you need to check where both da/dt and db/dt are zero. Since here both equations are the same, I have solved one and plotted the result at the end. You may put the title and other legend text as per required.
k1 = 3;
k2 = 1;
% Y is essentially [da/dt,db/dt]
f = @(t,Y) [-2*k1*(Y(1))^2 + k2*Y(2); 2*k1*(Y(1))^2 - k2*Y(2)];
% for plotting the vector field using quiver
y1 = linspace(0,1,21);
y2 = linspace(0,1,21);
[x,y] = meshgrid(y1,y2);
u = zeros(size(x));
v = zeros(size(x));
t=0;
for i = 1:numel(x)
% Calculating gradient value for a and b
Yprime = f(t,[x(i); y(i)]);
u(i) = Yprime(1);
v(i) = Yprime(2);
end
quiver(x,y,u,v,'r');
hold on
% Plotting the integrated value
for y10 = [0.3 0.5 0.7 0.8 1]
[ts,ys] = ode45(f,[0,50],[y10;0]);
plot(ys(:,1),ys(:,2),'k-')
end
for y20 = [0.35 0.65 0.8]
[ts,ys] = ode45(f,[0,50],[0;y20]);
plot(ys(:,1),ys(:,2),'k-')
end
grid on
% solve the equation for which gradient is zero and plotting it
syms a b
eqn = -2*k1*a^2 + k2*b == 0;
fimplicit(eqn, [0 1 0 1]);
hold off
Hope it is understandable.
More Answers (0)
See Also
Categories
Find more on Digital Filter Analysis 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!