Help with creating contour plot

7 ビュー (過去 30 日間)
Andrew
Andrew 2019 年 6 月 4 日
編集済み: Andrew 2019 年 6 月 8 日
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 件のコメント
Walter Roberson
Walter Roberson 2019 年 6 月 8 日
Not at the moment, but you can improve performance with
Vcont1 = double( subs( vpa(Vtot), [x,y,z], {0,Y1,Z1}) );
Andrew
Andrew 2019 年 6 月 8 日
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!

サインインしてコメントする。

回答 (0 件)

カテゴリ

Help Center および File ExchangeMatrix Indexing についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by