Graphing electric potential of a ring of charge using MATLAB help

Here is a summary of what I am trying to do:

Use MATLAB to compute the potential $$V$$ at any point $$(x, y, z)$$ in space due to a uniform ring of charge. Use a Riemann sum to compute the integral with increments, $$N$$, as a variable you can change. Plot your potential and field in the plane perpendicular to the area of the ring and passing through the center.

I have two versions of code that give me the same result: correct expression for $$V$$ ($$V_\text{tot}$$) and the correct vector field. My problem is getting my contour plot to fill the 2D space. The only difference between the two codes is that one uses a for loop to perform the summation and the other uses a sum command:

Version 1 (using for loop):

%% 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

intgrl = sum(integrand);
% unnecessary 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 = 1e-9; % linear charge density

Vtot = kC*rhol*R.*intgrl; % 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 because 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 perpendicular to a ring of charge (N = 200)')
str = {'Red line is side view', 'of ring of charge'};
text(-1,2,str)

hold on
% visually displaying line of charge on plot
circle = rectangle('Position',[-2 0 4 .1],'Curvature',[1,1]);
set(circle,'FaceColor',[1, 0, 0],'EdgeColor',[1, 0, 0]);

% 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
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


Version 2 (using sum):

syms x y z
R = 2; % radius of circle is 2 meters
N=100;
dphi = 2*pi/N; % discretizing the circular line of charge which spans 2pi

phiprime = 0:dphi:2*pi; %phiprime ranges from 0 to 2pi in increments of dphi

integrand = dphi./(sqrt(((x - R.*cos(phiprime) )).^2 + ((y - R.*sin(phiprime) ).^2) + z.^2));

phiprime = 0:dphi:2*pi;
intgrl = sum(integrand); % Riemann sum performed here

eps0 = 8.854e-12;
kC = 1/(4*pi*eps0);
rhol = 1e-9; % linear charge density

Vtot = kC*rhol*R.*intgrl; % symbolic expression for Vtot

%%
[Y1, Z1] = meshgrid(-4:.5:4,-4:.5:4);
Vcont1 = subs(Vtot, [x,y,z], {0,Y1,Z1});

contour(Y1,Z1,Vcont1)
xlabel('y - axis [m]')
ylabel('z - axis [m]')
title('V in a plane perpedicular to a ring of charge (N = 100)')
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 ring of charge on plot
set(circle,'FaceColor',[1, 0, 0],'EdgeColor',[1, 0, 0]);

g = gradient(-1.*(kC*rhol*R.*intgrl),[x,y,z]); % taking negative gradient of V and finding symbolic equations for Ex, Ey and Ez

% 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 doesn't accept symbolic arguments
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


Both versions of code produce the following graphs:

Note: the picture above this text should have the axes be $$y$$ and $$z$$ like the picture below, not $$x$$ and $$y$$. Also, the title of the picture below this text should be "$$E$$ in the plane..." not $$V$$.

As you can see, the vector field is correct while the contour plot seems to use only a few points around the ends of the ring and connect them with straight lines in a strange diamond shape. I can't get it to fill space.

As for my derivation of the formula for $$V$$, it is here:

• $$s$$ is the radial distance
• $$\phi$$ is the azimuthal angle
• An arbitrary point $$P(x,y,z)$$ has no superscript while the prime notation is used to identify a point on the ring of charge ($$\phi^\prime$$, $$x^\prime$$, etc)
• cursive $$r$$ indicates relative distance from a point on the ring to $$P(x,y,z)$$
• $$R$$ is the radius of the ring, I chose 2 meters
• $$\mathrm{d}s = R \mathrm{d}\phi$$ is a small amount of arc length

You don't need symbolic variables to compute the approximated potential for your Riemann sums. You can just use meshgrid to evaluate the potential in each point of interest. For the electric field, you can just compute the derivative analytically and then repeat the process for each component or compute a numerical gradient with gradient.

There are better ways to approximate the integral than just assume that you have straight segments with uniform charge distribution, but they are conceptually similar to this procedure.

The following snippet computes what you want without using symbolic variables.

%% Parameters
R = 2.0;
k = 1.0;
charge_density = 1.0;
ndiv_ring = 100;
ndiv_x = 50;
ndiv_y = 100;
ndiv_z = 100;

%% Computation
phi = linspace(0, 2*pi, ndiv_ring);
dphi = 2*pi/(ndiv_ring - 1);
x_ring = R*cos(phi);
y_ring = R*sin(phi);
x = linspace(0, 4, ndiv_x);
y = linspace(-4, 4, ndiv_y);
z = linspace(-4, 4, ndiv_z);
[X, Y, Z] = meshgrid(x, y, z);
V = zeros(size(X));
for cont = 1:ndiv_ring
Rx = X - R*cos(phi(cont));
Ry = Y - R*sin(phi(cont));
Rmag = sqrt(Rx.^2 + Ry.^2 + Z.^2);
V = V + dphi/Rmag;
end
V = k*R*charge_density*V;

hx = x(2) - x(1);
hy = y(2) - y(1);
hz = z(2) - z(1);
[Ex, Ey, Ez] = gradient(-V, hx, hy, hz);

%% Visualization
levels = 5;
isovalues = linspace(2, 9, levels);
for cont = 1:levels
[faces, verts, colors] = isosurface(X, Y, Z, V, isovalues(cont), V);
p = patch('Vertices', verts, 'Faces', faces, 'FaceVertexCData',...
colors, 'FaceColor','interp','EdgeColor','interp');
isonormals(X, Y, Z, V, p)
p.EdgeColor = 'none';
p.FaceAlpha = 0.4;
end
daspect([1 1 1])
view(3);
axis vis3d
camlight
lighting gouraud


The visualization part generate isosurfaces that are the equivalent to isocontours for 3D data. The resulting image should look like the following.

• I know it's been a long time since we chatted and you may have forgotten by now but I just wanted to say thank for posting the code. I didn't realize you could calculate the numerical gradient at specified divisions. I really appreciate it. Unfortunately had some personal stuff come up in my life and had to take some time off. This code is VERY efficient. You've helped a lot. It took me a while just to understand how the code worked. Thanks! Aug 10 '19 at 0:00
• @Andrew, I'm glad that it was useful to you. You can accept the answer if it solves your problem. Aug 10 '19 at 0:04

I received a solution to this question from MATLAB's community.

Essentially, I need to specify which contour lines to plot using the 'levels' spot on the 'contour()' command.

Levels allows you to not only choose how many but which lines to plot. If you define a vector such as

vlevel = linspace(20, 65, 10);


and then place it in the 'levels' spot of contour:

contour(Y1,Z1,Vcont1,vlevel)


Now the contour command will only plot 10 evenly spaced lines starting at the 20th line and ending at the 65th line (arbitrary parameters which you may choose based on your needs).

My contour plot now looks like:

Which is really beautiful imo. I should add that I increased my plot to show 15 lines between the 20th and the 60th so my image is not exactly representative of my example.