Problems with quiver plot

Hey!

I try to create a quiver plot with unequl axis length. I' like to have the arrows the same length, which somehow does not want to work.

Any ideas?

Thanks!

%% system paramters
eta = .1;
mu = .1;
nu = 1;
gamma = 2;

%% dependent variables 
roi = 2;

s = linspace(max([(1-roi)*s_0,0]),(1+roi)*s_0,10);

Unrecognized function or variable 's_0'.

p = linspace(max([(1-roi)*p_0,0]),(1+roi)*p_0,10);

[s,p] = meshgrid(s,p);
%% gradient flow
v = s.*p.^gamma ./ (1+(1+s).*p.^gamma);

ds = -v + eta;
dp = mu*(v - nu*p);

mag = sqrt(ds.^2 + dp.^2);
arrow_scale = 3E-1;
norm_ds = arrow_scale*ds./mag;
norm_dp = arrow_scale*dp./mag;
%% plot
q = quiver(s,p,norm_ds,norm_dp,'Autoscale','off', 'Color',.6*[1,1,1]);
q.ShowArrowHead = 'off';
q.Marker = '.';

[Matt J]

As shown above, your code doesn't provide all the necessary inputs to run.

[Sam Chak]

@Lukas, Do you truly want the arrows to have the same length? The relative lengths of the arrows indicate the magnitudes of the vectors (or the strengths of the gradient flow).

[David Goodmanson]

Often arrows of constant length, sans arrowheads, are used to show a direction field in the phase plane. Less clutter that way, and direction often provides the most crucial information. See e.g. Strogatz, Nonlinear Dynamics and Chaos 2nd ed., figs. 6.1.4 and 6.3.2 or what's quicker, 'slope field wiki '..

[Sam Chak]

Thanks @David Goodmanson. I agree that this approach will reduce clutter. The original values for s_0​ and p_0​ were not provided, and the original image posted was removed by the OP. In the original image, the range for the x-axis is [0, 35], while the range for the y-axis is [0, 0.35]. Thus, I was able to estimate the values for s_0​ and p_0​ that generate a plot very close to the original image.

By default, MATLAB's quiver function behaves like a mathematically-defined vector: the arrowhead indicates direction, and its length indicates relative magnitude. However, I don't know how to set a fixed length for all arrows in the phase plane plot. I attempted to use @Matt J's approach, but regardless of how I scaled the arrows, if scaled a little too large, the quivers on the left and right sides overlap, and if scaled a little too small, the quivers in the center region appear as dots.

Do you know how to maintain a constant length for the arrows as shown in the Slope field wiki?

s_0 = 10;   % estimated based on the original image posted by the OP (now removed)
p_0 = 0.1;  % estimated based on the original image posted by the OP (now removed)

%% system paramters
eta   = .1;
mu    = .1;
nu    = 1;
gamma = 2;

%% dependent variables 
roi = 2;

numArrX = 14;    % number of arrows per row     (OP used 10, which I find the grid relatively coarse)
numArrY = 14;    % number of arrows per column  (OP used 10, which I find the grid relatively coarse)
s = linspace(max([(1-roi)*s_0,0]), (1+roi)*s_0, numArrX);
p = linspace(max([(1-roi)*p_0,0]), (1+roi)*p_0, numArrY);

[s,p] = meshgrid(s,p);

%% gradient flow
v   = s.*p.^gamma./(1 + (1 + s).*p.^gamma);
ds  = -v + eta;
dp  = mu*(v - nu*p);

mag = sqrt(ds.^2 + dp.^2);
arrow_scale = 5E-1;     % <---- scaling proposed by Matt J
norm_ds = arrow_scale*ds./mag;
norm_dp = arrow_scale*dp./mag;

% Finding the equilibrium point
fun = @(x) [-(x(1).*x(2).^gamma./(1 + (1 + x(1)).*x(2).^gamma)) + eta;
            mu*(x(1).*x(2).^gamma./(1 + (1 + x(1)).*x(2).^gamma) - nu*x(2))];
x0  = [11, 1];          % initial guess
eq  = fsolve(fun, x0)   % equilibrium point

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

eq = 1×2

11.2222 0.1000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

%% plot
% streamslice(s, p, norm_ds, norm_dp);
q = quiver(s, p, norm_ds, norm_dp, 'off', 'Color', .6*[1,1,1]); % automatic scaling is disabled
q.ShowArrowHead = 'off';    % no arrowheads
q.Marker = '.';             % for the tails

% adding the equilibrium point to the slope field
hold on
plot(eq(1), eq(2), 'ro', 'markersize', 10, 'linewidth', 1.5)
hold off

xlabel('s')
ylabel('p')
xlim([0 30])
ylim([0 .3])

[David Goodmanson]

Hi Sam, here's the best I could do trying to reproduce the wikipedia plot, including the aspect ratio. The arrows are all normalized to the same value, but at least as importantly the plottting points for quiver are not equally spaced meshgrid values. Rather every x,y quiver point is changed slightly from what meshgrid has. I didn't use 'axis equal' so I guess the vectors are not quite constant length visually.

xx = linspace(-5,5,22);
yy = linspace(-10,10,22);
[x0 y0] = meshgrid(xx,yy);  
th = atan(x0.^2-x0-2);
sf = 1/2;        % factor to visually reduce the arrow length on the plot
u = sf*cos(th);
v = sf*sin(th);
x = x0 - u/2;    % move the center of the arrow to the constant-spaced points
y = y0 - v/2;
figure(1)
quiver(x,y,u,v,'showarrowhead','off','autoscale','off')
ylim([-10 10])
xlim([-10 10])
hold on
x1 = -5:.01:5
y1 = x1.^3/3 -x1.^2/2-2*x1+4;
y2 = x1.^3/3 -x1.^2/2-2*x1;
y3 = x1.^3/3 -x1.^2/2-2*x1-4;
plot(x1,y1,x1,y2,x1,y3)
hold off

[Sam Chak]

Hi @David Goodmanson

Thank you for your input. It appears that there is no specific parameter to set a constant length for all quiver objects without altering the original magnitudes of the directional components specified by u and v. However, you are absolutely correct that "constant length" representations are visually meaningless if the aspect ratios of the x- and y-axes are not equal.

s_0 = 10;   % estimated based on the original image posted by the OP (now removed)
p_0 = 0.1;  % estimated based on the original image posted by the OP (now removed)

%% system paramters
eta   = .1;
mu    = .1;
nu    = 1;
gamma = 2;

%% dependent variables 
roi = 2;

numArrX = 19;    % number of arrows per row
numArrY = 19;    % number of arrows per column
s = linspace(max([(1-roi)*s_0,0]), (1+roi)*s_0, numArrX);
p = linspace(max([(1-roi)*p_0,0]), (1+roi)*p_0, numArrY);

[s,p] = meshgrid(s,p);

%% gradient flow
v   = s.*p.^gamma./(1 + (1 + s).*p.^gamma);
ds  = -v + eta;
dp  = mu*(v - nu*p);

mag          = sqrt(ds.^2 + dp.^2);
Xarrow_scale = 4E-1;
Yarrow_scale = 1E-1;
norm_ds      = Xarrow_scale*ds./mag;
norm_dp      = Yarrow_scale*dp./mag;

% Finding the equilibrium point
fun = @(x) [-(x(1).*x(2).^gamma./(1 + (1 + x(1)).*x(2).^gamma)) + eta;
            mu*(x(1).*x(2).^gamma./(1 + (1 + x(1)).*x(2).^gamma) - nu*x(2))];
x0  = [11, 1];          % initial guess
eq  = fsolve(fun, x0)   % equilibrium point

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

eq = 1×2

11.2222 0.1000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

%% plot
l = streamslice(s, p, norm_ds, norm_dp, 0.5, 'noarrows');
set(l, 'Color', "#F63C4C");     % Red Salsa
hold on
q = quiver(s, p, norm_ds, norm_dp, 'off', 'Color', .6*[1,1,1]); % automatic scaling is disabled
q.ShowArrowHead = 'off';    % no arrowheads
q.Marker = '.';             % for the tails

% adding the equilibrium point to the slope field
plot(eq(1), eq(2), 'o', 'markersize', 10, 'linewidth', 1.5, 'Color', "#2F2CE0") % Palatinate Blue
hold off

title('Gradient flow')
xlabel('s')
ylabel('p')
xlim([0 30])
ylim([0 .3])

Answer 1

[Matt J]

I think you just need axis equal.

%% system paramters
eta = .1;
mu = .1;
nu = 1;
gamma = 2;

%% dependent variables 
roi = 2;
s_0=1; p_0=3;  %<---- Matt J chose randomly

s = linspace(max([(1-roi)*s_0,0]),(1+roi)*s_0,10);
p = linspace(max([(1-roi)*p_0,0]),(1+roi)*p_0,10);

[s,p] = meshgrid(s,p);
%% gradient flow
v = s.*p.^gamma ./ (1+(1+s).*p.^gamma);

ds = -v + eta;
dp = mu*(v - nu*p);

mag = sqrt(ds.^2 + dp.^2);
arrow_scale = 3E-1;
norm_ds = arrow_scale*ds./mag;
norm_dp = arrow_scale*dp./mag;
%% plot
q = quiver(s,p,norm_ds,norm_dp,'Autoscale','off', 'Color',.6*[1,1,1]);
q.ShowArrowHead = 'off';
q.Marker = '.'; 
axis equal %<---- Matt J added

●

[Sam Chak]

I believe that the missing values are

s_0 = 10;
p_0 = 0.1;

[Lukas]

Thanks for this suggestions. s_0 and p_0 have indeed a quite large aspect ration (ca. 100), which is why equal axis does not work to well...

[Matt J]

s_0 and p_0 have indeed a quite large aspect ration (ca. 100),

And yet s_0 and p_0 have the same units? If not, it makes more sense if they were made to...

Alternatively, I guess you could adjust the U,V data for the current aspect ratio. If you change the aspect ratio, though, this will have to be repeated.

%% system paramters
eta = .1;
mu = .1;
nu = 1;
gamma = 2;

%% dependent variables 
roi = 2;
s_0=0.1; p_0=10;  %<---- Matt J added

s = linspace(max([(1-roi)*s_0,0]),(1+roi)*s_0,10);
p = linspace(max([(1-roi)*p_0,0]),(1+roi)*p_0,10);

[s,p] = meshgrid(s,p);
%% gradient flow
v = s.*p.^gamma ./ (1+(1+s).*p.^gamma);

ds = -v + eta;
dp = mu*(v - nu*p);

mag = sqrt(ds.^2 + dp.^2);
arrow_scale = 2E-2; %<---- Matt J changed
norm_ds = arrow_scale*ds./mag;
norm_dp = arrow_scale*dp./mag;
%% plot
q = quiver(s,p,norm_ds,norm_dp,'off', 'Color',.6*[1,1,1]);
q.ShowArrowHead = 'off';
q.Marker = '.'; 
d=daspect; q.VData=q.VData*d(2);   %<---- Matt J added

●

[Lukas]

Perffect, this works. Many thanks!

[Matt J]

You're very welcome, but since it works, please Accept-click the answer to indicate so.