Effacer les filtres
Effacer les filtres

How to write a phase plane in matlab?

48 vues (au cours des 30 derniers jours)
Moira
Moira le 24 Avr 2020
Commenté : Moira le 30 Avr 2020
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:
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 commentaires
BALAJI KARTHEEK
BALAJI KARTHEEK le 24 Avr 2020
use quiver for arrows
Moira
Moira le 24 Avr 2020
The initial condition: a = 1, b = 0

Connectez-vous pour commenter.

Réponse acceptée

Raunak Gupta
Raunak Gupta le 30 Avr 2020
Modifié(e) : Raunak Gupta le 30 Avr 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.
  1 commentaire
Moira
Moira le 30 Avr 2020
It's so helpful!!! Thank you so much!!!

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

En savoir plus sur 2-D and 3-D Plots dans Help Center et File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by