Evaluating the error function by integrating $\mathrm e^{−t^2}$ with Simpson's rule

I am trying to evaluate an error function with Simpson's rule because there is no other way to integrate it. The function is

$${\rm erf}(x) = \frac{2}{\sqrt π}\int_0^x \mathrm e^{−t^2}\, \mathrm dt$$

and I want the value at 1.

I typed in the code:

Simpson[a_, b_, m_] :=
Module[{},
h = (b - a)/(2 m);
sum = 0;
s1 = 0;
For[k = 1, k <= m - 1, k++, s1 = s1 + f[a + h (2 k - 1)];];]
s2 = 0
For[k = 1, k <= m, k++, s2 = s2 + f[a + h (2 k - 1)];];
sum = h/3 (f[a] + f[b] + 2 s1 + 4 s2);
Return[sum]


and, of course, an error pops up in my output.

Does anyone know how to fix the error part so I can get a numerical answer?

• There are some serious syntax errors in your code. Also, are you aware that Erf is a built-in function and that Erf[1] // N returns 0.842701? – m_goldberg Apr 18 at 4:47

When I rewrite your code as

simpson[a_, b_, m_] :=
Module[{f, h, s1, s2, k},
f[t_] := E^-(t^2);
h = (b - a)/(2 m);
s1 = 0;
For[k = 1, k <= m - 1, k++, s1 = s1 + f[a + h (2 k - 1)]];
s2 = 0;
For[k = 1, k <= m, k++, s2 = s2 + f[a + h (2 k - 1)]];
h/3 (f[a] + f[b] + 2 s1 + 4 s2)]


then evaluating

2 N[simpson[0, 1, 500]]/Sqrt[Pi]


gives

0.842938

which is not a bad approximation. As to there being no better way than Simpson's rule, the built-in function NIntegrate, is certainly a better way.

Block[{f}, f[t_] := E^-(t^2); 2 NIntegrate[f[t], {t, 0, 1}]/Sqrt[Pi]]


gives the much better result

0.842701

Note that the built-in function when evaluated numerically

Erf[1] // N


gives

0.842701

which matches the NIntegrate result up to display precision.

If you want it, you can also use the following version:

n = 8;
a = 0;
b = 1;
h = (b - a)/n;
f[t_] := E^-(t^2);

simpson =
h/6 Sum[f[a + i h] + 4 f[a + (i + 1/2) h] + f[a + (i + 1) h], {i, 0, n - 1}] // N;

2*simpson/Sqrt[\[Pi]]
0.842701


The major point is the list of coefficients {1, 4, 2, 4, 2, ..., 1}, which I use a SparseArray to realize below. And the final sum can be conceived as the inner product of two vectors, realized by Dot:

Clear[nIntegrateBySimpson, erfintegrand]

nIntegrateBySimpson[func : (_Function | _Symbol), xrange_, n_Integer] :=
Module[{simpsoncoefficients, h, subd},
h = -Subtract @@ xrange/n;
subd = Subdivide[##, n] & @@ N[xrange];
simpsoncoefficients = SparseArray[{{1} -> 1., {i_?EvenQ} -> 4., {-1} -> 1.},
n + 1, 2.
];
h/3 simpsoncoefficients.(func /@ subd)
]

erfintegrand[t_] := 2/Sqrt[\[Pi]] Exp[-t^2]


So, the numerical integral can be obtained by following codes:

nIntegrateBySimpson[erfintegrand, {0, 1}, 12]

0.842701


One notices that the range on which the integral is performed need not be divided into too many parts, as shown above, only 12 samples have been collected.