3
$\begingroup$

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.

$\endgroup$
2
$\begingroup$

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

$\endgroup$
  • 1
    $\begingroup$ Thanks Henrik. I may have misstated my problem. I would need to print the integration over x for a specific value of b. For example, let assume that 1<b<10 with incrementation of 0.2. I would need to print the integration over x for b = 1,1.2,1.4,.1.6....10. $\endgroup$ – Jean-Hubert Olivier Aug 8 at 14:05

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.