2
$\begingroup$

I am trying to implement a 3D trapezoidal rule function, but I get wrong result and I do not understand what's wrong. Here is the code:

Int[xi_, yi_, zi_, xf_, yf_, zf_] := 
 Module[{gridr, gridth, gridph, int1, int2, int3, int4, int5, int6, 
   int7, int8, Int, dx = 0.2, dy = 0.2, dz = 0.2},
  gridr = Range[xi, xf, dx];
  gridth = Range[yi, yf, dy];
  gridph = Range[zi, zf, dz];
  int1 = f[xi, yi, zi] + f[xi, yi, zf] + f[xi, yf, zi] + 
    f[xi, yf, zf] + f[xf, yi, zi] + f[xf, yi, zf] + f[xf, yf, zi] + 
    f[xf, yf, zf];
  int2 = Apply[Plus, 
    Apply[Plus, f[xi, yi, #]] + Apply[Plus, f[xi, yf, #]] + 
       Apply[Plus, f[xf, yi, #]] + Apply[Plus, f[xf, yf, #]] & /@ 
     Rest[gridph]];
  int3 = Apply[Plus, 
    Apply[Plus, f[xi, #, zi]] + Apply[Plus, f[xi, #, zf]] + 
       Apply[Plus, f[xf, #, zi]] + Apply[Plus, f[xf, #, zf]] & /@ 
     Rest[gridth]];
  int4 = Apply[Plus, 
    Apply[Plus, f[#, yi, zi]] + Apply[Plus, f[#, yf, zi]] + 
       Apply[Plus, f[#, yi, zf]] + Apply[Plus, f[#, yf, zf]] & /@ 
     Rest[gridr]];
  int5 = Plus @@ 
     Flatten@Outer[f[xi, #1, #2] &, Rest[gridth], Rest[gridph]] + 
    Plus @@ Flatten@Outer[f[xf, #1, #2] &, Rest[gridth], Rest[gridph]];
  int6 = Plus @@ 
     Flatten@Outer[f[#1, yi, #2] &, Rest[gridr], Rest[gridph]] + 
    Plus @@ Flatten@
      Outer[f[#1, yf, #2] &, Rest[gridr], Rest[gridph]];
  int7 = Plus @@ 
     Flatten@Outer[f[#1, #2, zi] &, Rest[gridr], Rest[gridth]] + 
    Plus @@ Flatten@
      Outer[f[#1, #2, zf] &, Rest[gridr], Rest[gridth]];
  int8 = Plus @@ 
    Flatten@Outer[f[#1, #2, #3] &, Rest[gridr], Rest[gridth], 
      Rest[gridph]];
  Int = int1 + 2 (int2 + int3 + int4) + 4 (int5 + int6 + int7) + 
    8 int8;
  (dx dy dz )/8 Int]

The function f is given by:

f[x_, y_, z_] := x y z;

So the integral should be zero, but I get 4.096. What am I doing wrong?

Improve this question
$\endgroup$
3
  • $\begingroup$ The function is taken from here: utkstair.org/clausius/docs/che505/pdf/IE_eval_N-Dints.pdf page 4 $\endgroup$
    – mattiav27
    Apr 23 at 8:34
  • $\begingroup$ You have not stated which integral you are trying to compute. $\endgroup$
    – Henrik Schumacher
    Apr 23 at 11:30
  • $\begingroup$ @HenrikSchumacher I am trying to integrate the function f at the end of the post. The integral should give the same result of Integrate[f[x,y,z],{x,-8,8},{y,-8,8},{z,-8,8}] $\endgroup$
    – mattiav27
    Apr 23 at 11:39

1 Answer 1

Reset to default
5
$\begingroup$

I don't know what you are doing wrong, but the fact that you can hardly read it yourself should ring a bell: The more code you have, the more can go wrong.

Here is a shorter and more idiomatic implementation (hopefully correct).

ClearAll[Int];
Int[f, {xi_, xf_, xn_}, {yi_, yf_, yn_}, {zi_, zf_, zn_}] := 
  Module[{x, y, z, xw, yw, zw},
   x  = Subdivide[N[xi], N[xf], xn];
   y  = Subdivide[N[yi], N[yf], yn];
   z  = Subdivide[N[zi], N[zf], zn];
   xw = Normal[SparseArray[{{1} -> 0.5, {xn + 1} -> 0.5}, {xn + 1}, 1.]];
   yw = Normal[SparseArray[{{1} -> 0.5, {yn + 1} -> 0.5}, {yn + 1}, 1.]];
   zw = Normal[SparseArray[{{1} -> 0.5, {zn + 1} -> 0.5}, {zn + 1}, 1.]];
   
   (N[xf] - N[xi]) (N[yf] - N[yi]) (N[zf] - N[zi]) / N[xn yn zn] Dot[Outer[f, x, y, z], zw, yw, xw]
   ];

Int[f, {-8, 8, 90}, {-8, 8, 90}, {-8, 8, 90}]

2.45397*10^-11

Improve this answer
$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.