Implementing Simpson's rule with two variables

Question

For one of my lab project I need to integrate a given function (in this example, f1, using Simpson's ruls. When only one variable x exists, I was able to find a code that calculate the integration:

f1[x_] := 3 + Cos[x + 3]
LEFT1[f_, a_, b_, n_] := 
  Module[{h, i},
    h = (b - a)/n;
    N[h Sum[f1[a + i h], {i, 1, n - 1, 2}]]]
RIGHT1[f_, a_, b_, n_] := 
  Module[{h, i},
    h = (b - a)/n;
    N[h Sum[f[a + i h], {i, 2, n - 2, 2}]]]
TRAP1[f_, a_, b_, n_] := 4 LEFT1[f, a, b, n] + 2 RIGHT1[f, a, b, n]
MID1[f_, a_, b_, n_] := 
  Module[{h, i},
    h = (b - a)/n; 
    N[h (f[a] + f[b])]]
SIMP1[f_, a_, b_, n_] := (MID1[f, a, b, n] + TRAP1[f, a, b, n])/3.0

SIMP1[f1, -1, 10, 10]

32.5062

Instead of using f1[x_] := 3 + Cos[x + 3], I would like to use f1[x_] = 3 + Cos[x + 3 b] and calculate the integration for value of b ranging from 1 to 50. I am not sure how I should start tackling this problem.

Answer 1

How about this? It implements the tensor product rule associated to the quadrature weights of Simpson 1D-rule.

Simpson2D[f_, {xa_, xb_, xn_}, {ya_, yb_, yn_}] := 
 Module[{Δx, xω, Δy, yω},
  Δx = (xb - xa)/(2. xn);
  xω = ConstantArray[2./3. Δx, 2 xn + 1];
  xω[[2 ;; ;; 2]] = 4./3. Δx;
  xω[[1]] = 1./3. Δx;
  xω[[-1]] = 1./3. Δx;

  Δy = (yb - ya)/(2. yn);
  yω = ConstantArray[2./3. Δy, 2 yn + 1];
  yω[[2 ;; ;; 2]] = 4./3. Δy;
  yω[[1]] = 1./3. Δy;
  yω[[-1]] = 1./3. Δy;
  xω.Outer[f, Subdivide[N[xa], N[xb], 2 xn], Subdivide[N[ya], N[yb], 2 yn]].yω
  ]

Test:

f = {x, y} \[Function] 3. + Cos[x + y];
a = Simpson2D[f, {-1, 10, 100}, {1, 50, 100}];
b = Integrate[f[x, y], {x, -1, 10}, {y, 1, 50}];
Abs[1 - a/b]

3.21721*10^-9