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.