In [1]:
import numpy as np
import sympy as sym
import matplotlib.pyplot as pp
import IPython.display as ipd
In [2]:
def euler( f, t_0, y_0, t_M, M ):
    
    delta_t = ( t_M - t_0 ) / M
    
    y = [ y_0 ]
    t = [ t_0 ]
    
    for i in range( M ):
        
        y_k = y[ -1 ]
        t_k = t[ -1 ]
        y.append( y_k + delta_t * f( t_k, y_k ) )
        t.append( t_k + delta_t )
    
    return ( t, y )

$$\Large \theta''( t ) = -\frac{g}{L}\sin\theta( t ) $$ $$\Large \theta''( t ) = f( t, \theta( t ), \theta'( t ) ) $$ $$\Large f( t, y_2, y_1 ) = -\frac{g}{L}\sin y_2 $$ $g = L$

$$ \begin{pmatrix}% f( t, y_n, y_{n - 1}, \dots, y_1 )\\ y_1\\ y_2\\ \vdots\\ y_{n - 1} \end{pmatrix} $$

In [3]:
def f( t, y ):
    
    return np.array( [ -np.sin( y[ 1 ] ), y[ 0 ] ] )

def Jf( t, y ):
    
    return np.array(
        [
            [ 0.0, 0.0, -np.cos( y[ 1 ] ) ],
            [ 0.0, 1.0, 0.0 ]
        ]
    )
In [4]:
A = 1.0

y_0 = np.array( [ 0.0, A ] )
t_0 = 0.0
In [5]:
( t, y ) = euler( f, t_0, y_0, 4 * np.pi, 100000 )
In [6]:
pp.plot( t, [ v[ 1 ] for v in y ] )
pp.plot( t, [ A * np.cos( x ) for x in t ] )
pp.show()
In [7]:
T = sym.symbols( 't' )
F = sym.Function( 'f' )
Y = sym.Function( 'y' )
In [8]:
def show_expr( expr ):
    
    ipd.display( ipd.Math( sym.latex( expr ) ) )
In [9]:
show_expr( sym.diff( F( T, Y( T ) ), T ) )
$$\frac{d}{d t} y{\left (t \right )} \left. \frac{\partial}{\partial \xi_{2}} f{\left (t,\xi_{2} \right )} \right|_{\substack{ \xi_{2}=y{\left (t \right )} }} + \left. \frac{\partial}{\partial \xi_{1}} f{\left (\xi_{1},y{\left (t \right )} \right )} \right|_{\substack{ \xi_{1}=t }}$$

$$\Large \frac{\mathrm d}{\mathrm dt}f\bigl( t, y( t ) \bigr) = y'(t)f_y\bigl( t, y( t ) \bigr) + f_x\bigl( t, y( t ) \bigr) $$

$$\Large \frac{\mathrm d}{\mathrm dt}f\bigl( t, y( t ) \bigr) = f\bigl( t, y( t ) \bigr)f_y\bigl( t, y( t ) \bigr) + f_x\bigl( t, y( t ) \bigr) $$

$$\Large \frac{\mathrm d}{\mathrm dt}f\bigl( t, y( t ) \bigr) = ff_y + f_x $$

In [10]:
def taylor2( f, Jf, t_0, y_0, t_M, M ):
    
    delta_t = ( t_M - t_0 ) / M
    
    y = [ y_0 ]
    t = [ t_0 ]
    
    for i in range( M ):
        
        y_k = y[ -1 ]
        t_k = t[ -1 ]
        f_k = f( t_k, y_k )
        Jf_k = Jf( t_k, y_k )
        
        tmp = y_k + delta_t * f_k
        tmp += delta_t ** 2 * 0.5 * ( Jf_k[ :, 0 ] + np.matmul( Jf_k[ :, 1 :  ], f_k ) )
        
        y.append( tmp )
        t.append( t_k + delta_t )
    
    return ( t, y )
In [11]:
( t_T, y_T ) = taylor2( f, Jf, t_0, y_0, 4 * np.pi, 100 )
In [12]:
pp.plot( t, [ v[ 1 ] for v in y ] )
pp.plot( t_T, [ v[ 1 ] for v in y_T ] )
pp.show()
In [13]:
def heun( f, t_i, y_i, t_f, M ):
    
    h = ( t_f - t_i ) / M
    hh = h / 2.0
    
    y = [ y_i ]
    t = [ t_i ]
    
    for i in range( M ):
        
        t_k = t[ -1 ]
        y_k = y[ -1 ]
        kappa_1 = f( t_k, y_k )
        kappa_2 = f( t_k + h, y_k + h * kappa_1 )
        
        y.append( y_k + hh * ( kappa_1 + kappa_2 ) )
        t.append( t_k + h )
        
    return ( t, y )
In [14]:
( t_H, y_H ) = heun( f, t_0, y_0, 4 * np.pi, 100 )
pp.plot( t, [ v[ 1 ] for v in y ] )
pp.plot( t_H, [ v[ 1 ] for v in y_H ] )
pp.show()
In [15]:
print( max( [ abs( y_T[ i ][ 1 ] - y_H[ i ][ 1 ] ) for i in range( len( y_H ) ) ] ) )

0.004908643623912606
In [16]:
def ponto_medio( f, t_i, y_i, t_f, M ):
    
    h = ( t_f - t_i ) / M
    hh = h / 2.0
    
    y = [ y_i ]
    t = [ t_i ]
    
    for i in range( M ):
        
        t_k = t[ -1 ]
        y_k = y[ -1 ]
        kappa_1 = f( t_k, y_k )
        kappa_2 = f( t_k + hh, y_k + hh * kappa_1 )
        
        y.append( y_k + h * kappa_2 )
        t.append( t_k + h )
        
    return ( t, y )
In [17]:
( t_M, y_M ) = ponto_medio( f, t_0, y_0, 4 * np.pi, 100 )
pp.plot( t, [ v[ 1 ] for v in y ] )
pp.plot( t_M, [ v[ 1 ] for v in y_M ] )
pp.show()
In [18]:
print( max( [ abs( y_T[ i ][ 1 ] - y_M[ i ][ 1 ] ) for i in range( len( y_M ) ) ] ) )

0.0023648442315396256
In [19]:
b_2 = 0.0
b_3 = 3.0 / 4.0
b_1 = 1.0 - ( b_2 + b_3 )

$$ \begin{split} b_1 + b_2 + b_3 & {}= 1\\ b_2c_2 + b_3c_3 & {}= \frac12\\ b_2c_2^2 + b_3c_3^2 & {}= \frac13\\ b_3a_{3,2}c_2 & {}= \frac16\\ a_{2,1} & {}= c_2 \\ a_{3,1} + a_{3,2} & {}= c_3 \end{split} $$

$$ c_3 = \frac{\frac12 - b_2c_2}{b_3} $$

$$ b_2c_2^2 + b_3\left(\frac{\frac12 - b_2c_2}{b_3}\right)^2 = \frac13 $$

$$ b_2c_2^2 + \frac{(\frac12 - b_2c_2)^2}{b_3} = \frac13 $$

$$ b_3b_2c_2^2 + \frac14 - b_2c_2 + (b_2c_2)^2 = \frac {b_3}3 $$

$$ (b_2 + b_3)b_2c_2^2 - b_2c_2 + \frac14 - \frac{b_3}3 = 0 $$

Como $b_2 = 0$: $$ \frac14 - \frac{b_3}3 = 0 $$ ou seja $$ b_3 = \frac34. $$

In [20]:
c_3 = 2.0 / 3.0
c_2 = 1.0 / 3.0
In [21]:
a_21 = c_2
a_31 = 0.0
a_32 = c_3 - a_31
In [22]:
print( b_1 + b_2 + b_3 )

1.0
In [23]:
print( b_2 * c_2 + b_3 * c_3 )

0.5
In [24]:
print( b_2 * c_2 ** 2 + b_3 * c_3 ** 2 )

0.3333333333333333
In [25]:
print( b_3 * a_32 * c_2 )

0.16666666666666666
In [26]:
print( a_21 - c_2 )

0.0
In [27]:
print( a_31 + a_32 - c_3 )

0.0
In [28]:
def rk3( f, t_i, y_i, t_f, M ):
    
    h = ( t_f - t_i ) / M
    hh = h / 2.0
    
    y = [ y_i ]
    t = [ t_i ]
    
    for i in range( M ):
        
        t_k = t[ -1 ]
        y_k = y[ -1 ]
        
        kappa_1 = f( t_k, y_k )
        kappa_2 = f( t_k + c_2 * h, y_k + h * a_21 * kappa_1 )
        kappa_3 = f( t_k + c_3 * h, y_k + h * ( a_31 * kappa_1 + a_32 * kappa_2 ) )
        
        y.append( y_k + h * ( b_1 * kappa_1 + b_2 * kappa_2 + b_3 * kappa_3 ) )
        t.append( t_k + h )
        
    return ( t, y )
In [29]:
M = 50
( t_3, y_3 ) = rk3( f, t_0, y_0, 4 * np.pi, M )
pp.plot( t, [ v[ 1 ] for v in y ] )
pp.plot( t_3, [ v[ 1 ] for v in y_3 ] )
( t_HH, y_HH ) = heun( f, t_0, y_0, 4 * np.pi, M )
pp.plot( t_HH, [ v[ 1 ] for v in y_HH ], 'k' )
pp.show()
In [32]:
( h, x ) = sym.symbols( 'h, x' )
print( sym.integrate( ( x - h ) / ( 0 - h ), ( x, 0, h ) ) )
print( sym.integrate( ( x - 0 ) / ( h - 0 ), ( x, 0, h ) ) )

h/2
h/2
In [33]:
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 [43]:
def par_ie( n, b, h = sym.symbols( 'h' ) ):
    
    nodes_e = [ -i * h for i in range( n + 1 ) ]
    nodes_i = [ -i * h for i in range( -1, n + 1 ) ]
    
    ells_e = [ create_ell( i, nodes_e ) for i in range( len( nodes_e ) ) ]
    ells_i = [ create_ell( i, nodes_i ) for i in range( len( nodes_i ) ) ]
    
    x_I = -b * h
    x_F = h
    
    As_e = [
        sym.factor(
            sym.integrate( ell, ( x, x_I, x_F ) ) / h
        ) for ell in ells_e
    ]

    As_i = [
        sym.factor(
            sym.integrate( ell, ( x, x_I, x_F ) ) / h
        ) for ell in ells_i
    ]
    
    return ( As_i, As_e )
In [44]:
print( par_ie( 0, 0 ) )

([1/2, 1/2], [1])
In [45]:
print( par_ie( 1, 0 ) )

([5/12, 2/3, -1/12], [3/2, -1/2])
In [56]:
def pc( f, t_i, y_i, t_f, M, As_e, As_i, b_e, b_i, p ):
    
    As_e = [ float( A ) for A in As_e ]
    As_i = [ float( A ) for A in As_i ]
    
    n_e = len( As_e )
    n_i = len( As_i )
    
    h = ( t_f - t_i ) / M
    
    Ys = y_i
    Ts = [ t_i + i * h for i in range( len( Ys ) ) ]
    Fs = [ f( Ts[ i ], Ys[ i ] ) for i in range( len( Ys ) ) ]
    
    for k in range( M - len( Ys ) + 1 ):
        
        # Método explícito
        it = 0.0
        for i in range( n_e ):
            
            it += As_e[ i ] * Fs[ -( i + 1 ) ]
            
        it *= h
        it += Ys[ -( b_e + 1 ) ]
        
        # Método implícito
        c = 0.0
        for i in range( n_i - 1 ):
            
            c += As_i[ i ] * Fs[ -( i + 1 ) ]
            
        c *= h
        c += Ys[ -( b_i + 1 ) ]
        
        t_next = Ts[ -1 ] + h
        
        for i in range( p ):
            
            F = f( t_next, it )
            it = c + h * As_i[ -1 ] * F
            
        Fs.append( F )
        Ts.append( t_next )
        Ys.append( it )
        
    return ( Ts, Ys )
In [59]:
M = 150

pp.plot( t, [ v[ 1 ] for v in y ] )

( t_HH, y_HH ) = heun( f, t_0, y_0, 4 * np.pi, M )
pp.plot( t_HH, [ v[ 1 ] for v in y_HH ], 'k' )

( t_E, y_E ) = euler( f, t_0, y_0, 4 * np.pi, M )
pp.plot( t_E, [ v[ 1 ] for v in y_E ], 'g' )

( t_PC, y_PC ) = pc( f, t_0, [ y_0 ], 4 * np.pi, M, [ 1.0 ], [ 0.5, 0.5 ], 0, 0, 3 )
pp.plot( t_PC, [ v[ 1 ] for v in y_PC ] )

pp.show()