# Simple 3D trapzoidal rule gone wrong

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?

• The function is taken from here: utkstair.org/clausius/docs/che505/pdf/IE_eval_N-Dints.pdf page 4 Apr 23 at 8:34
• You have not stated which integral you are trying to compute. Apr 23 at 11:30
• @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}] Apr 23 at 11:39

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