How would I manually perform an integral in MATLAB?

My professor has told us not to use the integral3 function in matlab to solve a triple integral. He explained to me today that I should use a nested loop in order to do what we normally would do with an integral. I can't seem to figure out how to do this in matlab. Below I've provided the integral in question that I'm trying to solve as well as my code that I have currently. When I run the code below I get back 18y+54z+6 when the answer should be 16.

∭∇·F dV where dV is dx*dy*dz and the bounds for all three are from -1 to 1.

syms x y z

f = 2*x+y^2+z^2;

fx = diff(f,x)
fy = diff(f,y)
fz = diff(f,z)

F = fx+fy+fz

sum = 0;
dx = 1;
for x = -1:dx:1
    sum = sum + fx*dx
    for y = -1:dx:1
        sum = sum + fy*dx
        for z = -1:dx:1
            sum = sum+fz*dx
        end
    end
end

Thank you for any help you might be able to provide!

[Walter Roberson]

Please do not name any variable "sum". You interfere with the use of the important MATLAB function "sum" when you do that, and even if you are consistent about it in your code you are going to confuse other people who read the code. Remember, your assignment will be marked so other people need to read your code, and if they find it confusing you aren't going to get the best marks.

[William Downing]

Thanks for the tip. I hadn't realized it was a matlab function already.

Answer 1

[Mehdi M]

Hi William, I think the function F should probably be like this. F= (2*x) i + (y^2) j + (z^2) k; So ∇·F= 2+2y+2z. now we integrate numerically, that is

clc
clear
h=.005;
dx=h;
dy=h;
dz=h;
x=-1:dx:1;
y=-1:dy:1;
z=-1:dz:1;
dv=dx*dy*dz;

m=0;
I=0;
for i=1:length(x)
    for j=1:length(y)
        for k=1:length(z)

              m=m+1;
              I=I+(2+2*y(j)+2*z(k))*dv;

       end
    end
 end
 I

To gain more accurate results you can decrease the value of "h"!

[William Downing]

So, how would I go about performing divergence with the way you have the code set up?

    for i = 1:3
        for j = 1:3
            for k = 1:3
                f = divergence((2*x)*i + (y.^2)*j + (z.^2)*k,x,y,z);
                %div = divergence(f,x,y,z);
                sum = sum+div*dv;
            end
        end
    end

Thanks again for the help, I'm trying to figure out how to write it so that I'm not just copying your bit of code without understanding what is going on.