Help with creating contour plot
7 ビュー (過去 30 日間)
古いコメントを表示
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
2019 年 6 月 8 日
Not at the moment, but you can improve performance with
Vcont1 = double( subs( vpa(Vtot), [x,y,z], {0,Y1,Z1}) );
回答 (0 件)
参考
カテゴリ
Help Center および File Exchange で Matrix Indexing についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!