kaushal.singh_12 asked . 05/08/2018

How can I use Matlab to evaluate ∫10∫10∫101(xyz)xyz dx dy dz?

Matlab

Expert Answer

Kshitij Singh ( Staff ) answered . 06/08/2018

There are many choices! Here are a few:
Naive Numerical Integration

%Areas of rectangles:
f=@(x,y,z) (x.*y.*z).^(-x.*y.*z); %the function you want to integrate
dx = 1e-2; %discretization step
[X,Y,Z]=ndgrid(dx/2:dx:1); %I omit the left endpoints to avoid 0^0 issuses
F = f(X,Y,Z);
my_integral = sum(F(:))*dx^3; %the answer is about 1.215

The problem with this numerical approach is that it runs into memory issues very quickly. For example, the storage requirement for what I wrote is proportional to dx^(-3) = 10^6. But if I had made dx=1e-3, I'd have needed 1000 times the storage which would exceed memory limitations. Also, if we want a quadruple integral rather than a triple integral, the storage scales like dx^-4 = 10^8 for dx = 1e-2. That is at the upper limit of what can be done on my machine.
Monte Carlo

f=@(x,y,z) (x.*y.*z).^(-x.*y.*z); %the function you want to integrate
n_trials = 1e7; %number of iid points to sample
r = rand(n_trials, 3); %iid samples uniform on unit cube
my_integral2 = mean( f( r(:,1), r(:,2), r(:,3) ) ); %the answer is about 1.215

This approach is nice because it scales nicely with dimension. If you integral had been 4 dimensional, we could use the same number of points and still get a good approximation.
Symbolic
This approach requires the symbolic toolbox and it still fails on this problem, but I'll write it anyway since it would work for problems that the symbolic integrator can understand.

syms x y z
my_integral3 = int(int(int((x*y*z)^(-x*y*z),x,0,1),y,0,1),z,0,1); %fails

Built-in Triplequad

f=@(x,y,z) (x.*y.*z).^(-x.*y.*z); %the function you want to integrate
tol = 1e-6; %tolerance
my_integral4 = triplequad(f, 0,1,0,1,0,1,tol) %returns 1.214832