Integrating the volume of a sphere using a step function as an integrand

My task is to calculate the volume of a sphere using the integral3 module in matlab.

My integration is over the cartezian (x,y,z) coordinates of a cube that contains my sphere. To only integrate over the x,y,z coordinates of the sphere inside the cube, I use a step function called "sphere" that returns the integrand 1 if the coordinates are inside or on the sphere, and 0 if they are outside the sphere.

However, I haven't been able to code this into something that does not return an error. Please help!

Here is my step function:

q = integral3(sphere,-1,1,-1,1,-1,1)

function t = sphere(x,y,z)
    r0 = 1;
    r1 = (x.^2+y.^2+z.^2).^0.5;
    if r1 <= r0
        t = 1;
    else 
        t = 0; 
    end
end

However this returns the error:

Not enough input arguments.
Error in untitled>sphere (line 3)
r1 = (x.^2+y.^2+z.^2).^0.5;

I've also tried to use an anonymous function:

fun = @(x,y,z) sphere(x,y,z);
q = integral3(fun,-1,1,-1,1,-1,1)

But this also returns a bunch of errors:

Error using integral2Calc>integral2t/tensor (line 241)
Integrand output size does not match the input size.

Error in integral2Calc>integral2t (line 55)
[Qsub,esub] = tensor(thetaL,thetaR,phiB,phiT);

Error in integral2Calc (line 9)
[q,errbnd] = integral2t(fun,xmin,xmax,ymin,ymax,optionstruct);

Error in integral3>innerintegral (line 128)
Q1 = integral2Calc( ...
    
Error in integral3>@(x)innerintegral(x,fun,yminx,ymaxx,zminxy,zmaxxy,integral2options) (line 111)
f = @(x)innerintegral(x, fun, yminx, ymaxx, ...
    
Error in integralCalc/iterateScalarValued (line 314)
fx = FUN(t);

Error in integralCalc/vadapt (line 132)
[q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 75)
[q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in integral3 (line 113)
Q = integralCalc(f,xmin,xmax,integralOptions);

Answer 1

[Matt J]

You won't be able to get very good accuracy integrating a step function because of its discontinuity. It would be better if you integrated in spherical coordinates where a step function is unnecessary. Regardless, here is how the code should look if you must integrate in Cartesian space.

r0=1;
fun=@(x,y,z)double(sqrt(x.^2+y.^2+z.^2)<=r0);

q = integral3(fun,-r0,r0,-r0,r0,-r0,r0,'AbsTol',1e-2)

q = 4.1887

[Ilan Shumilin]

Thanks for your help!

[Ilan Shumilin]

Actually, Id appreciate some help with my main question, where my volume integration has some meaning. I followed the same method, but now my integrand is a function of (x,y,z). However, my script returns some error about matrixes not being the same size. Do you have any idea what the problem is?

Parameters of the system

k_b = 1.380649*10^(-23); 
h = 6.62607004*10^(-34);
Nav = 6.02214076*10^23;
e = 1.60217662*10^(-19);
ec = 8.8541878176*10^(-12);
T = 298;
C = 1;
di_cof = 81.0-C.*6.6;
eps = di_cof*ec;  
z_ion = 1;
r0 = 7.8*10^(-10); 
x0 = 3*r0;
d = 500;
kappa = (2*e^2*1000*Nav*C/(eps*k_b*T))^(0.5);
constant = z_ion*e/(4*pi*eps*(1+kappa*r0));
dis = 4*10^(-9); %distance between the two spheres

Parameters of the integration

d_x = 100*10.^(-9); %indicates the distance of the integration
d_y = 100*10.^(-9); %indicates the distance of the integration
d_z = 100*10.^(-9); %indicates the distance of the integration

Integration

fun = @(x,y,z)double(and(sqrt((x-x0).^2+y.^2+z.^2)>=r0,sqrt((x+x0).^2+y.^2+z.^2)>=r0).*constant.*((exp(-kappa*(sqrt((x+x0).^2+y.^2+z.^2)-r0))/sqrt((x+x0).^2+y.^2+z.^2) + exp(-kappa*(sqrt((x-x0).^2+y.^2+z.^2)-r0))/sqrt((x-x0).^2+y.^2+z.^2))));
q = integral3(fun,-d_x,d_x,-d_y,d_y,-d_z,d_z,'AbsTol',1e-1)

Error:

Warning: Matrix is singular to working precision.
Warning: Matrix is singular to working precision.
Arrays have incompatible sizes for this operation.

Error in DLVO (line 21)
cof = @(x,y,z)double(and(sqrt((x-x0).^2+y.^2+z.^2)>=r0,sqrt((x+x0).^2+y.^2+z.^2)>=r0).*constant.*((exp(-kappa*(sqrt((x+x0).^2+y.^2+z.^2)-r0))/sqrt((x+x0).^2+y.^2+z.^2) + exp(-kappa*(sqrt((x-x0).^2+y.^2+z.^2)-r0))/sqrt((x-x0).^2+y.^2+z.^2))));

Error in integral3>@(y,z)fun(x(1)*ones(size(z)),y,z) (line 129)
@(y,z)fun(x(1)*ones(size(z)),y,z), ...
    
Error in integral2Calc>integral2t/tensor (line 237)
Z1 = FUN(X(VTSTIDX),Y(VTSTIDX)); NFE = NFE + 1;

Error in integral2Calc>integral2t (line 55)
[Qsub,esub] = tensor(thetaL,thetaR,phiB,phiT);

Error in integral2Calc (line 9)
[q,errbnd] = integral2t(fun,xmin,xmax,ymin,ymax,optionstruct);

Error in integral3>innerintegral (line 128)
Q1 = integral2Calc( ...
    
Error in integral3>@(x)innerintegral(x,fun,yminx,ymaxx,zminxy,zmaxxy,integral2options) (line 111)
f = @(x)innerintegral(x, fun, yminx, ymaxx, ...
    
Error in integralCalc/iterateScalarValued (line 314)
fx = FUN(t);

Error in integralCalc/vadapt (line 132)
[q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 75)
[q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in integral3 (line 113)
Q = integralCalc(f,xmin,xmax,integralOptions);

[Matt J]

You need to be doing element-wise division ./

Everything should be element-wise, for that matter.

[Ilan Shumilin]

Wow! Thanks so much! again!