What value of $N$ to use in Simpson's rule to reach desired accuracy?

Question

I need to calculate Simpsons rule for the integral of $$\frac{e^x-1}{\sin x}$$ from $0$ to $\pi/2$ with minimum number of intervals $N$ up to $10^{-6}$ accuracy. Wolfram alpha seems to be giving me a different answer than what I'm getting. In short, what should the value of $N$ be?

public static double simpsonsRuleFunction2(double valueN, double valueA, double valueB, double valueDx) {

    double e = 0.0;
    double simpsonsRule = 0.0;
    double valueHolder = 0.0;
    valueN = 2;
    valueA = 0;
    valueB = (Math.PI)/2; 

    for(int i = 1; i<=valueN + 1; i++) {
        valueDx = (valueB-valueA)/valueN;
        e = valueA + ((i-1)*valueDx);

        if (i==1) {
        // Limit as x -> 0
            simpsonsRule += 0;

        } else if (i % 2 == 0) {
            simpsonsRule += 4*((Math.exp(e)-1)/Math.sin(e));
        } else if ((i % 2 != 0) && ( i > 1) && (i < valueN + 1)) {
            simpsonsRule += 2*((Math.exp(e)-1)/Math.sin(e));   
        } else if (i == valueN + 1) {
            simpsonsRule += ((Math.exp(e)-1)/Math.sin(e));
        }
        System.out.println("e: " + e);
        System.out.println("simpsonsRule2: " + simpsonsRule);

    }
    simpsonsRule = simpsonsRule *((valueDx)/3);

    System.out.println("simpsonsRule2: " + simpsonsRule);

    while(Math.abs(valueHolder - simpsonsRule) > Math.pow(10,-6)) {
        System.out.println("\nValueHolder2" + valueHolder);
        valueHolder = simpsonsRule;
        valueN +=2;
        valueDx = (valueB-valueA)/valueN;
        simpsonsRule = 0;
        for(int i = 1; i<=valueN + 1; i++) {
            e = valueA + ((i-1)*valueDx);

            if (i==1) {
            // Limit as x -> 0
                simpsonsRule += 1;
            } else if (i % 2 == 0) {
                simpsonsRule += 4*((Math.exp(e)-1)/Math.sin(e));
            } else if ((i % 2 != 0) && ( i > 1) && (i < valueN + 1)) {
                simpsonsRule += 2*((Math.exp(e)-1)/Math.sin(e));   
            } else if (i == valueN + 1) {
                simpsonsRule += ((Math.exp(e)-1)/Math.sin(e));
            }
            System.out.println("e: " + e);
            System.out.println("simpsonsRule2: " + simpsonsRule);
            System.out.println("valueB " + valueB);
            System.out.println("valueDx" + valueDx);
        }
        simpsonsRule = simpsonsRule *((valueDx)/3);
        System.out.println("simpsonsRule2: " + simpsonsRule);
    }
    return valueN;

Answer 1

If you use Simpson's Rule with $N$ subintervals to estimate an integral, there is a well-know upper bound on the error of your estimation: $$E\leq\frac{1}{180}\frac{(b-a)^5}{N^4}M$$ where $M$ is the maximum absolute value obtained by the fourth derivative of $f$. You can always improve upon this by calculating Simspon's Rule over two or more subintervals and adding the results, so that for some of the subintervals, you could theoretically have smaller values of $M$. But putting that aside, you need: $$\frac{1}{180}\frac{(\pi/2-0)^5}{N^4}M\leq10^{-6}$$ $$N\geq\left[\frac{10^6\pi^5}{180(32)}M\right]^{1/4}$$

Now what is $M$? Well, you have to find the fourth derivative of $f$ and then maximize that, which may involve calculating the 5th derivative of $f$. Of course, you can also allow yourself a weaker bound (and therefore ask yourself to use a larger $N$) and just find some value for $M$ that you can establish is definitely larger than the actual maximum absolute value of $f^{(4)}$. But start by seeing if you can find the actual value for $M$.

I didn't investigate formally, but visually it appears that $f^{(4)}$ here is an increasing function on $[0,\pi/2]$. (Maybe you can prove that by looking at $f^{(5)}$.) So then it's easy to get a value for $M$ by just taking $M=f^{(4)}(\pi/2)$.

Lstly I'd just stress that this provides you with an $N$ that is guaranteed to give the desired accuracy. It should be noted that possibly, a smaller $N$ would also give the desired accuracy. Basically, the $M$ we would use is only doing a good job for the estimate on the far right of the interval. If the whole thing were cut up into smaller subintervals to perform Simpson on, you could use smaller $M$ for all but the last subinterval.

Answer 2

Pragmatic view: Use Richardson extrapolation to get error estimates and estimates for the constants.

The error is by theory $E_N=S_N-S_*=C·N^{-4}+O(N^{-5})$. Use $$ S_N-S_{2N}=\frac{15}{16}·C·N^{-4}+O(N^{-5}) $$ to get an estimate for $E_N\approx \frac{16}{15}(S_N-S_{2N})$. Then the new step number should be around $N\cdot\sqrt[4]{|E_N|·10^6}$.

Start with $N=8$ (intervals each with its midpoint, so $2N+1=17$ nodes), this should already give an error in the right magnitude.

---

General remarks on the code:

• Use more structure in the code. Use separate functions for function evaluation and Simpson rule.
• Investigate the use of expm1 for a better evaluation of $e^x-1$ for small values, if such are used in the sampling nodes.

---

Carrying out these ideas gives

a,b = 1e-20, pi/2 # vectorized evaluation does not play nice with branches for singularities
f=lambda x: expm1(x)/sin(x)
def test(N): # N = number of segments with midpoint
    x = np.linspace(a,b, 4*N+1);
    fx = f(x);
    SN = (0.5*(fx[0]+fx[-1]) + 2*sum(fx[2::4]) + sum(fx[4:-4:4]))*(b-a)/(3*N);
    S2N = (0.5*(fx[0]+fx[-1]) + 2*sum(fx[1::2]) + sum(fx[2:-2:2]))*(b-a)/(6*N);

    EN = abs(SN-S2N)*16/15
    print("used N=",N)
    print("S_N  = %.16f, error estimate E_2N = %.6g\nS_2N = %.16f"%(SN, EN, S2N))
    print("proposed N=",N*(EN*1e+6)**0.25)
test(8);

returning

used N= 8
S_N  = 3.0016050095733235, error estimate E_2N = 9.43905e-06
S_2N = 3.0015961604675852
proposed N= 14.0223891023

so that using $N=15$ should give the desired accuracy.

used N= 15
S_N  = 3.0015963354073461, error estimate E_2N = 7.69719e-07
S_2N = 3.0015956137958382
proposed N= 14.0499345967