3D integration over a region bounded by some planes.

I want to perform 3D integrations of some functions over a region defined by the following 14 planes. As an example, consider the function . The region is bounded by:

• 2 planes parallel to yz-plane at and .
• 2 planes parallel to xz-plane at and .
• 2 planes parallel to xy-plane at and .
• 8 planes which are penpendicular to the vectors - given in below code at mid point of these vector ().

Code to visulize the last 8 planes. First 6 planes are quite eays to imagine.

clear; clc;

figure; hold on;
axis([-1 1 -1 1 -1 1])

% function to calculate the 8 vectors:
f = @(vec) vec(1)*[-1 1 1] + vec(2)*[1 -1 1] + vec(3)*[1 1 -1];


% 8 vectors:
v1 = f([1, 1, 1]);
v2 = f([-1, -1, -1]);
v3 = f([0, 1, 0]);
v4 = f([0, -1, 0]);
v5 = f([0, 0, 1]);
v6 = f([0, 0, -1]);
v7 = f([1, 0, 0]);
v8 = f([-1, 0, 0]);


% draw planes perpendicular to v vectors at v/2 point:
draw_plane(v1)
draw_plane(v2)
draw_plane(v3)
draw_plane(v4)
draw_plane(v5)
draw_plane(v6)
draw_plane(v7)
draw_plane(v8)


xlabel('x');
ylabel('y');
zlabel('z');
box on;
set(gca,'fontname','times','fontsize',16)
view([-21   19])
hold off;

function [] = draw_plane(v)

% I took this algorithem from ChatGPT to draw a plane perpendicular to v

plane_size = 3;
midpoint = v./2;
normal_vector = v / norm(v);
perpendicular_vectors = null(normal_vector)';

[t1, t2] = meshgrid(linspace(-plane_size, plane_size, 100));
plane_points = midpoint + t1(:) * perpendicular_vectors(1,:) + t2(:) * perpendicular_vectors(2,:);
X = reshape(plane_points(:,1), size(t1));
Y = reshape(plane_points(:,2), size(t1));
Z = reshape(plane_points(:,3), size(t1));
surf(X, Y, Z, 'FaceAlpha', 0.5, 'EdgeColor', 'none');
plot3(midpoint(1), midpoint(2), midpoint(3), 'go', 'MarkerSize', 10, 'MarkerFaceColor', 'g');
text(midpoint(1), midpoint(2), midpoint(3), ['(',num2str(v(1)),',',num2str(v(2)),',',num2str(v(3)),')']);

end

[Torsten]

Either use a Delaugnay triangulation of the region of interest and sum the integrals over the triangles or use Monte- Carlo Integration.

[Luqman Saleem]

@Torsten Thank you so much for your comments. You are always here for me whenever I come here with an integration problem.

Could you please provide me with a simple example of any of these integration methods for a 3D integral? It would be a great help.

[Torsten]

Can you provide the region you want to integrate over as a system of linear inequalities ?

Or do you have another method to decide whether a point in [-1,1]x[-1,1]x[-1,1] lies in the region you want to integrate over ?

E.g. [-1,1]x[-1,1]x[-1,1] could be given as

x>=-1
y>=-1
z>=-1
x<=1
y<=1
z<=1

[Luqman Saleem]

Yes I think the region bounded by these planes is:

equ1 = x <= 1;
equ2 = x >= -1;
equ3 = y <= 1;
equ4 = y >= -1;
equ5 = z <= 1;
equ6 = z >= -1;


equ7  = x + y + z <= 3/2; 
equ8  = x + y + z >= -3/2; 
equ9  = x - y + z <= 3/2;
equ10 = x - y + z >= -3/2;
equ11 = x + y - z <= 3/2;
equ12 = x + y - z >= -3/2;
equ13 = -x + y + z <= 3/2;
equ14 = -x + y + z >= -3/2;

Answer 1

[Torsten]

f = @(x,y,z)x.^2+y+2*z;
N = [100,1000,10000,100000,1000000,10000000,100000000];
Fi = zeros(numel(N),1);
Fo = Fi;
for i = 1:numel(N)
  n = N(i);  
  x = -1+2*rand(n,1);
  y = -1+2*rand(n,1);
  z = -1+2*rand(n,1);
  idx = x+y+z<=1.5 & x+y+z>=-1.5 & x-y+z<=1.5 & x-y+z>=-1.5 & ...
        x+y-z<=1.5 & x+y-z>=-1.5 & -x+y+z<=1.5 & -x+y+z>=-1.5;
  Fi(i) = 8*sum(f(x(idx),y(idx),z(idx)))/n;  % 8 is the volume of [-1,1]x[-1,1]x[-1,1]
  idx = ~idx;
  Fo(i) = 8*sum(f(x(idx),y(idx),z(idx)))/n;  % 8 is the volume of [-1,1]x[-1,1]x[-1,1]
end
Fi+Fo

ans = 7x1

0.1365 2.9963 2.7342 2.7236 2.6703 2.6643 2.6667

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

plot(1:numel(N),[Fi,Fo,Fi+Fo])
grid on

integral3(f,-1,1,-1,1,-1,1)

ans = 2.6667

f = @(x,y,z)(x.^2+y+2*z).*...
            (x+y+z<=1.5).*(x+y+z>=-1.5).*(x-y+z<=1.5).*(x-y+z>=-1.5).*...
            (x+y-z<=1.5).*(x+y-z>=-1.5).*(-x+y+z<=1.5).*(-x+y+z>=-1.5);
integral3(f,-1,1,-1,1,-1,1)

Warning: Reached the maximum number of function evaluations (10000). The result fails the global error test.

Warning: The integration was unsuccessful.

ans = NaN

For the details, take a look at

https://en.wikipedia.org/wiki/Monte_Carlo_integration

[Luqman Saleem]

worked like a charm. thank you very much.