# Numerical Integration Via Adaptive Simpson's Method

I have written a code that uses the Adaptive Simpson's method to approximate integration. For those who are unaware of this Adaptive Simpson's method; Adaptive Simpson's method

In my code, I count the number of function evaluations are needed. I am wondering if there is a way to reduce the number of function evaluations needed but still retaining the same outputs and error. Here is my code:

Simpson[f_, a_, b_, er_] := Module[{s1, s2, h2, h = (b - a)/2},

s1 = h (f[a] + 4 f[a + h] + f[b])/3;
h2 = h/2;
s2 = h2 (f[a] + 4 f[a + h2] + 2 f[a + h] + 4 f[b - h2] + f[b])/3;

If[Abs[s1 - s2]/15 < er, s2 + (s2 - s1)/15,
Simpson[f, a, a + h, er/2] + Simpson[f, a + h, b, er/2]]];
count = 0;
h[x_] := Module[{}, count++; 3 Sqrt[x]];
N[Simpson[h, 0, 1, 10.^-5] - 2]
count


For example, when approximating the integral of $3*sqrt(x)$ from 0 to 1, my algorithm makes 312 function calls of the function h, but I suspect some of these evaluations are repeated.

With a simple modification to your code

Simpson[f_, a_, b_, er_] :=
Module[{s1, s2, h2, h = (b - a)/2},
s1 = h (f[a] + 4 f[a + h] + f[b])/3;
h2 = h/2;
s2 = h2 (f[a] + 4 f[a + h2] + 2 f[a + h] + 4 f[b - h2] + f[b])/3;
If[Abs[s1 - s2]/15 < er, s2 + (s2 - s1)/15,
Simpson[f, a, a + h, er/2] + Simpson[f, a + h, b, er/2]]];
count = 0;
h[x_] := Module[{}, Sow[x]; count++; 3 Sqrt[x]];
{result, xvalues} = Reap[N[Simpson[h, 0, 1, 10.^-5] - 2]];
{Length[Union[Flatten[xvalues]]], count}
(* {81, 312} *)


we can see that indeed your code evaluates h 312 times but at only 81 distinct values.

If the execution time of the integration is dominated by evaluating h this is clearly undesirable. A simple fix is to use "memoisation". This avoids re-evaluating h for values of x where it has already been computed. This can easily be done by using a definition for h such as

h[x_]:=h[x]= 3 Sqrt[x]


There is a very clear discussion of this in the Mathematica documentation https://reference.wolfram.com/language/tutorial/FunctionsThatRememberValuesTheyHaveFound.html

• Thank you for responding. I am a little confused to what you did with the last 3 lines of the code. I was thinking more of along the lines of altering the algorithm so repeated divisions of 3 aren't done so less function calls are needed. – larry Oct 9 '16 at 19:33
• I am just not sure how to alter the algorithm to remove the repeated divisions of 3 without changing the answer/error. – larry Oct 9 '16 at 19:34
• I think my answer is a simple solution to your problem - I've updated it and hope it is now clearer. – mikado Oct 9 '16 at 20:05

An alternative to memoization is to pass along the function values already obtained to the Simpson's routine. This data will be stored in the stack until no longer needed, since the algorithm is recursive. It should not be too bad in practice though, unless you set a very small error.

I would also suggest converting exact input to approximate real numbers. Otherwise, the expression could possibly grow to be an extremely large symbolic representation of an exact number. This could be slow, eat up lots of memory and even crash Mathematica (it did for me, when I pushed it too far). Uncomment the (*N@*) to do this.

I kept the OP's Simpson[] user-interface and used an "internal" version iSimpson for the computation. The function iSimpson[] has for arguments, the function f, the sample points {a, c, b} with c == (a+b)/2, the corresponding function values {fa, fc, fb} and the error tolerance er. The function values are reused.

Simpson[f_, a_, b_, er_] :=
iSimpson[f, (*N@*){a, (a + b)/2, b}, (*N@*){f[a], f[(a + b)/2], f[b]}, er];
iSimpson[f_, {a_, c_, b_}, {fa_, fc_, fb_}, er_] :=
Module[{c1, c2, f1, f2, s1, s2, h2},
h2 = (b - a)/4;                             (* next half step size *)
c1 = a + h2;                                (* next midpoints *)
c2 = b - h2;
f1 = f[c1];                                 (* function values at next midpoints *)
f2 = f[c2];
s1 = 2 h2 (fa + 4 fc + fb)/3;               (* integrals *)
s2 = h2 (fa + 4 f1 + 2 fc + 4 f2 + fb)/3;

If[Abs[s1 - s2]/15 < er,
s2 + (s2 - s1)/15,
iSimpson[f, {a, c1, c}, {fa, f1, fc}, er/2] +
iSimpson[f, {c, c2, b}, {fc, f2, fb}, er/2]]
];

count = 0;
h[x_] := Module[{}, count++; 3 Sqrt[x]];
N[Simpson[h, 0, 1, 10.^-5] - 2]
count
(*
-5.73618*10^-7
81
*)


You can also automatically memoize (@mikado's solution) as follows. This way the values are only temporarily stored during the computation of the integral. Again I use an outer function Simpson2 that defines a memoized ff and calls the inner iSimpson2, which has the same code as the OP's Simpson.

ClearAll[Simpson2, iSimpson2];
Simpson2[f_, a_, b_, er_] := Module[{ff},
ff[x_] := ff[x] = f[x];
iSimpson2[ff, a, b, er]
];
iSimpson2[f_, a_, b_, er_] :=
Module[{s1, s2, h2, h = (b - a)/2},
s1 = h (f[a] + 4 f[a + h] + f[b])/3;
h2 = h/2;
s2 = h2 (f[a] + 4 f[a + h2] + 2 f[a + h] + 4 f[b - h2] + f[b])/3;
If[Abs[s1 - s2]/15 < er, s2 + (s2 - s1)/15,
Simpson2[f, a, a + h, er/2] + Simpson2[f, a + h, b, er/2]]];

count = 0;
h[x_] := Module[{}, count++; 3 Sqrt[x]];
N[Simpson2[h, 0, 1, 10.^-5] - 2]
count
(*
-5.73618*10^-7
81
*)


Update: Performance comparison

The difference in performance is greater than I imagined and so is worth reporting.

Timing

h[x_] := Module[{}, count++; Sin[x^2]];
count = 0;
{Simpson[h, 0., 12., 10.^-12], count} // AbsoluteTiming
count = 0;
{Simpson2[h, 0., 12., 10.^-12], count} // AbsoluteTiming
(*
{1.75486, {0.590432, 90077}}
{7.98656, {0.590432, 90077}}
*)


Memory use

MemoryConstrained[Simpson[h, 0., 12., 10.^-12], 100000]
MemoryConstrained[Simpson[h, 0., 12., 10.^-12], 110000]
(*
$Aborted 0.590432 *) MemoryConstrained[Simpson2[h, 0., 12., 10.^-12], 35000000] MemoryConstrained[Simpson2[h, 0., 12., 10.^-12], 36000000] (*$Aborted
0.590432
*)


The difference in memory use is predictable, although I hadn't investigated just how much of a difference it could be. Evaluation in the recursive Simpson[] makes a depth-first traversal of the evaluation tree, and when evaluation of a sub-branch finishes, the memory storing the sample points and function values is released. In the memoized version, all of the downvalues of ff are accumulated until the entire computation is performed. This overhead evidently adds considerably to the evaluation time as well.