How to compute a double Integral using two Riemann sums and graph f(x,y,z) = c
조회 수: 2 (최근 30일)
이전 댓글 표시
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 doesn't give solutions that match analytical solutions).
Basically, I'm trying to compute the potential and field of a uniformly charged spherical shell and plot the results in the space outside the shell. From my research, it seems to best way to plot this would be to plot isosurface plots for different isovalues (graph f(x,y,z) = c for multiple values of c) 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?
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 symbols
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) - R.*cos(phiprime)).^2)
end
dInt2(e) = dInt2 + sum(dInt1);
end
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
Right now I am currently 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', '_plus');
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) - R.*cos(phiprime)).^2)
Thanks in advance!
댓글 수: 0
답변 (0개)
참고 항목
카테고리
Help Center 및 File Exchange에서 Calculus에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!