3
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – m_goldberg Apr 18 at 4:47
4
$\begingroup$

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.

$\endgroup$
1
$\begingroup$

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
$\endgroup$
0
$\begingroup$

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.

$\endgroup$

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.