Plot dynamics of an ODE

8 次查看(过去 30 天)
I have an ODE and I need to plot with respect to and varying p. Following is the sample plot required.
I managed to draw one plot for and following is the sample code.
[xi,x] = meshgrid(0:0.001:1,0:0.002:1);
p =1
f2 = @(x,a,b) a*b*(x)./(b-2+x));
f1 = @(x,a,b) a*b*(1-x)./(2*b-x);
rs = @(x1,x2,a,b,p) x1.*f1(x,a,b) - p.*(1-x1).*f2(x,a,b);
figure
contour(x1,2,rs(x1,x2,1,2,1),[0 0]);
Next I added a loop for p as follows:
np = 10;
phi = linspace(0,1,np);
for p = 1:np
f2 = @(x,a,b) a*b*(x)./(b-2+x));
f1 = @(x,a,b) a*b*(1-x)./(2*b-x);
rs = @(x1,x2,a,b,p) x1.*f1(x,a,b) - p.*(1-x1).*f2(x,a,b);
R = rs(x1,x2,1,2, p(p));
contour(x1,x2,R);
hold on
end
hold off
But, instead of the required ourput of 10 plots in one graph, I get many plots (>10) in my output. Can someone please help me with this?
  1 个评论
Cris LaPierre
Cris LaPierre 2022-9-21
编辑:Cris LaPierre 2022-9-21
You have not copied your code correctly here, as both examples do not run.

请先登录,再进行评论。

采纳的回答

Cris LaPierre
Cris LaPierre 2022-9-21
I'm assuming your goal is to plot a line for when the function rs is equal to zero for various values of p. You could do that using contour, but I don't think that is the best way. However, you could still use contour to get the (x,y) coordinates of the 0 contour line for each, and then plot those.
np = 10;
[x1,x2] = meshgrid(0:0.001:1,0:0.002:1);
for p = 1:np
f2 = @(x,a,b) a*b*(x)./(b-2+x);
f1 = @(x,a,b) a*b*(1-x)./(2*b-x);
rs = @(xi,x,a,b,p) xi.*f1(x,a,b) - p.*(1-xi).*f2(x,a,b);
R = rs(x1,x2,1,2, p);
figure(1)
M=contour(x1,x2,R,[0 0]);
figure(2)
plot(M(1,2:end),M(2,2:end),'DisplayName',"p="+p)
axis([0 1 0 1])
hold on
end
hold off
legend(Location="northwest")
  4 个评论
Cris LaPierre
Cris LaPierre 2022-9-22
I would keep the meshgrid approach. That will be faster than looping.
I'm not sure of a best way, but one way to do this without contour is to look for the min(abs(R)) in each row. This won't be the exact zero point, but will be close enough for the visualization.
np = 10;
X = 0.001:0.001:1;
Y = 0.002:0.002:1;
[x1,x2] = meshgrid(X,Y);
for p = 1:np
f2 = @(x,a,b) a*b*(x)./(b-2+x);
f1 = @(x,a,b) a*b*(1-x)./(2*b-x);
rs = @(xi,x,a,b,p) xi.*f1(x,a,b) - p.*(1-xi).*f2(x,a,b);
R = rs(x1,x2,1,2, p);
[xval,xind] = min(abs(R),[],2);
X_z0(:,p) = X(xind);
end
plot(X_z0,Y)
axis([0 1 0 1])
legend(Location="northwest")

请先登录,再进行评论。

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 2-D and 3-D Plots 的更多信息

产品


版本

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by