[Scilab-users] Int3D / Triple integration

[arctica1963]

Hello all,

An update on a solution. Following e-mail correspondence with a fellow
Scilab user (Javier Domingo), he has worked out a general solution to
vectorising the X,Y,Z arrays for the tetrahedrons required by int3d. So
taking this part of the code and adding a call to int3d within a function we
get a simpler route to doing triple integrals given a function (f) and lower
and upper limits of integration defined by x1,x2,y1,y2,z1,z2:

function [Integral, Error] = Integral_3d (f, x1, x2, y1, y2, z1, z2)
// Divide prism (given by abscissa: x1, x2, ordinate: y1 to y2 and z1 to z2)
into
// 12 tetrahedra (not regular), starting from the center of the prism, cover
all its volume;
// providing the array IX (abscissa of the vertices of the triangles),
// and the array IY (ordinate of the vertices of the triangles).
   xc = (x1 + x2) / 2; yc = (y1 + y2) / 2; zc = (z1 + z2) / 2; // center of
prism
// coordinates of the prism tips (2 prisms on each face)
// bottom -top- right -left- front -rear-
   LX = [xc, xc, xc, xc, xc, xc, xc, xc, xc, xc, xc, xc;
       x1, x1, x1, x1, x2, x2, x1, x1, x1, x1, x1, x1;
       x2, x1, x2, x1, x2, x2, x1, x1, x2, x1, x2, x1;
       x2, x2, x2, x2, x2, x2, x1, x1, x2, x2, x2, x2];
   LY = [yc, yc, yc, yc, yc, yc, yc, yc, yc, yc, yc, yc;
       y1, y1, y1, y1, y1, y2, y1, y2, y1, y1, y2, y2;
       y1, y2, y1, y2, y2, y2, y2, y2, y1, y1, y2, y2;
       y2, y2, y2, y2, y1, y1, y1, y1, y1, y1, y2, y2];
   LZ = [zc, zc, zc, zc, zc, zc, zc, zc, zc, zc, zc, zc;
       z1, z1, z2, z2, z1, z2, z1, z2, z1, z1, z1, z1;
       z1, z1, z2, z2, z1, z1, z1, z1, z1, z2, z1, z2;
       z1, z1, z2, z2, z2, z2, z2, z2, z2, z2, z2, z2];
       
[Integral, Error] = int3d (LX, LY, LZ, f, 1, [0,100000,1.d-5,1.d-7]);

endfunction

As a simple test one can define a function v=x^2 + y^2 + z^2 with limits of
0 to 1 - which simplifies to 1.0 as the sum of the iterated integrals.

deff('v=f(xyz,numfun)','v=xyz(1)^2+xyz(2)^2+xyz(3)^2')
x1=0;x2=1;y1=0;y2=1;z1=0;z2=1;

--> [Integral, Error] = Integral_3d (f, x1, x2, y1, y2, z1, z2)
 Integral  = 

   1.
 Error  = 

   1.110D-14

Thanks to Javier for his work on defining/clarifying the X,Y,Z arrays and
logic for defining the equation in Scilab. As a suggestion it would seem
reasonable to have this aspect either built into the function (int3d) or for
a separate mesh3d function to build tetrahedrons in a format compatible with
int3d.

Always great to exchange ideas to formulate a solution to a problem.

Code tested under Scilab version 6.1.0

Lester




--
Sent from: http://mailinglists.scilab.org/Scilab-users-Mailing-Lists-Archives-f2602246.html