In [1]:
import numpy as np
import sympy as sym
import IPython.display as ipd
import matplotlib.pyplot as pp
In [2]:
def f( x ):
    
    return np.sin( x )
In [3]:
def trapezio( f, a, b, N ):
    
    x = np.linspace( a, b, N + 1 )
    fx = f( x )
    
    retval = fx[ 0 ] + 2.0 * fx[ 1 : -1 ].sum() + fx[ -1 ]
    retval = retval * ( b - a ) / ( 2.0 * N )
    
    return retval
In [4]:
def If( a, b ):
    
    return -np.cos( b ) + np.cos( a )
In [5]:
print( If( 0, 1 ) )

0.45969769413186023
In [6]:
print( trapezio( f, 0, 1, 10 ) )

0.4593145488579764

$$ \left|\frac{(b - a)^3}{12N^2}f''( \eta )\right| \leq \epsilon $$

$$ \left|\frac{(b - a)^3}{12N^2}\right||f''( \eta )| \leq \epsilon $$

$$ M = \max_{x \in [ a, b ]}|f''(x)| $$

$$ \left|\frac{(b - a)^3}{12N^2}\right||f''( \eta )| \leq \left|\frac{(b - a)^3}{12N^2}\right|M \leq \epsilon $$

$$ \frac{(b - a)^3}{12N^2}M \leq \epsilon $$

$$ \frac{12N^2}{(b - a)^3M} \geq \frac 1\epsilon $$

$$ N^2 \geq \frac {(b - a)^3M}{12\epsilon} $$

$$ N \geq \sqrt{\frac {(b - a)^3M}{12\epsilon}} $$

In [7]:
eps = 1e-12
M = abs( np.sin( 1.0 ) )
b = 1.0
a = 0.0

N = ( ( ( b - a ) ** 3 * M ) / ( 12.0 * eps ) ) ** 0.5
In [8]:
print( N )

264806.68810912746
In [9]:
print( trapezio( f, 0, 1, 264807 ) - If( 0.0, 1.0 ) )

-5.461742169643458e-13
In [10]:
def create_ell( i, points, x = sym.symbols( 'x' ) ):
    
    ell = 1
    for j in range( len( points ) ):
        
        if j != i:
            ell = ell * \
               ( x - points[ j ] ) / ( points[ i ] - points[ j ] )
    
    return ell
In [11]:
def show_expr( expr ):
    
    ipd.display( ipd.Math( sym.latex( expr ) ) )
In [12]:
( x, a, b ) = sym.symbols( 'x a b' )

def create_points( n, a = sym.symbols( 'a' ), b = sym.symbols( 'b' ) ):
    
    delta = ( b - a ) / n
    retval = []
    for i in range( n + 1 ):
        retval.append( a + delta * i )
        
    return retval

points = create_points( 2 )
print( points )

[a, a/2 + b/2, b]
In [13]:
ells = [ create_ell( i, points ) for i in range( len( points ) ) ]
In [14]:
As = [
    sym.factor( sym.integrate( ell, ( x, a, b ) ) ) for
    ell in ells
]
In [15]:
print( As )

[-(a - b)/6, -2*(a - b)/3, -(a - b)/6]
In [16]:
print(
    sym.factor(
        sym.integrate(
            ( x - a ) * ( x - ( a + b ) / 2 ) * ( x - b ),
            ( x, a, b )
        )
    )
)

0
In [17]:
show_expr(
    sym.factor(
        sym.integrate(
            ( x - a ) * ( ( x - ( a + b ) / 2 ) ** 2 ) * ( x - b ),
            ( x, a, b )
        )
    )
)
$$\frac{1}{120} \left(a - b\right)^{5}$$
In [18]:
24 * 120 == 90 * 2 ** 5
Out[18]:
True

$$ \frac{(b - a)^5}{180N^4}f^{(4)}( \eta ) $$

$$ M := \max_{x \in [ a, b ]}|f^{( 4 )}( x )| $$

$$ \left|\frac{(b - a)^5}{180N^4}f^{(4)}( \eta )\right| = \left|\frac{(b - a)^5}{180N^4}\right||f^{(4)}( \eta )| $$

$$ \left|\frac{(b - a)^5}{180N^4}f^{(4)}( \eta )\right| = \left|\frac{(b - a)^5}{180N^4}\right||f^{(4)}( \eta )| \le \frac{(b - a)^5}{180N^4}M \le \epsilon $$

$$ \frac{(b - a)^5}{180N^4}M \le \epsilon $$

$$ \frac{(b - a)^5}{180\epsilon}M \le N^4 $$

$$ N \ge \sqrt[4]{\frac{(b - a)^5}{180\epsilon}M} $$

In [19]:
M = np.sin( 1 )
eps = 1e-12
a = 0.0
b = 1.0

N = ( ( b - a ) ** 5 * M / ( 180 * eps ) ) ** 0.25
In [20]:
print( N )

261.48191690325257
In [21]:
def simpson( f, a, b, N ):
    
    x = np.linspace( a, b, N + 1 )
    fx = f( x )
    
    retval =       fx[ 0 ] + \
             4.0 * fx[ 1 : -1 : 2 ].sum() + \
             2.0 * fx[ 2 : -1 : 2 ].sum() + \
                   fx[ -1 ]
    retval = retval * ( b - a ) / ( 3.0 * N )
    
    return retval
In [22]:
print( simpson( f, a, b, 262 ) - If( a, b ) )

5.420108806220014e-13
In [23]:
X = sym.symbols( 'x' )
n = 2
polys = [ X ** i for i in range( n + 2 ) ]
print( polys )

[1, x, x**2, x**3]
In [24]:
def inner_prod( f, g, X = sym.symbols( 'x' ) ):
    
    return sym.integrate( f * g, ( X, -1, 1 ) )
In [25]:
def gram_schmidt( x, inner_prod ):
    
    n = len( x )
    
    v = []
    for k in range( n ):
        
        tmp = x[ k ]
        for i in range( k - 1 ):
            tmp -= inner_prod( x[ k ], v[ i ] ) * v[ i ]
        
        v.append( tmp / ( sym.sqrt( inner_prod( tmp, tmp ) ) ) )
        
    return v
In [26]:
qs = gram_schmidt( polys, inner_prod )
print( qs )

[sqrt(2)/2, sqrt(6)*x/2, 3*sqrt(10)*(x**2 - 1/3)/4, 5*sqrt(14)*(x**3 - 3*x/5)/4]
In [27]:
points_g = sym.solve( qs[ -1 ] )
points_g.sort()
print( points_g )

[-sqrt(15)/5, 0, sqrt(15)/5]
In [28]:
ells_g = [ create_ell( i, points_g ) for i in range( len( points_g ) ) ]
In [29]:
print( ells_g )

[5*x*(x - sqrt(15)/5)/6, -5*(x - sqrt(15)/5)*(x + sqrt(15)/5)/3, 5*x*(x + sqrt(15)/5)/6]
In [30]:
As_g = [
    sym.factor( sym.integrate( ell, ( X, -1, 1 ) ) ) for
    ell in ells_g
]
In [31]:
print( As_g )

[5/9, 8/9, 5/9]
In [32]:
sym.factor(
    sym.integrate( 
        ( ( X - points_g[ 0 ] ) * ( X - points_g[ 1 ] ) * ( X - points_g[ 2 ] ) ) ** 2,
        ( X, -1, 1 )
    )
)
Out[32]:
8/175
In [33]:
print( inner_prod( qs[ 3 ], X ** 2 ) )

0
In [34]:
def gauss_quad( n ):
    
    x = sym.symbols( 'x' )
    
    polys = [ x ** i for i in range( n + 2 ) ]
    
    q = gram_schmidt( polys, inner_prod )[ -1 ]
    
    nodes = sym.solve( q )
    nodes.sort()
    
    As =[
        sym.factor( sym.integrate( create_ell( i, nodes ), ( x, -1, 1 ) ) )
        for i in range( len( nodes ) )
    ]
    
    return ( nodes, As )
In [35]:
( nodes, As ) = gauss_quad( 2 )
print( nodes )
print( As )

[-sqrt(15)/5, 0, sqrt(15)/5]
[5/9, 8/9, 5/9]
In [36]:
def gauss_rep( f, a, b, M, nodes, As ):
    
    nodes = np.array( [ float( node ) for node in nodes ] )
    As = np.array( [ float( A ) for A in As ] )
    
    h = ( b - a ) / ( 2.0 * M )
    
    retval = 0.0
    
    for i in range( M ):
        
        c = i * 2.0 * h + h
        
        retval += ( f( nodes * h + c ) * As ).sum()
        
    return retval * h
In [37]:
print( gauss_rep( f, a, b, 9, nodes, As ) - If( a, b ) )

4.2932324362254803e-13

$$ \left( \frac{b - a}2 \right)^{7}\frac{f^{(6)}( \eta )}{15750M^6} $$

$$ \left| \left( \frac{b - a}2 \right)^{7}\frac{f^{(6)}( \eta )}{15750M^6} \right| = \left( \frac{b - a}2 \right)^{7}\frac{|f^{(6)}( \eta )|}{15750M^6} $$

$$ U := \max_{x \in [ a, b ]}|f^{(6)}( x )| $$

$$ \left| \left( \frac{b - a}2 \right)^{7}\frac{f^{(6)}( \eta )}{15750M^6} \right| = \left( \frac{b - a}2 \right)^{7}\frac{|f^{(6)}( \eta )|}{15750M^6} \leq \left( \frac{b - a}2 \right)^{7}\frac{U}{15750M^6} $$

$$ \left( \frac{b - a}2 \right)^{7}\frac{U}{15750M^6} \leq \epsilon $$

$$ M^6 \geq \left( \frac{b - a}2 \right)^{7}\frac{U}{15750\epsilon} $$

$$ M \geq \sqrt[6]{\left( \frac{b - a}2 \right)^{7}\frac{U}{15750\epsilon}} $$

In [38]:
U = np.sin( b )
print( ( ( ( b - a ) / 2.0 ) ** 7 * U / ( 15750 * eps )  ) ** ( 1.0 / 6.0 )  )

8.644862756922969
In [39]:
b = 2.0
M = int( np.ceil( ( ( ( b - a ) / 2.0 ) ** 7 * U / ( 15750 * eps )  ) ** ( 1.0 / 6.0 ) ) )
print( M )

20
In [40]:
print( gauss_rep( f, a, b, M, nodes, As ) - If( a, b ) )

7.027711745877241e-13