[Scilab-users] Int3D / Triple integration

classic Classic list List threaded Threaded
8 messages Options
arctica1963 arctica1963
Reply | Threaded
Open this post in threaded view
|

[Scilab-users] Int3D / Triple integration

Hello,

Following on from a previous query I had on this matter, I am not having
much success with int3d for triple integration, largely due to the need to
build tetrahedrons. CGLab does have this option (delauney3d) but is not
available for 6.1.0, and it is not entirely clear how it would interface
with int3d.

Monte-Carlo integration is one route - located this online (added xmin xmax
etc):

function mcInt3d = intg3d(f,xmin,xmax,ymin,ymax,zmin,zmax,N)
vol = (xmax-xmin)*(ymax-ymin)*(zmax-zmin)        
xr = (xmax-xmin)*rand(1,N)+xmin        
yr = (ymax-ymin)*rand(1,N)+ymin        
zr = (zmax-zmin)*rand(1,N)+zmin        
sumf= sum(f(xr,yr,zr))  
mcInt3d = vol*sumf/N
endfunction

It does work but in order to get close to a good approximation, N needs to
be quite large. I am sure if we could get int3d to be more user friendly
that would be great.

More of a shout out to anyone who has used int3d using limits defined for
x,y,z and how to vectorise the tetrahedrons. Would love to be able to crack
this problem.

Thanks

Lester



--
Sent from: http://mailinglists.scilab.org/Scilab-users-Mailing-Lists-Archives-f2602246.html
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
arctica1963 arctica1963
Reply | Threaded
Open this post in threaded view
|

Re: Int3D / Triple integration

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
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
mottelet mottelet
Reply | Threaded
Open this post in threaded view
|

Re: Int3D / Triple integration

Hi Lester,

If I understand well, you are only interested by integrating on

[x1,x2] x [y1,y2] x [z1,z2]

and not a general volume, that's it ?

S.

Le 02/04/2021 à 11:53, arctica1963 a écrit :

> 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: https://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/mailinglists.scilab.org/Scilab-users-Mailing-Lists-Archives-f2602246.html
> _______________________________________________
> users mailing list
> [hidden email]
> https://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/lists.scilab.org/mailman/listinfo/users

--
Stéphane Mottelet
Ingénieur de recherche
EA 4297 Transformations Intégrées de la Matière Renouvelable
Département Génie des Procédés Industriels
Sorbonne Universités - Université de Technologie de Compiègne
CS 60319, 60203 Compiègne cedex
Tel : +33(0)344234688
http://www.utc.fr/~mottelet

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
arctica1963 arctica1963
Reply | Threaded
Open this post in threaded view
|

Re: Int3D / Triple integration

Hi Stephane,

At the moment I am just trying to understand how Scilab works with triple
integration of f(x,y,z) with limits for xyz.

Lester



--
Sent from: http://mailinglists.scilab.org/Scilab-users-Mailing-Lists-Archives-f2602246.html
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
mottelet mottelet
Reply | Threaded
Open this post in threaded view
|

Re: Int3D / Triple integration


Le 02/04/2021 à 12:53, arctica1963 a écrit :
Hi Stephane,

At the moment I am just trying to understand how Scilab works with triple
integration of f(x,y,z) with limits for xyz. 

Ok, when you say "limits" for xyz you mean that each variable varies in a given constant interval, that's what I meant by the rectangular parallelepiped [x1,x2] x [y1,y2] x [z1,z2]. In fact it is a pity that Scilab does not handle this case but only the more general case of a collection of (eventually disconnected) tetrahedrons. However, cutting your parallepiped in 5  (https://www.geogebra.org/m/C3TjXxFY)  is enough to use int3d, since they will be recursively divided to attain the required precision:

deff('v=f(xyz,numfun)','v=xyz(1)^2+xyz(2)^2+xyz(3)^2')
xlim=[0 1];
ylim=[0 1];
zlim=[0 1];
[x,y,z]=ndgrid(xlim,ylim,zlim);
i = [5 8 2 3
     5 8 2 6
     5 8 3 7
     5 2 3 1
     2 3 8 4]';
[result,err] = int3d(x(i),y(i),z(i),f)


--> result
 result  = 

   1.0000000

--> err
 err  = 

   1.110D-14

S.

Lester



--
Sent from: https://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/mailinglists.scilab.org/Scilab-users-Mailing-Lists-Archives-f2602246.html
_______________________________________________
users mailing list
[hidden email]
https://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/lists.scilab.org/mailman/listinfo/users
-- 
Stéphane Mottelet
Ingénieur de recherche
EA 4297 Transformations Intégrées de la Matière Renouvelable
Département Génie des Procédés Industriels
Sorbonne Universités - Université de Technologie de Compiègne
CS 60319, 60203 Compiègne cedex
Tel : +33(0)344234688
http://www.utc.fr/~mottelet

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
arctica1963 arctica1963
Reply | Threaded
Open this post in threaded view
|

Re: Int3D / Triple integration

Hi Stephane,

Thanks for the information and methodology, useful to know. Learn something
new all the time with Scilab!

Kind regards

Lester



--
Sent from: http://mailinglists.scilab.org/Scilab-users-Mailing-Lists-Archives-f2602246.html
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
arctica1963 arctica1963
Reply | Threaded
Open this post in threaded view
|

Re: Int3D / Triple integration

Quick query re: your code,

How is the index (i) defined?
i = [5 8 2 3
     5 8 2 6
     5 8 3 7
     5 2 3 1
     2 3 8 4]';

Just trying to fully understand your method.

Lester



--
Sent from: http://mailinglists.scilab.org/Scilab-users-Mailing-Lists-Archives-f2602246.html
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
mottelet mottelet
Reply | Threaded
Open this post in threaded view
|

Re: Int3D / Triple integration

They are the number of the vertices of the cube, as generated by ndgrid:

--> [(1:8)' x(:) y(:) z(:)]
  ans  =

    1.   0.   0.   0.
    2.   1.   0.   0.
    3.   0.   1.   0.
    4.   1.   1.   0.
    5.   0.   0.   1.
    6.   1.   0.   1.
    7.   0.   1.   1.
    8.   1.   1.   1.

S.

Le 02/04/2021 à 15:49, arctica1963 a écrit :

> Quick query re: your code,
>
> How is the index (i) defined?
> i = [5 8 2 3
>       5 8 2 6
>       5 8 3 7
>       5 2 3 1
>       2 3 8 4]';
>
> Just trying to fully understand your method.
>
> Lester
>
>
>
> --
> Sent from: https://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/mailinglists.scilab.org/Scilab-users-Mailing-Lists-Archives-f2602246.html
> _______________________________________________
> users mailing list
> [hidden email]
> https://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/lists.scilab.org/mailman/listinfo/users

--
Stéphane Mottelet
Ingénieur de recherche
EA 4297 Transformations Intégrées de la Matière Renouvelable
Département Génie des Procédés Industriels
Sorbonne Universités - Université de Technologie de Compiègne
CS 60319, 60203 Compiègne cedex
Tel : +33(0)344234688
http://www.utc.fr/~mottelet

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users