Contour Graph/Plot displays edgy function

Hi all,

I am pretty new to MATLAB and I tried many many times to search online for the answer, but I probably failed to formulate my reserch properly.

I am trying to contour plot the equipotentials and streamlines of a doublet flow for my hydrodynamics course. But I cannot undersatnd why every time I run the script, the function is displayed sort of "edgy". And if I increase the steps from like [-10:1:10] to [-10:0.1:10] in the axis paramenters it just shows me a very tiny version (still edgy) of what I was trying to plot before.

Here's my script:

x = -10:1:10;
y = -10:1:10;
[X,Y] = meshgrid(x,y);
mu = 0.001;
phi = mu*(X./(X.^2+Y.^2));
figure
contour(X,Y,phi, '--r', 'DisplayName', 'Equipotentials')
hold on
psi = mu*(Y./(X.^2+Y.^2));
contour(X,Y,psi, '-b', 'DisplayName', 'Streamlines')
legend
title('Plot of Equipotentials and Streamlines of a Doublet Flow')
xlabel('x-axis')
ylabel('y-axis')

As you can see the displayed funtion isn't really linear, could someone help?

Thank you so much in advance and apologize for my ignorance.

PS: also, is there a way to show the direction of the flow (arrows to show in which direction it's rotating)?

Answer 1

[Dave B]

It looks like it's getting smaller, but actually MATLAB is just picking 10 (different) linearly space levels. The increased resolution is changing (I think) both the min and max of phi and psi, so 10 levels are closer to the middle. You could either add more levels, or be explicit about where you want them:

x = -10:.1:10;
y = -10:.1:10;
[X,Y] = meshgrid(x,y);
mu = 0.001;
phi = mu*(X./(X.^2+Y.^2));
figure
lvls = (-1:.2:1)/1e3; % these are the levels that the -10:1:10 contour picked

contour(X,Y,phi,lvls, '--r', 'DisplayName', 'Equipotentials')
hold on
psi = mu*(Y./(X.^2+Y.^2));
contour(X,Y,psi,lvls, '-b', 'DisplayName', 'Streamlines')
legend
title('Plot of Equipotentials and Streamlines of a Doublet Flow')
xlabel('x-axis')
ylabel('y-axis')

Answer 2

[VBBV]

x = -10:.1:10;
y = -10:.1:10;
[X,Y] = meshgrid(x,y);
mu = 0.1;
phi = mu*(X./(X.^2+Y.^2));
[U,V] = gradient(phi,0.1,0.1);
figure
contour(X,Y,phi,100, '--r', 'DisplayName', 'Equipotentials');
psi = mu*(Y./(X.^2+Y.^2));
[A,B] = gradient(psi,0.1,0.1);
hold on
contour(X,Y,psi,100, '-b', 'DisplayName', 'Streamlines')
hold on
l = streamslice(X,Y,U,V);
m = streamslice(X,Y,A,B);
set(m,'Color','r')

Use the streamslice function to get the flow direction

[Cesare Primultini]

Thank you so much! I really did not know which answer to accept, they were all helpful

Answer 3

[Sulaymon Eshkabilov]

The answer to your posed question lies in the formulation of phi and psi. By decreassing the step size by a factor of 10 for instance the max and min values of phi and psi are also changing by the same factor. Thus, in order to better visualize the formulations, you may consider axis limits, e.g.:

axis([xmin xmax ymin ymax]) % for contour

You can also try to plot it using meshc(), e.g.:

meshc(X,Y, phi), colorbar, colormap jet;
axis([xmin xmax ymin ymax zmin zmax])

Answer 4

[Chris]

phi goes to +/- infinity near the origin.

x = -10:1:10;
y = -10:1:10;
[X,Y] = meshgrid(x,y);
mu = 0.001;
phi = mu*(X./(X.^2+Y.^2));
max(phi(:))

ans = 1.0000e-03

figure
surf(X,Y,phi)
view(-10,20)

At lower resolutions, the maximum calculated distance from 0 is not that high. In the above plot, it's 0.001. As you increase the number of points, points get closer to 0 and the maximum magnitude on the z axis increases. Below, it's 0.01. (I've kept the z limits the same in this plot, but you can see the surface extends far beyond the limits)

x = -10:.1:10;
y = -10:.1:10;
[X,Y] = meshgrid(x,y);
mu = 0.001;
phi = mu*(X./(X.^2+Y.^2));
max(phi(:))

ans = 0.0100

figure
surf(X,Y,phi)
zlim([-.001,.001])
view(-10,20)

contour is scaling its levels based on those limits. One way to address this is by increasing the number of levels. Edit: or use a vector, as Dave suggested.

x = -10:.1:10;
y = -10:.1:10;
[X,Y] = meshgrid(x,y);
mu = 0.001;
phi = mu*(X./(X.^2+Y.^2));
figure
contour(X,Y,phi,50, '--r', 'DisplayName', 'Equipotentials')
hold on
psi = mu*(Y./(X.^2+Y.^2));
contour(X,Y,psi,50, '-b', 'DisplayName', 'Streamlines')
legend
title('Plot of Equipotentials and Streamlines of a Doublet Flow')
xlabel('x-axis')
ylabel('y-axis')

[Cesare Primultini]

Thanks for the detailed explanation! I’ll play around with matlab to understand better this matter