I'm trying to compute the potential and electric field of a uniformly charged spherical shell and plot the results in the space outside the shell using MATLAB. I've tried everything from for loops to sum commands and more and I just can't even get the right expression for the double integral (because I can calculate that my result gives solutions that don't match analytical solutions).

From my research, it seems to best way to plot this would be to plot isosurface plots for different isovalues on the same figure to show how the potential and electric field change with distance from the sphere. But that is my second problem. I can't even compute the double integral properly. So, my question is twofold:

  1. how can I compute a double integral using Riemann sums that contain symbolic variables acting as constants? I am purposely avoiding using a double integral command such as integral2() or a Reimann command like trapz.

  2. how can I plot the resulting expression in 3 dimensions considering I would have a function of three variables?

The current state of my code to calculate the potential (which I define first and then later will take the gradient to get the electric field) is as follows:

syms r phi theta % symbolic values representing a random point in space 
outside a uniformly charged spherical shell

% constants
N = 3; % total number of increments to divide the riemann sum into
Q = +20 .* 1e-9; % total charge on sphere
R = 1; % radius of sphere is 1 meter
sigma = Q/(4*pi*(R^2)); % uniform surface charge denisty
eps0 = 8.854e-12;
kC = 1/(4*pi*eps0); % coulomb's constant
dphi = 2*pi/N; % discretizing the interval over which I sum
dtheta = pi/N;

% Performing two riemann sums
% dA = R^2.*cos(thetaprime).*dphi.*dtheta; small area on surface of sphere
% in spherical coordinates
phiprime = linspace(0,2*pi,N);
thetaprime = linspace(0,pi,N);
dInt1 = sym(zeros(0,N)); % preallocating space allowing the presence of 
dInt2 = sym(zeros(0,N));

for e = 1: length(phiprime)
for m = 1: length(thetaprime)
  dInt1(m) = dInt1 + 
((R^2).*cos(thetaprime).*dphi.*dtheta)./sqrt((r.*sin(theta).*cos(phi) - 
R.*sin(thetaprime).*cos(phiprime)).^2 + (r.*sin(theta).*sin(phi) - 
R.*sin(thetaprime).*sin(phiprime)).^2 + (r.*cos(theta) - 

  dInt2(e) = dInt2 + sum(dInt1);

Int = sum(Int2);

Vtot = kC*sigma.*Int; % total symbolic expression for electric potential

%% Computing Values for Vtot at various points in space and plotting them

rinterval = linspace(1,4,12);
phiinterval = linspace (0,2*pi,12);
thetainterval = linspace(0,pi,12);
[X,Y,Z] = meshgrid(rinterval,phiinterval,thetainterval);

Vcont = subs(Vtot, [r,phi,theta], {X,Y,Z}); % substituting 2D coordinates 
now for symbolic values

For this code I am getting the error:

Error using symengine
Array sizes must match.

Error in sym/privBinaryOp (line 1022)
        Csym = mupadmex(op,args{1}.s, args{2}.s, varargin{:});

Error in  +  (line 7)
X = privBinaryOp(A, B, 'symobj::zipWithImplicitExpansion', 

Error in Numerical_Integration_of_2D_Surfaces (line 27)
     dInt1(m) = dInt1 +    
     ((R^2).*cos(thetaprime).*dphi.*dtheta)./sqrt((r.*sin(theta).*cos(phi) -
     R.*sin(thetaprime).*cos(phiprime)).^2 + (r.*sin(theta).*sin(phi) -
     R.*sin(thetaprime).*sin(phiprime)).^2 + (r.*cos(theta) - 
  • $\begingroup$ Comments are not for extended discussion; this conversation has been moved to chat. $\endgroup$ – nicoguaro Jun 20 at 16:14

Your Answer

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

Browse other questions tagged or ask your own question.