Approximate This Integral!

The midpoint rule is simple to code and has a reasonable error bound. First, we code the midpoint rule in Python:

def midpoint_integral(f, a, b, n):
x = a+(b-a)/(2*n)
dx = (b-a)/n
ans = 0
while x < b:
    ans += f(x)*dx
    x += dx
return ans

We then solve for the minimum n required to achieve error less than 0.0005. The error of the midpoint rule is known to be $k \frac{(b-a)^{3}}{24n^{2}}$ , where k is the maximum value over the interval [a, b] of $|f^{``}(x)|$ . In this case, k = 100, as can be determined through basic calculus. We then solve the inequality: $\frac{100(\frac{\pi}{4})^{3}}{24n^{2}} < 0.0005$ and arrive at n being at least 45. Running the command

int(1000*midpoint_integral(lambda x: math.cos(x) ** 100, 0, math.pi/4, 45))

gives a result of $\textbf{125}$ .

$\displaystyle \int\limits_{0}^{\pi/4} \left(\cos x\right)^{100} dx\\ \displaystyle\approx \int\limits_{0}^{\pi/4} \left(1-\frac{x^2}{2}\right)^{100} dx \\ \displaystyle\approx \int\limits_{0}^{\pi/4} e^{-50x^2} dx \\ \displaystyle\approx \int\limits_{0}^{\infty} e^{-50x^2} dx \\ \displaystyle= \frac12 \sqrt{\frac{\pi}{50}} \approx 0.125331\ldots$

We approximate the function for rectangles in the interval between the limits of integration. Since the integral is the area below the curve of the function, we gotta add the areas of the rectangles.

Code in C:

#include <stdio.h>
#include <math.h>

#define PI acos(-1)

int main(void) {
    double interval = PI/400000; /* small sub-interval */
    double next;
    double integral = 0;
    int i;
    double cosseno;

    for (next = interval; next <= PI/4.0; next += interval) {
        cosseno = cos(next);
        for (i = 0; i < 99; i++) { /* really inefficient multiplication*/
            cosseno *= cos(next);
        }

        integral += cosseno * interval;
    }

    printf("integral = %lf\n", integral);

    return 0;
}

With a small enough $\epsilon$ , simply calculating the sum of $\epsilon\cdot\cos^{100}(i\epsilon)$ with $i$ from 0 to $\dfrac{\pi}{4\epsilon}$ will give the answer by the definition of integration.

Well, this is a CS problem, so I solved it the CS way using Python. The program is quite simple.

for i in range(1,100): # Just iterate the problem 99 times

 A = 0                         # A is the area under the curve y=cos(x)**100 which is the integral N

 n = 200*i

 d = pi/4/n                # Setting the slice dx

 for j in range(n):     # Approximation of the integration

      x = d*i

      y = cos(x)**100

      A+= y*d

 print n, A

After the computation finish just inspect the results for the first three digits after the decimal, which was found to be $\boxed{125}$ .

Python code:

import math

l=.0000001

x=0

b=float(math.pi)/4

s=0.0 while x<'b: y=l (math.cos(x)) *100

    x=x+l

    s=s+y

print s

I used Riemann sum rectangle approximations to estimate the value of the integral.

def integralValue():
    """Approximates the value of an integral"""
    from math import cos, pi

    # We use the classic Riemann sum rectangle approximation method (left-hand intervals)
    # Using 1000 rectangles should give enough accuracy

    integralValue = 0

    for x in range(1000):
        x *= pi/4000
        # Multiplies x by the width of the rectangle to obtain the correct x-value

        rectangleArea = (cos(x) ** 100) * (pi / 4000)

        # The width of each rectangle is pi/4000 and the height is cos^100(x)

        integralValue += rectangleArea

    return integralValue

Since only three digits of this are needed, we can approximate, rather than solve, this integral. To do this, we can break the function into some number of small pieces (~10000 is easily computable and gives enough precision).

The area of each piece can be computed by calculating the value of the function at that point (the height of the rectangle) and multiplying by 1/10000th of the width of the integration domain. By approximating the Riemann sum, we can approximate a difficult integral quite simply.

Let f(x) = $cos(x)^{100}$ and F(x) integral of f(x). Than N = F(pi/4) - F(0).

Let F(0) = c

Let e = $10^{-4}$

Knowing that f(x) = $\frac{F(x + e) - F(x)}{e}$ => F(x + e) = F(x) + e*f(x) we can write the series:

• F(0) = c
• F(e) = c + f(0)*e
• F(2e) = F(e) + f(e) e = c + f(0) e + f(e)*e

....

• F(pi/4) = c + e*(f0 + f(e) + f(2e) +...+f(pi/4 - e))

Than N = F(pi/4) - F(0) = e (f(0) + f(e) + f(2e) +...+ f(pi/4 - e)). Here is the code:

    double N = 0;
    double e = 0.0001;

    for(double x = 0; x < Math.PI / 4; x += e) {
        N += f(x) * e;
    }

    System.out.println(N);

Simpson's Rule

from math import *

def x(i):
    global a, h
    return a+h*i

def f(x):
    return cos(x)**100

a = 0
b = pi/4
n = 10000
h = (b-a)/(2*n)

res = f(x(0)) + f(x(2*n))
for i in range(1, 2*n):
    res += ((i%2 == 1)*2 + 2)*f(x(i))
res *= h/3

print res

Despite that this problem is solvable in literally three seconds in Mathematica, I wrote up a neat little numerical integrator that utilizes the trapezoid method. I am too lazy to comment it out, but the variable names are pretty self-explanatory. It seems to work on 100000 ordinates just fine.

import math

def trapezoid_integrator(lower_limit, upper_limit, ordinates):
    c = float(upper_limit - lower_limit) / ordinates

    integrand_values = []

    for i in range(ordinates + 1):
        integrand_val = (math.cos(i))**100
        integrand_values.append(integrand_val)

    trapezoid_sum = 0

    for i in range(len(integrand_values)):
        if (i == 0) or (i == len(integrand_values)-1):
            trapezoid_sum += integrand_values[i]
        else: trapezoid_sum += 2*integrand_values[i] 

     integral = c * trapezoid_sum
     return integral


print trapezoid_integrator(0, math.pi / 4, 100000)

Simple Simpson's rule solution:

import math 
sum = 0
inc = (math.pi/4)
for i in range(1000):
    cos1 = math.cos(i*inc)**100
    cos2 = math.cos((i+1) * inc)**100
    cos = (cos1 + cos2)/2
    sum += cos * inc
print(sum)