Help with creating contour plot

I am trying to graph the electric potential lines in a plane perpendicular to a ring of charge in the x-y plane. Deriving the formula for off axes points is tedious so I set up a the integrand and then summed it up over a discretized interval in the form of a reimann sum using a for loop. I know my expression for V (Vtot) is correct because it approximately matches analytical solutions for V and -grad(Vtot) appears as it should in the quiver plot to the right below. The problem is probably somewhere in my substitution for the symbols x,y and z (perhaps the meshgrid).
%% Computing a symbolic expression for V for anywhere in space
syms x y z % phiprime is angle that an elemental dq of the circular charge is located at, x,y and z are arbitrary
% points in space outside the charge distribution
N = 200; % number of increments to sum
R = 2; % radius of circle is 2 meters
dphi = 2*pi/N; % discretizing the circular line of charge which spans 2pi
integrand = 0;
for phiprime = 0:dphi:2*pi %phiprime ranges from 0 to 2pi in increments of dphi
integrand = integrand + dphi./(sqrt(((x - R.*cos(phiprime) )).^2 + ((y - R.*sin(phiprime) ).^2) + z.^2));
end
integral = sum(integrand); % uncessary but harmless step that I leave to show that I am using the
% sum of the above expression for each dphi
eps0 = 8.854e-12;
kC = 1/(4*pi*eps0);
rhol = 1*10^-9; % linear charge density
Vtot = kC*rhol*R.*integral; % symbolic expression for Vtot
%% Graphing V & E in plane perpedicular to ring & passing through center
[Y1, Z1] = meshgrid(-4:.5:4, -4:.5:4);
Vcont1 = subs(Vtot, [x,y,z], {0,Y1,Z1}); % Vcont1 stands for V contour; 1 is becuase I do the plane of the ring next
contour(Y1,Z1,Vcont1)
xlabel('y - axis [m]')
ylabel('z - axis [m]')
title('V in a plane perpedicular to a ring of charge (N = 200)')
str = {'Red line is side view', 'of ring of charge'};
text(-1,2,str)
hold on
circle = rectangle('Position',[-2 0 4 .1],'Curvature',[1,1]); % visually displaying line of charge on plot
set(circle,'FaceColor',[1, 0, 0],'EdgeColor',[1, 0, 0]);
g = gradient(-1.*(kC*rhol*R.*integral),[x,y,z]); % taking negative gradient of V and finding symbolic equations
% for Ex, Ey and Ez
% now substituting all the values of the 2D coordinate system for the symbolic x and y variables to get
% numeric values for Ex and Ey because the gradient command can't differentiate a scalar
Ey1 = subs(g(2), [x y z], {0,Y1,Z1});
Ez1 = subs(g(3), [x y z], {0,Y1,Z1});
E1 = sqrt(Ey1.^2 + Ez1.^2); % full numeric magnitude of E in y-z plane
Eynorm1 = Ey1./E1; % This normalizes the electric field lines
Eznorm1 = Ez1./E1;
quiver(Y1,Z1,Eynorm1,Eznorm1);
hold off
***NOTE FOR BELOW FIGURES: The axes are labeled wrong. My mistake! ***** The x axis is (y) and the y axis is (z) ******
As you can see, my contour plot is not filling space. When I lower the value of N, suddenly my contour plot begins to fill space as seen below which makes no sense. N should be totally unrelated; it only determines the accuracy of Vtot (NOTE: I changed the increments of the meshgrid of the plot below, that's why the vectors are less dense).
I've used the method of symbols and a for loop and then substitution for a LINE of charge and everything worked perfectly. I don't know why this contour plot is being uncooperative for a ring of charge. Thanks in advance!

8 Comments

Read about contour you can specify the levels you want.
Andrew
Andrew on 7 Jun 2019
Edited: Andrew on 8 Jun 2019
Thanks for responding!
I'm aware of the levels option for contour plots. It's been suggested by one other person already and I'm surprised people suggest it because levels only makes the contour plot more dense with lines, it doesn't plot lines further out. Essentially, you will see the same small diamond shape I posted but the diamond will be more dense with lines. My problem is that the '*subs*' command doesnt seem to actually be substituting all of the 2D meshgrid for Vtot
Thanks again!
"because levels only makes the contour plot more dense with lines, it doesn't plot lines further out."
contour() without being passed levels determines the levels dynamically according to the range of data. Passing levels does not necessarily make the contour plot more dense with lines: it permits you to control which lines are plotted. Suppose for example that there is an inf in the data, then without being passed levels, linspace(0,inf,10) puts the levels at 0 Inf Inf Inf Inf Inf Inf Inf Inf Inf which is not useful.
Just to show you, when I specify levels as say 15, I get this:
Zoomed in, the left "diamond" looks like this:
I believe there might be inf in the data as the denominator is:
Where the apostrophi is my "prime" notation as in my code:
phiprime
and dϕ' in my derivation is represented in the code as:
dphi
I'm only pointing this out so that hopefully you can help me avoid the inf values somehow. Obnviously, these inf would occur at points (-2,0) and (2,0). I'm almost certain you're right about the problem being related to inf values. I guess my question is now, how can I create a contour plot that includes all of the 2D grid and near but not including points where the the function explodes to infinity?
The problem I see with your example that would show specifying levels is necessary is that you have a linspace which goes to inf. My plot goes a finite distance from -2 to 2 along both dimensions with specified steps that should be plugged into Vtot. It seems these steps aren't being plugged in and the plot is stuck at points near infinity despite me telling it to use only -2:.5:2 . It's as if it's only using -2.5:.5:-1.5
Sorry if I misunderstood you. Thank you for helping
vlevel = linspace(20, 65, 10);
contour(Y1,Z1,Vcont1,vlevel)
Andrew
Andrew on 8 Jun 2019
Edited: Andrew on 8 Jun 2019
OOOOOOHHHHH. I understand what you mean now. I didn't realize you could pick out which lines were plotted by using another vector defined as a linspace and plug that vector into the 'levels' spot of the contour() command.
Now that I re-read the section on 'contour' I realize what "To draw the contour lines at specific heights, specify levels as a vector of monotonically increasing values" means.
Wow you don't how much you've helped me. I've spent hours trying to figure out ways around this and the solution was right under my nose. Thank you Walter. And thanks KSSV. Really appreciate and wish there was some way to repay you besides accepting the answer (which I don't know how to do by the way, sorry I'm new).
Not that you care but here is my beautiful contour plot now:That's exactly what I needed.
I don't know if I should ask this or create a new question but,
is there a way I could make a contour plot of a riemann sum WITHOUT using syms and the 'subs' command to sub values for the variables which are acting as constants in the riemann sum?
Essentially, what I already did here but a method which avoids the 'subs' command since it's my understanding that 'subs' takes a long time to execute for many values.
Not at the moment, but you can improve performance with
Vcont1 = double( subs( vpa(Vtot), [x,y,z], {0,Y1,Z1}) );
I've never heard of Variable-precision arithmetic. After some cursory research, it seems very interesting and useful; and you expressed it in terms of my code. Wow. THANK YOU. Thank you so much. Your responses were prompt (not that they needed to be) and very helpful. Thanks and have a nice day!

Sign in to comment.

Answers (0)

Asked:

on 4 Jun 2019

Edited:

on 8 Jun 2019

Community Treasure Hunt

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

Start Hunting!