Charge Densities

Triple integrals are used for densities other than mass densities.  In electricity and magnetism, we often consider a "charge cloud" of positively and/or negatively charged particles.  Because the number of charges being considered is very large, counting or otherwise keeping track of individual charges is not practical.

For example, let us suppose that we place a rather high number of "charges" into a box.  Clearly, as you can see below, counting the number of "charges" in the box is not practical.

[Maple Plot]

(commands used to generate figure above)

Instead, we consider a charge density  that describes the distribution of charges in the charge cloud per unit volume.  Often, a charge density is defined

rho(x,y,z)  = number of charges in Coulombs per unit volume (in cubic meters)

Let's look at an example.  The charge cloud above corresponds to the charge density

rho(x,y,z) = 5000*(x+y)  "charges" per unit volume

over box between z= 0 and z= 2 over the unit square.  Let's look at it again.

>    chargerand:= rand(1..100)/100:
for i from 1 to 10000 by 1 do
   ppp[i]:=[1-chargerand()^2,1-chargerand()^2,2*chargerand()]:
end do:
pcloud:=pointplot3d({seq(ppp[i],i=1..10000)},symbol=point,axes=boxed,scaling=constrained,color=blue):
pcloud;

>   

The charge density represents the number of "charges" per unit of volume near a given point. That is, the small number of charges dQ  in a small "box" with width dx , length dy , and height dz is given by

dQ = rho(x,y,z)*dx*dy*dz

Let's illustrate this fact below:

>    xmin:=0:
xmax:=1:
ymin:=0:
ymax:=1:
ztop:=2:
zbot:=0:

n:=10: #total number of subintervals in x direction
m:=10: #total number of subintervals in y direction
r:=10:  #total number of subintervals in z direction
x_subinterval:=5: #enter an integer between 1 and n
y_subinterval:=4: #enter an integer between 1 and m
z_subinterval:=5: #enter an integer between 1 and r

# Calculate point, box, surface
deltx:=(xmax-xmin)/n:
x_at := xmin+deltx*x_subinterval:
delty:=(ymax-ymin)/m:
y_at := ymin+delty*y_subinterval:
deltz:= abs(subs(x=x_at,y=y_at,ztop)-subs(x=x_at,y=y_at,zbot) )/r:
z_at := subs(x=x_at,y=y_at,zbot)+deltz*z_subinterval:

U_at:=subs(x=x_at-deltx/2,y=y_at-delty/2,z=z_at-deltz/2,U):

#surfaces
p1:=plot3d({ztop,zbot},x=xmin..xmax,y=ymin..ymax,
    color=magenta,style=wireframe,axes=boxed,grid=[n+1,m+1]):
p2:=plot3d([xmin,y,v],y=ymin..ymax,v=subs(x=xmin,zbot)..subs(x=xmin,ztop),
    color=blue,style=wireframe,grid=[m+1,r+1]):
p3:=plot3d([xmax,y,v],y=ymin..ymax,v=subs(x=xmax,zbot)..subs(x=xmax,ztop),
    color=blue,style=wireframe,grid=[m+1,r+1]):
p4:=plot3d([x,ymin,v],x=xmin..xmax,v=subs(y=ymin,zbot)..subs(y=ymin,ztop),
    color=green,style=wireframe,grid=[m+1,r+1]):
p5:=plot3d([x,ymax,v],x=xmin..xmax,v=subs(y=ymax,zbot)..subs(y=ymax,ztop),
    color=green,style=wireframe,grid=[m+1,r+1]):

#singlebox
side1:=polygon([[x_at,y_at,z_at-deltz],[x_at,y_at,z_at],[x_at,y_at-delty,z_at],
                [x_at,y_at-delty,z_at-deltz]],color=magenta):
side2:=polygon([[x_at-deltx,y_at,z_at-deltz],[x_at-deltx,y_at,z_at],
                [x_at-deltx,y_at-delty,z_at],
                [x_at-deltx,y_at-delty,z_at-deltz]],color=magenta):
side3:=polygon([[x_at,y_at,z_at-deltz],[x_at,y_at,z_at],[x_at-deltx,y_at,z_at],
                [x_at-deltx,y_at,z_at-deltz]],color=magenta):
side4:=polygon([[x_at,y_at-delty,z_at-deltz],[x_at,y_at-delty,z_at],                                 [x_at-deltx,y_at-delty,z_at],
                [x_at-deltx,y_at-delty,z_at-deltz]],color=magenta):
atop:= polygon([[x_at,y_at,z_at],[x_at-deltx,y_at,z_at],[x_at-deltx,y_at-delty,z_at],
                [x_at,y_at-delty,z_at]],color=magenta):
bot:=  polygon([[x_at,y_at,z_at-deltz],[x_at-deltx,y_at,z_at-deltz],
                [x_at-deltx,y_at-delty,z_at-deltz],
                [x_at,y_at-delty,z_at-deltz]],color=magenta):

t1:=textplot3d([1.01*x_at,1.01*y_at,z_at-deltz/2,"dQ = rho(x,y,z)dxdydz"],color=black,
               font=[TIMES,ROMAN,14],align=RIGHT):

display(side1,side2,side3,side4,atop,bot,p1,p2,p3,p4,p5,t1,pcloud,orientation=[-45,60]);

>   

If we now partition the entire box and calculate the limit as the volume of the boxes approaches 0, then the result will be an exact estimate of the number of "charges" in the container.

>    #surfaces
   p1:=plot3d({ztop,zbot},x=xmin..xmax,y=ymin..ymax,
       color=grey,style=wireframe,axes=boxed,grid=[15,15]):
   p2:=plot3d([xmin,y,v],y=ymin..ymax,v=subs(x=xmin,zbot)..subs(x=xmin,ztop),
       color=grey,style=wireframe,grid=[15,15]):
   p3:=plot3d([xmax,y,v],y=ymin..ymax,v=subs(x=xmax,zbot)..subs(x=xmax,ztop),
       color=grey,style=wireframe,grid=[15,15]):
   p4:=plot3d([x,ymin,v],x=xmin..xmax,v=subs(y=ymin,zbot)..subs(y=ymin,ztop),
       color=grey,style=wireframe,grid=[15,15]):
   p5:=plot3d([x,ymax,v],x=xmin..xmax,v=subs(y=ymax,zbot)..subs(y=ymax,ztop),
       color=grey,style=wireframe,grid=[15,15]):

jj:=0:
for n from 2 to 6 by 2 do
   ii:=0:
  
   m := n:
   r := n:

   for j from 1 to n do
     deltx:=(xmax-xmin)/n:
     x_at := xmin+j*deltx:
     for k from 1 to m do
        delty:=(ymax-ymin)/m:
        y_at := ymin+k*delty:
        deltz:= abs(subs(x=x_at,y=y_at,ztop)-subs(x=x_at,y=y_at,zbot) )/r:
        for l from 1 to r do
            z_at := subs(x=x_at,y=y_at,zbot)+l*deltz:

   #singlebox
   side1:=polygon([[x_at,y_at,z_at-deltz],[x_at,y_at,z_at],[x_at,y_at-delty,z_at],
                [x_at,y_at-delty,z_at-deltz]],color=magenta):
   side2:=polygon([[x_at-deltx,y_at,z_at-deltz],[x_at-deltx,y_at,z_at],
                [x_at-deltx,y_at-delty,z_at],
                [x_at-deltx,y_at-delty,z_at-deltz]],color=magenta):
   side3:=polygon([[x_at,y_at,z_at-deltz],[x_at,y_at,z_at],[x_at-deltx,y_at,z_at],
                [x_at-deltx,y_at,z_at-deltz]],color=magenta):
   side4:=polygon([[x_at,y_at-delty,z_at-deltz],[x_at,y_at-delty,z_at],                                                                                                                                 [x_at-deltx,y_at-delty,z_at],
                [x_at-deltx,y_at-delty,z_at-deltz]],color=magenta):
   atop:= polygon([[x_at,y_at,z_at],[x_at-deltx,y_at,z_at],[x_at-deltx,y_at-delty,z_at],
                [x_at,y_at-delty,z_at]],color=magenta):
   bot:=  polygon([[x_at,y_at,z_at-deltz],[x_at-deltx,y_at,z_at-deltz],
                [x_at-deltx,y_at-delty,z_at-deltz],
                [x_at,y_at-delty,z_at-deltz]],color=magenta):
   ii:=ii+1:
   pp[ii]:=display(side1,side2,side3,side4,atop,bot):
    end do: # l
    end do: # k
    end do: # j
jj:=jj+1:
qq[jj]:=display(seq(pp[i],i=1..ii),insequence=false):
end do:  # n
ani:=display(seq(qq[j],j=1..jj),insequence=true):
display(ani,p1,p2,p3,p4,p5,axes=boxed,pcloud,orientation=[-45,60]);

>   

Letting the volumes of the individual boxes approach 0 implies a triple integral of the charge density over the box.  Thus, the total charge Q  is given by

Q = Limit(sum(sum(sum(rho(s[j],t[k],u[l])*Delta*x[j]*Delta*y[k]*Delta*z[l],l = 1 .. r),k = 1 .. m),j = 1 .. n),h = 0)   =   int(int(int(rho(x,y,z),z = 0 .. 2),y = 0 .. 1),x = 0 .. 1)

Since rho(x,y,z) = 5000*(x+y) , we can calculate this triple integral.

>    Q:=Int(Int(Int(5000*(x+y),z=0..2),y=0..1),x=0..1) = int(int(int(5000*(x+y),z=0..2),y=0..1),x=0..1);

>   

Indeed, we began this section by placing exactly Q  charges into the box.