Help with creating contour plot

Question

Andrew 2019-6-4

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/465471-help-with-creating-contour-plot

编辑： 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 个评论
显示 6更早的评论隐藏 6更早的评论

Andrew 2019-6-8

编辑：Andrew 2019-6-8

在 MATLAB Online 中打开

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

Andrew 2019-6-8

编辑：Andrew 2019-6-8

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.

Walter Roberson 2019-6-8

在 MATLAB Online 中打开

Not at the moment, but you can improve performance with

Vcont1 = double( subs( vpa(Vtot), [x,y,z], {0,Y1,Z1}) );

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!

请先登录，再进行评论。

请先登录，再回答此问题。

Help with creating contour plot

8 个评论
显示 6更早的评论隐藏 6更早的评论

回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

Help with creating contour plot

8 个评论 显示 6更早的评论隐藏 6更早的评论

回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

8 个评论
显示 6更早的评论隐藏 6更早的评论