3D integration over a region bounded by some planes.

3 次查看（过去 30 天）
Luqman Saleem
Luqman Saleem 2024-8-8
编辑： Torsten 2024-8-8
Ran in:
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
  4 个评论
显示 2更早的评论
Torsten
Torsten 2024-8-8
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
Luqman Saleem 2024-8-8
编辑：Luqman Saleem 2024-8-8
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;

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2024-8-8
编辑：Torsten 2024-8-8
Ran in:
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

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Logical 的更多信息

标签

产品


版本

R2024a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by