3
$\begingroup$

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
g = gradient(-1.*(kC*rhol*R.*intgrl),[x,y,z]);

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

Contour Plot

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$.

Electric Field Plot

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: Derivation of V

  • $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
$\endgroup$
2
$\begingroup$

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;

%% Compute the gradient
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.

enter image description here

$\endgroup$
  • $\begingroup$ 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! $\endgroup$ – Andrew Aug 10 at 0:00
  • $\begingroup$ @Andrew, I'm glad that it was useful to you. You can accept the answer if it solves your problem. $\endgroup$ – nicoguaro Aug 10 at 0:04
2
$\begingroup$

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: enter image description here

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.

A link to our full discussion here: https://www.mathworks.com/matlabcentral/answers/465471-help-with-creating-contour-plot?s_tid=prof_contriblnk

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.