In [8]:
%load_ext Cython
In [12]:
def f(x):
    return cos(2*x) / e^x

def Simpson(f, a, b, eps, level, level_max):
    level = level + 1
    h = b - a
    c = (a + b)/2
    
    print h, c
    print f(a), f(b), f(c)
    one_simpson = h * (f(a) + 4*f(c) + f(b))/6
    d = (a + c)/2
    e = (c + b)/2
    two_simpson = h * (f(a) + 4*f(d) + 2*f(c) + 4*f(e) + f(b))/12
    print level, eps, one_simpson, two_simpson

    if level >= level_max:
        raise ValueError("maximum level reached: " + two_simpson.str())
    else:
        if abs(two_simpson - one_simpson) < 15 * eps:
            return two_simpson + (two_simpson - one_simpson)/15
        else:
            left_simpson = Simpson(f, a, c, eps/2, level, level_max)
            right_simpson = Simpson(f, c, b, eps/2, level, level_max)
            return left_simpson + right_simpson
In [25]:
%%cython

# encoding: utf-8
# cython: profile=True
# filename: simpson.pyx

#from libc.math cimport cos, exp as e
cdef extern from "math.h":
    double cos(double)
    double exp(double)

cdef inline double f(double x):
    return cos(2*x) / exp(x)

cpdef simpson(f, double a, double b, double eps, int level, int level_max):
    level = level + 1
    cdef double h = b - a
    cdef c = (a + b)/2
    cdef double one_simpson = h * (f(a) + 4*f(c) + f(b))/6
    cdef double d = (a + c)/2
    cdef double e = (c + b)/2
    cdef double two_simpson = h * (f(a) + 4*f(d) + 2*f(c) + 4*f(e) + f(b))/12
    print level, eps, one_simpson, two_simpson

    if level >= level_max:
        raise ValueError("maximum level reached: ", level_max, two_simpson)
    else:
        if abs(two_simpson - one_simpson) < 15 * eps:
            return two_simpson + (two_simpson - one_simpson)/15
        else:
            left_simpson = simpson(f, a, c, eps/2, level, level_max)
            right_simpson = simpson(f, c, b, eps/2, level, level_max)
            return left_simpson + right_simpson
In [26]:
f
Out[26]:
<function f at 0x7f894c51c668>
In [27]:
Simpson
Out[27]:
<function Simpson at 0x7f7d2e107938>
In [28]:
simpson(f, 0, 5*pi/4, 0.5 * 10^-3, 0, 5)

1 0.0005 0.394651050365 0.138209666247
2 0.00025 0.107091467975 0.175825580865
3 0.000125 0.36273356355 0.366881289759
4 6.25e-05 0.33542953724 0.335554211958
4 6.25e-05 0.0314517525188 0.0315602938143
3 0.000125 -0.186907982685 -0.186985849085
2 0.00025 0.0311181982722 0.0281471391512
Out[28]:
0.20808808187162123
In [6]:
Simpson(f, 0, 5*pi/4, 0.5 * 10^-3, 0, 400).n()

5/4*pi 5/8*pi
1 0 -1/2*sqrt(2)*e^(-5/8*pi)
1 0.000500000000000000 -5/24*pi*(2*sqrt(2)*e^(-5/8*pi) - 1) 5/48*pi*(4*cos(15/8*pi)*e^(-15/16*pi) - 2*sqrt(-sqrt(2) + 2)*e^(-5/16*pi) - sqrt(2)*e^(-5/8*pi) + 1)
5/8*pi 5/16*pi
1 -1/2*sqrt(2)*e^(-5/8*pi) -1/2*sqrt(-sqrt(2) + 2)*e^(-5/16*pi)
2 0.000250000000000000 -5/96*pi*(4*sqrt(-sqrt(2) + 2)*e^(-5/16*pi) + sqrt(2)*e^(-5/8*pi) - 2) 5/192*pi*(8*cos(5/16*pi)*e^(-5/32*pi) + 8*cos(15/16*pi)*e^(-15/32*pi) - 2*sqrt(-sqrt(2) + 2)*e^(-5/16*pi) - sqrt(2)*e^(-5/8*pi) + 2)
5/16*pi 5/32*pi
1 -1/2*sqrt(-sqrt(2) + 2)*e^(-5/16*pi) cos(5/16*pi)*e^(-5/32*pi)
3 0.000125000000000000 5/192*pi*(8*cos(5/16*pi)*e^(-5/32*pi) - sqrt(-sqrt(2) + 2)*e^(-5/16*pi) + 2) 5/384*pi*(8*cos(5/32*pi)*e^(-5/64*pi) + 4*cos(5/16*pi)*e^(-5/32*pi) + 8*cos(15/32*pi)*e^(-15/64*pi) - sqrt(-sqrt(2) + 2)*e^(-5/16*pi) + 2)
5/32*pi 5/64*pi
1 cos(5/16*pi)*e^(-5/32*pi) cos(5/32*pi)*e^(-5/64*pi)
4 0.0000625000000000000 5/192*pi*(4*cos(5/32*pi)*e^(-5/64*pi) + cos(5/16*pi)*e^(-5/32*pi) + 1) 5/384*pi*(4*cos(5/64*pi)*e^(-5/128*pi) + 2*cos(5/32*pi)*e^(-5/64*pi) + 4*cos(15/64*pi)*e^(-15/128*pi) + cos(5/16*pi)*e^(-5/32*pi) + 1)
5/32*pi 15/64*pi
cos(5/16*pi)*e^(-5/32*pi) -1/2*sqrt(-sqrt(2) + 2)*e^(-5/16*pi) cos(15/32*pi)*e^(-15/64*pi)
4 0.0000625000000000000 5/384*pi*(2*cos(5/16*pi)*e^(-5/32*pi) + 8*cos(15/32*pi)*e^(-15/64*pi) - sqrt(-sqrt(2) + 2)*e^(-5/16*pi)) 5/768*pi*(2*cos(5/16*pi)*e^(-5/32*pi) + 8*cos(25/64*pi)*e^(-25/128*pi) + 4*cos(15/32*pi)*e^(-15/64*pi) + 8*cos(35/64*pi)*e^(-35/128*pi) - sqrt(-sqrt(2) + 2)*e^(-5/16*pi))
5/16*pi 15/32*pi
-1/2*sqrt(-sqrt(2) + 2)*e^(-5/16*pi) -1/2*sqrt(2)*e^(-5/8*pi) cos(15/16*pi)*e^(-15/32*pi)
3 0.000125000000000000 5/192*pi*(8*cos(15/16*pi)*e^(-15/32*pi) - sqrt(-sqrt(2) + 2)*e^(-5/16*pi) - sqrt(2)*e^(-5/8*pi)) 5/384*pi*(8*cos(25/32*pi)*e^(-25/64*pi) + 4*cos(15/16*pi)*e^(-15/32*pi) + 8*cos(35/32*pi)*e^(-35/64*pi) - sqrt(-sqrt(2) + 2)*e^(-5/16*pi) - sqrt(2)*e^(-5/8*pi))
5/8*pi 15/16*pi
-1/2*sqrt(2)*e^(-5/8*pi) 0 cos(15/8*pi)*e^(-15/16*pi)
2 0.000250000000000000 5/96*pi*(8*cos(15/8*pi)*e^(-15/16*pi) - sqrt(2)*e^(-5/8*pi)) 5/192*pi*(8*cos(25/16*pi)*e^(-25/32*pi) + 4*cos(15/8*pi)*e^(-15/16*pi) + 8*cos(35/16*pi)*e^(-35/32*pi) - sqrt(2)*e^(-5/8*pi))
Out[6]:
0.208088081871621
In [ ]: