In [1]:
import sympy as sym
import numpy as np
import matplotlib.pyplot as pp
import IPython.display as ipd
import scipy.interpolate as sci
import scipy.linalg
In [2]:
def show_expr( expr, reps = [ "_{i}", "_{i + 1}", "_{i + 2}", "_{i + 3}", "_{i + 4}" ] ):
    
    ltx = sym.latex( expr )
    for i in range( len( reps ) ):
       ltx = ltx.replace( "_{%i}" % ( i, ), reps[ i ] ) 
    
    ipd.display(
        ipd.Math(
            ltx
        )
    )
In [3]:
# Ordem da spline:

k = 3

# Vamos criar as variáveis simbólicas:

( x, x_i, h_i, y_i, v_i ) = sym.symbols( 'x, x_0, h_0, y_0, v_0' )
In [4]:
# Coeficientes dos polinômios

As = []

for i in range( k ):
    
    As.append( [] )    
    for j in range( k + 1 ):
        
        As[ i ].append( sym.symbols( '%c' % ( chr( 65 + j ), ) + '_%i' % ( i, ) ) )
In [5]:
Ps = [ 0, 0 ]

for j in range( len( As[ 0 ] ) ):
    
    Ps[ 0 ] += As[ 0 ][ j ] * ( ( x - x_i ) ** j )
    Ps[ 1 ] += As[ 1 ][ j ] * ( ( x - ( x_i + h_i ) ) ** j )
In [6]:
interp = [
    sym.Equality( Ps[ 0 ].subs( x, x_i ), y_i ),
    sym.Equality( Ps[ 1 ].subs( x, x_i + h_i ), y_i + v_i )
]

for eq in interp:
    show_expr( eq )
$$A_{i} = y_{i}$$
$$A_{i + 1} = v_{i} + y_{i}$$
In [7]:
conti = []

for j in range( k ):
    
    conti.append( 
        sym.Equality(
            Ps[ 0 ].diff( x, j ).subs( x, x_i + h_i ),
            Ps[ 1 ].diff( x, j ).subs( x, x_i + h_i )
        )
    )

for eq in conti:
    show_expr( eq )
$$A_{i} + B_{i} h_{i} + C_{i} h_{i}^{2} + D_{i} h_{i}^{3} = A_{i + 1}$$
$$B_{i} + 2 C_{i} h_{i} + 3 D_{i} h_{i}^{2} = B_{i + 1}$$
$$2 C_{i} + 6 D_{i} h_{i} = 2 C_{i + 1}$$
In [8]:
tmp = conti[ 0 ]

tmp = sym.Equality(
    tmp.lhs - As[ 0 ][ 0 ],
    tmp.rhs - As[ 0 ][ 0 ]
)

tmp = tmp.subs( interp[ 0 ].lhs, interp[ 0 ].rhs )

tmp = tmp.subs( interp[ 1 ].lhs, interp[ 1 ].rhs )

tmp = sym.Equality(
    tmp.lhs / h_i,
    tmp.rhs / h_i
).expand()

conti[ 0 ] = tmp

for eq in conti:
    show_expr( eq )
$$B_{i} + C_{i} h_{i} + D_{i} h_{i}^{2} = \frac{v_{i}}{h_{i}}$$
$$B_{i} + 2 C_{i} h_{i} + 3 D_{i} h_{i}^{2} = B_{i + 1}$$
$$2 C_{i} + 6 D_{i} h_{i} = 2 C_{i + 1}$$
In [9]:
alpha = conti[ 0 ].lhs.coeff( As[ 0 ][ 3 ] ) / \
        conti[ 1 ].lhs.coeff( As[ 0 ][ 3 ] )
EqC = sym.Equality(
    conti[ 0 ].lhs - alpha * conti[ 1 ].lhs,
    conti[ 0 ].rhs - alpha * conti[ 1 ].rhs
)
EqC = sym.Equality(
    EqC.lhs * 3 / h_i,
    EqC.rhs * 3 / h_i
).expand()
EqC = sym.Equality(
    EqC.lhs - ( EqC.lhs - As[ 0 ][ 2 ] ),
    EqC.rhs - ( EqC.lhs - As[ 0 ][ 2 ] )
)

show_expr( EqC )
$$C_{i} = - \frac{2 B_{i}}{h_{i}} - \frac{B_{i + 1}}{h_{i}} + \frac{3 v_{i}}{h_{i}^{2}}$$
In [10]:
alpha = conti[ 0 ].lhs.coeff( As[ 0 ][ 2 ] ) / \
        conti[ 1 ].lhs.coeff( As[ 0 ][ 2 ] )
EqD = sym.Equality(
    conti[ 0 ].lhs - alpha * conti[ 1 ].lhs,
    conti[ 0 ].rhs - alpha * conti[ 1 ].rhs
)
EqD = sym.Equality(
    EqD.lhs * -2 / ( h_i ** 2 ),
    EqD.rhs * -2 / ( h_i ** 2 )
).expand()
EqD = sym.Equality(
    EqD.lhs - ( EqD.lhs - As[ 0 ][ 3 ] ),
    EqD.rhs - ( EqD.lhs - As[ 0 ][ 3 ] )
)

show_expr( EqD )
$$D_{i} = \frac{B_{i}}{h_{i}^{2}} + \frac{B_{i + 1}}{h_{i}^{2}} - \frac{2 v_{i}}{h_{i}^{3}}$$
In [11]:
EqCNext = EqC.subs( As[ 1 ][ 1 ], As[ 2 ][ 1 ] )\
             .subs( As[ 0 ][ 1 ], As[ 1 ][ 1 ] )\
             .subs( As[ 0 ][ 2 ], As[ 1 ][ 2 ] )\
             .subs( h_i, sym.symbols( 'h_1' ) )\
             .subs( v_i, sym.symbols( 'v_1' ) )


show_expr( EqC )
show_expr( EqCNext )
show_expr( EqD )
show_expr( conti[ -1 ] )
$$C_{i} = - \frac{2 B_{i}}{h_{i}} - \frac{B_{i + 1}}{h_{i}} + \frac{3 v_{i}}{h_{i}^{2}}$$
$$C_{i + 1} = - \frac{2 B_{i + 1}}{h_{i + 1}} - \frac{B_{i + 2}}{h_{i + 1}} + \frac{3 v_{i + 1}}{h_{i + 1}^{2}}$$
$$D_{i} = \frac{B_{i}}{h_{i}^{2}} + \frac{B_{i + 1}}{h_{i}^{2}} - \frac{2 v_{i}}{h_{i}^{3}}$$
$$2 C_{i} + 6 D_{i} h_{i} = 2 C_{i + 1}$$
In [12]:
EqB = conti[ -1 ].subs( EqC.lhs, EqC.rhs )\
                 .subs( EqCNext.lhs, EqCNext.rhs )\
                 .subs( EqD.lhs, EqD.rhs )\
                 .expand()
            
EqB = sym.Equality(
    EqB.lhs + 4 * As[ 1 ][ 1 ] / sym.symbols( 'h_1' ),
    EqB.rhs + 4 * As[ 1 ][ 1 ] / sym.symbols( 'h_1' )
)
EqB = sym.Equality(
    EqB.lhs + 2 * As[ 2 ][ 1 ] / sym.symbols( 'h_1' ),
    EqB.rhs + 2 * As[ 2 ][ 1 ] / sym.symbols( 'h_1' )
)
EqB = sym.Equality(
    EqB.lhs + 6 * v_i / ( h_i ** 2 ),
    EqB.rhs + 6 * v_i / ( h_i ** 2 )
)
EqB = sym.Equality(
    EqB.lhs / 2,
    EqB.rhs / 2
)
EqB = sym.Equality(
    EqB.lhs.collect( As[ 1 ][ 1 ] ),
    EqB.rhs
)

show_expr( conti[ -1 ] )
show_expr( EqB )
$$2 C_{i} + 6 D_{i} h_{i} = 2 C_{i + 1}$$
$$\frac{B_{i}}{h_{i}} + B_{i + 1} \left(\frac{2}{h_{i + 1}} + \frac{2}{h_{i}}\right) + \frac{B_{i + 2}}{h_{i + 1}} = \frac{3 v_{i + 1}}{h_{i + 1}^{2}} + \frac{3 v_{i}}{h_{i}^{2}}$$

$$ \begin{pmatrix} \displaystyle\frac 1{h_0} & \displaystyle\frac2{h_0} + \displaystyle\frac2{h_1} & \displaystyle\frac1{h_1}\\ \\ & \displaystyle\frac 1{h_1} & \displaystyle\frac2{h_1} + \displaystyle\frac2{h_2} & \displaystyle\frac1{h_2}\\ \\ && \displaystyle\frac 1{h_2} & \displaystyle\frac2{h_2} + \displaystyle\frac2{h_3} & \displaystyle\frac1{h_3}\\ \\ &&\qquad\ddots&\qquad\ddots & \qquad\ddots\\ \\ &&& \displaystyle\frac 1{h_{n - 2}} & \displaystyle\frac2{h_{n - 2}} + \displaystyle\frac2{h_{n - 1}} & \displaystyle\frac1{h_{n - 1}} \end{pmatrix} \begin{pmatrix} B_0\\ B_1\\ B_2\\ \vdots\\ B_n \end{pmatrix} = \begin{pmatrix} \displaystyle\frac{3v_0}{h_0^2} + \displaystyle\frac{3v_1}{h_1^2}\\ \\ \displaystyle\frac{3v_1}{h_1^2} + \displaystyle\frac{3v_2}{h_2^2}\\ \\ \vdots\\ \\ \displaystyle\frac{3v_{n - 2}}{h_{n - 2}^2} + \displaystyle\frac{3v_{n - 1}}{h_{n - 1}^2} \end{pmatrix} $$

In [13]:
EqNE = sym.Equality(
    -Ps[ 0 ].diff( x, 2 ).subs( x, x_i ) + 6 * v_i / h_i ** 2,
                                         + 6 * v_i / h_i ** 2
).subs( EqC.lhs, EqC.rhs )

show_expr( EqNE, [] )
$$\frac{4 B_{0}}{h_{0}} + \frac{2 B_{1}}{h_{0}} = \frac{6 v_{0}}{h_{0}^{2}}$$

$$ \begin{pmatrix} \displaystyle\frac 4{h_0} & \displaystyle\frac2{h_0}\\ \\ \displaystyle\frac 1{h_0} & \displaystyle\frac2{h_0} + \displaystyle\frac2{h_1} & \displaystyle\frac1{h_1}\\ \\ & \displaystyle\frac 1{h_1} & \displaystyle\frac2{h_1} + \displaystyle\frac2{h_2} & \displaystyle\frac1{h_2}\\ \\ && \displaystyle\frac 1{h_2} & \displaystyle\frac2{h_2} + \displaystyle\frac2{h_3} & \displaystyle\frac1{h_3}\\ \\ &&\qquad\ddots&\qquad\ddots & \qquad\ddots\\ \\ &&& \displaystyle\frac 1{h_{n - 2}} & \displaystyle\frac2{h_{n - 2}} + \displaystyle\frac2{h_{n - 1}} & \displaystyle\frac1{h_{n - 1}} \end{pmatrix} \begin{pmatrix} B_0\\ B_1\\ B_2\\ \vdots\\ B_n \end{pmatrix} = \begin{pmatrix} \displaystyle\frac{6v_0}{h_0^2}\\ \\ \displaystyle\frac{3v_0}{h_0^2} + \displaystyle\frac{3v_1}{h_1^2}\\ \\ \displaystyle\frac{3v_1}{h_1^2} + \displaystyle\frac{3v_2}{h_2^2}\\ \\ \vdots\\ \\ \displaystyle\frac{3v_{n - 2}}{h_{n - 2}^2} + \displaystyle\frac{3v_{n - 1}}{h_{n - 1}^2} \end{pmatrix} $$

In [14]:
EqND = sym.Equality(
    Ps[ 0 ].diff( x, 2 ).subs( x, x_i + h_i ) + 6 * v_i / h_i ** 2,
                                              + 6 * v_i / h_i ** 2
).subs( EqC.lhs, EqC.rhs ).subs( EqD.lhs, EqD.rhs )
EqND = sym.Equality(
    EqND.lhs.expand(),
    EqND.rhs
)

show_expr( EqND, [ '_{n - 1}', '_{n}' ] )
$$\frac{2 B_{n - 1}}{h_{n - 1}} + \frac{4 B_{n}}{h_{n - 1}} = \frac{6 v_{n - 1}}{h_{n - 1}^{2}}$$

$$ \begin{pmatrix} \displaystyle\frac 4{h_0} & \displaystyle\frac2{h_0}\\ \\ \displaystyle\frac 1{h_0} & \displaystyle\frac2{h_0} + \displaystyle\frac2{h_1} & \displaystyle\frac1{h_1}\\ \\ & \displaystyle\frac 1{h_1} & \displaystyle\frac2{h_1} + \displaystyle\frac2{h_2} & \displaystyle\frac1{h_2}\\ \\ && \displaystyle\frac 1{h_2} & \displaystyle\frac2{h_2} + \displaystyle\frac2{h_3} & \displaystyle\frac1{h_3}\\ \\ &&\qquad\ddots&\qquad\ddots & \qquad\ddots\\ \\ &&& \displaystyle\frac 1{h_{n - 2}} & \displaystyle\frac2{h_{n - 2}} + \displaystyle\frac2{h_{n - 1}} & \displaystyle\frac1{h_{n - 1}}\\ \\ &&&& \displaystyle\frac 2{h_{n - 1}} & \displaystyle\frac 4{h_{n - 1}} \end{pmatrix} \begin{pmatrix} B_0\\ B_1\\ B_2\\ \vdots\\ B_n \end{pmatrix} = \begin{pmatrix} \displaystyle\frac{6v_0}{h_0^2}\\ \\ \displaystyle\frac{3v_0}{h_0^2} + \displaystyle\frac{3v_1}{h_1^2}\\ \\ \displaystyle\frac{3v_1}{h_1^2} + \displaystyle\frac{3v_2}{h_2^2}\\ \\ \vdots\\ \\ \displaystyle\frac{3v_{n - 2}}{h_{n - 2}^2} + \displaystyle\frac{3v_{n - 1}}{h_{n - 1}^2}\\ \\ \displaystyle\frac{6v_{n - 1}}{h_{n - 1}^2} \end{pmatrix} $$

In [15]:
x_nodes = np.linspace( -6.0, 6.0, 8 )
y_nodes = np.sin( x_nodes )

cs = sci.CubicSpline( x_nodes, y_nodes, bc_type = 'natural' )
# cs = sci.CubicSpline( x_nodes, y_nodes, bc_type = ( ( 2, 0.0 ), ( 2, 0.0 ) ) )
In [16]:
pp.figure( figsize = ( 12, 8 ) )

x_plot = np.linspace( -6.0, 6.0, 500 )

pp.plot( x_plot, np.sin( x_plot ) )
pp.plot( x_plot, cs( x_plot ) )

pp.show()
In [17]:
err = []
n_points = []

a = -1.0
b = 1.0

def runge( x ):
    
    return np.exp( -np.abs( x ) )
#     return 1.0 / ( 1.0 + 25.0 * ( x ** 2 ) )

x_eval = np.linspace( a, b, 500 )
runge_x_eval = runge( x_eval )

for n in range( 7, 50 ):
    
    x_nodes = np.linspace( a, b, n )
    y_nodes = runge( x_nodes )
    
    cs = sci.CubicSpline( x_nodes, y_nodes, bc_type = 'natural' )
    
    cs_x_eval = cs( x_eval )
    curr_err = np.max( np.abs( cs_x_eval - runge_x_eval ) )
    
    err.append( curr_err )
    n_points.append( n )
In [18]:
pp.plot( n_points, err )
Out[18]:
[<matplotlib.lines.Line2D at 0x7f40b0c58780>]

$$ \epsilon_n \approx C h^q = \frac C{n^q} $$

$$ \log \epsilon_n \approx \log C - q\log n $$

In [19]:
def q_estimate( n_points, err ):
    
    n_points = np.array( n_points, dtype = 'float' )
    err = np.array( err )
    
    y = np.log( err )
    
    A = np.empty( ( len( n_points ), 2 ) )
    A[ :, 0 ] = 1.0
    A[ :, 1 ] = np.log( n_points )
    AA = np.dot( A.T, A )
    
    b = np.dot( A.T, y )
    
    x = np.linalg.solve( AA, b )
    
    return -x[ 1 ]
In [20]:
print( q_estimate( n_points[ 0 : : 2 ], err[ 0 : : 2 ] ) )
print( n_points[ 0 : : 2 ] )

1.0620814922553918
[7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 33, 35, 37, 39, 41, 43, 45, 47, 49]
In [21]:
pp.figure( figsize = ( 12, 8 ) )

x_plot = np.linspace( a, b, 500 )

pp.plot( x_plot, runge( x_plot ) )
pp.plot( x_plot, cs( x_plot ) )

pp.show()
In [22]:
def cs_rep( x_nodes, y_nodes ):
    
    x = np.array( x_nodes )
    y = np.array( y_nodes )
    
    h = x[ 1 :  ] - x[ : -1 ]
    v = y[ 1 :  ] - y[ : -1 ]
    
    b = np.empty( ( len( x ) ) )
    b[  0 ] = 6.0 * v[  0 ] / ( h[  0 ] ** 2 )
    b[ -1 ] = 6.0 * v[ -1 ] / ( h[ -1 ] ** 2 )
    
    h2 = h ** 2
    tmp = v / h2
    b[ 1 : -1 ] = 3.0 * ( tmp[ : -1 ] + tmp[ 1 :  ] )
    
    M = np.empty( ( 3, len( x ) ) )

    M[ 1, 0 ] = 4.0 / h[ 0 ]
    M[ 2, 0 ] = 1.0 / h[ 0 ]
    
    for i in range( 1, len( x ) - 1 ):
        
        M[ 0, i ] = 1.0 / h[ i - 1 ]
        M[ 1, i ] = 2.0 / h[ i - 1 ] + 2.0 / h[ i ]
        M[ 2, i ] = 1.0 / h[ i ]
    
    M[  0,  1 ] += 1.0 / h[  0 ]
    M[ -1, -2 ] += 1.0 / h[ -1 ]
    
    M[ 0, -1 ] = 1.0 / h[ -1 ]
    M[ 1, -1 ] = 4.0 / h[ -1 ]
    
    B = scipy.linalg.solve_banded( ( 1, 1 ), M, b )
    
    C = -( 2.0 * B[ : -1 ] + B[ 1 : ] ) / h + 3.0 * v / h2
    D = ( B[ : -1 ] + B[ 1 : ] ) / h2 - 2.0 * v / ( h ** 3 )
    
    return ( x, y, B, C, D )
In [23]:
def cs_eval( rep, x ):
    
    # Índice do polinômio atual:
    i = 0
    
    y = []
    
    for j in range( len( x ) ):
        
        while ( x[ j ] >= rep[ 0 ][ i ] ) and ( i < len( rep[ 0 ] ) - 1 ):
            i += 1
        i -= 1
        if ( i < 0 ):
            i = 0
        
        h = x[ j ] - rep[ 0 ][ i ]
        
        tmp = rep[ 4 ][ i ]
        tmp = h * tmp + rep[ 3 ][ i ]
        tmp = h * tmp + rep[ 2 ][ i ]
        tmp = h * tmp + rep[ 1 ][ i ]
        
        y.append( tmp )
        
    return np.array( y )
In [24]:
x_nodes = np.linspace( -6.0, 6.0, 8 )
y_nodes = np.sin( x_nodes )

rep = cs_rep( x_nodes, y_nodes )
cs = sci.CubicSpline( x_nodes, y_nodes, bc_type = 'natural' )

pp.figure( figsize = ( 12, 8 ) )

x_plot = np.linspace( -6.0, 6.0, 500 )

pp.plot( x_plot, np.sin( x_plot ) )
pp.plot( x_plot, cs( x_plot ) )
pp.plot( x_plot, cs_eval( rep, x_plot ) )
pp.show()