Fixemos $\boldsymbol c = ( c_1, c_2, c_3 )^T$ e $l > 0$, então nossa superfície será formada pelos pontos $\bigl( x_1, x_2, \gamma( x_1, x_2 ) \bigr)^T$ onde $$ \gamma( x_1, x_2 ) = l \bigl( ( x_1 - c_1 )^2 + ( x_2 - c_2 )^2 \bigr) + c_3. $$
O gradiente dessa funcão é dado por $$ \nabla \gamma( x_1, x_2 ) = 2l \begin{pmatrix} x_1 - c_1\\ x_2 - c_2 \end{pmatrix} $$ e sua hessiana por $$ \nabla^2 \gamma( x_1, x_2 ) = 2l \begin{pmatrix} 1 & 0\\ 0 & 1 \end{pmatrix}. $$
# É útil importar as bibliotecas de antemão.
import numpy as np
import matplotlib.pyplot as pp
from mpl_toolkits.mplot3d import *
# Vamos implementar as funções relacionadas à curva que serão utilizadas posteriormente.
c = [ 1.0, 1.0, 2.0 ]
l = 0.25
def gamma( x ):
return c[ 2 ] + l * ( ( x[ 0 ] - c[ 0 ] ) ** 2 + ( x[ 1 ] - c[ 1 ] ) ** 2 )
def grad_gamma( x ):
retval = np.empty( x.shape )
retval[ 0 ] = l * 2.0 * ( x[ 0 ] - c[ 0 ] )
retval[ 1 ] = l * 2.0 * ( x[ 1 ] - c[ 1 ] )
return retval
def hess_gamma( x ):
retval = np.empty( ( x.shape[ 0 ], x.shape[ 0 ] ) )
retval[ 0, 0 ] = 2.0 * l
retval[ 0, 1 ] = 0.0
retval[ 1, 0 ] = 0.0
retval[ 1, 1 ] = 2.0 * l
return retval
# Vamos visualizar nossa superfície
n_pts = 30
x = np.outer( np.linspace( -1.0 + c[ 0 ], 1.0 + c[ 0 ], n_pts ), np.ones( n_pts ) )
y = np.outer( np.linspace( -1.0 + c[ 1 ], 1.0 + c[ 1 ], n_pts ), np.ones( n_pts ) ).T
z = np.empty( ( n_pts, n_pts ) )
for i in range( n_pts ):
for j in range( n_pts ):
z[ i, j ] = gamma( [ x[ i, j ], y[ i, j ] ] )
fig = pp.figure( figsize = ( 12, 12 ) )
ax = fig.gca( projection = '3d' )
ax.view_init( elev = 10.0, azim = 195.0 )
surf = ax.plot_surface( x, y, z )
pp.show()
Lembremo-nos que nossa função de iteração será $$ \boldsymbol \phi( \boldsymbol x ) = -\gamma( \boldsymbol x )\nabla \gamma( \boldsymbol x ). $$ Ou seja, cada uma de suas componentes é $$ \phi_i( \boldsymbol x ) = -\gamma( \boldsymbol x )\frac{\partial \gamma}{\partial x_i}( \boldsymbol x ). $$
Desta forma, podemos calcular cada componente da jacobiana $$ J\boldsymbol \phi( \boldsymbol x )_{i, j} = -\frac{\partial \gamma}{\partial x_j}( \boldsymbol x )\frac{\partial \gamma}{\partial x_i}( \boldsymbol x ) - \gamma( \boldsymbol x )\frac{\partial^2 \gamma}{\partial x_i\partial x_j}( \boldsymbol x ). $$
Em notação matricial temos: $$ J\phi( \boldsymbol x ) = -\nabla \gamma( \boldsymbol x )\nabla \gamma( \boldsymbol x )^T - \gamma( \boldsymbol x )\nabla^2 \gamma( \boldsymbol x ) $$
Nota: na aula o método acelerado não convergiu por causa da troca de sinal $\boldsymbol \phi( \boldsymbol x ) = \gamma( \boldsymbol x )\nabla \gamma( \boldsymbol x )$, corrigindo esse erro o resultado foi o que podem ver a seguir. É uma situação que mostra que errar "apenas um sinal" pode ser bem problemático.
# Implementamos a função de iteração e sua Jacobiana:
def phi( x, **kwargs ):
return -gamma( x, **kwargs ) * grad_gamma( x, **kwargs )
def jaco_phi( x, **kwargs ):
gg = grad_gamma( x, **kwargs )
return -np.outer( gg, gg ) - gamma( x, **kwargs ) * hess_gamma( x, **kwargs )
x0 = np.array( [ 0.1, 0.1 ] )
x = x0
for i in range( 10 ):
x = phi( x )
print( x )
Note que obtivemos divergência, mas talvez o método acelerado funcione.
# Implementação da função de iteração modificada. Observe que o
# produto da inversa por um vetor é computado através da solução
# de um sistema linear. Para problemas em duas dimensões a diferença
# é pequena, mas em dimensões mais elevadas o cálculo dessa forma
# economiza um tempo considerável de computação.
def psi( x ):
M = jaco_phi( x ) - np.eye( x.shape[ 0 ] )
d = x - phi( x )
s = np.linalg.solve( M, d )
return x + s
x = [ x0 ]
stop = False
while not stop:
x.append( psi( x[ -1 ] ) )
stop = ( ( x[ -2 ] - x[ -1 ] ) ** 2 ).sum() ** 0.5 <= 1e-14
print( x[ -1 ] )
print( len( x ) )
Como vimos, convergimos em apenas 6 iteraçoes mostrando que a função de iteração modificada pode ser útil mesmo quando a função de iteração original não convergiria, corrigindo um problema ainda pior do que a baixa velocidade de convergência.
Agora sabemos qual o ponto fixo onde devemos computar a jacobiana da função original para verificar o critério de convergência.
#Calculamos a jacobiana nesse ponto
print( jaco_phi( x[ -1 ] ) )
Como esperado, uma vez que não houve convergência da sequência gerada pela função de iteração $\phi$, nem o critério das linhas nem o das colunas é satisfeito no ponto fixo para o qual desejamos convergir.
Agora vamos visualizar o resultado que acabamos de obter. Faremos uma imagem da função distância $$ d( x_1, x_2 ) = \sqrt{ x_1^2 + x_2^2 + \gamma( x_1, x_2 )^2 } $$ e plotaremos o ponto obtido acima.
# Implementamos a função distância:
def d( x ):
return ( x[ 0 ] ** 2 + x[ 1 ] ** 2 + gamma( x ) ** 2 ) ** 0.5
n_pts = 512
X = np.outer( np.linspace( -1.0 + c[ 0 ], 1.0 + c[ 0 ], n_pts ), np.ones( n_pts ) )
Y = np.outer( np.linspace( -1.0 + c[ 1 ], 1.0 + c[ 1 ], n_pts ), np.ones( n_pts ) ).T
Z = np.empty( ( n_pts, n_pts ) )
for i in range( n_pts ):
for j in range( n_pts ):
Z[ i, j ] = d( [ X[ i, j ], Y[ i, j ] ] )
fig = pp.figure( figsize = ( 12, 12 ) )
pp.imshow( np.flipud( Z ), extent = ( X[ 0, 0 ], X[ -1, 0 ], Y[ 0, 0 ], Y[ 0, -1 ] ) )
pp.colorbar()
pp.scatter( x[ -1 ][ 0 ], x[ -1 ][ 1 ], c = 'r' )
pp.show()
Exercício 1¶
Refaça o exemplo acima com diferentes valores para $\boldsymbol c$ e $l$. Em particular, mantendo os outros parâmetros mas fazendo $c_3 = 0$ as iterações de ponto fixo convergem.
Para este caso estude a velocidade de convergência obtida e compare com a norma da matriz jacobiana no ponto limite
Exercício 2¶
Refaça o exemplo acima com superfícies mais complicadas. Seja criativo na hora de escolher essas superfícies e aproveite o código acima para resolver o problema (basta reimplementar as funções da segunda célula de código).