How to write a phase plane in matlab?

Question

Moira 2020-4-24

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/520450-how-to-write-a-phase-plane-in-matlab

评论： Moira 2020-4-30

采纳的回答： Raunak Gupta

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 个评论
显示 1更早的评论隐藏 1更早的评论

BALAJI KARTHEEK 2020-4-24

use quiver for arrows

Moira 2020-4-24

The initial condition: a = 1, b = 0

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Raunak Gupta 2020-4-30

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/520450-how-to-write-a-phase-plane-in-matlab#answer_429329

编辑：Raunak Gupta 2020-4-30

在 MATLAB Online 中打开

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