How to write a phase plane in matlab?
11 次查看(过去 30 天)
显示 更早的评论
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 个评论
采纳的回答
Raunak Gupta
2020-4-30
编辑:Raunak Gupta
2020-4-30
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.
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Digital Filter Analysis 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!