O método de Newton-Raphson em várias dimensões consiste na seguinte iteração:
x(k+1) = x(k) - J(x(k))⁻¹f(x(k))
Onde J é o jacobiano. Para demonstrar a convergência deste método, basta mostrar que o mapa:
g(x) = x - J(x)⁻¹f(x)
é uma contração para valores próximos o suficiente de algum zero da função f e aplicar o Teorema do Ponto Fixo de Banach.
O seguinte código implementa o método e aplica ele ao exercício 2 da página 193 do livro Numerical Methods in Economics de Keneth Judd:
# -*- coding: utf-8 -*-
"""
Created on Fri Jul 8 02:28:33 2016
@author: Leonardo Marinho
Código que implementa o método de Newton-Raphson em duas variáveis e aplica-o
para o exercício 2 da página 193 do livro Numerical Methods in Economics, de
Keneth Judd. É a solução para a questão 8 da Aula 25 de Métodos Computacionais
em Economia - Aprendizagem de Máquinas: Aprendizagem por Reforço.
"""
import numpy as np
import matplotlib.pyplot as plt
from numpy.linalg import inv
from numpy.linalg import norm
from math import isnan
def jacobian_1(x):
"""Jacobiano do sistema 5.7.2 da pagina 175 do livro
Numerical Methods in Economics, de Keneth Judd."""
return np.array([[.2 * x[0] ** (-.8), .2 * x[1] ** (-.8)],
[.1 * x[0] ** (-.9), .4 * x[1] ** (-.6)]])
def jacobian_2(x):
"""Jacobiano do sistema 5.7.3 da pagina 175 do livro
Numerical Methods in Economics, de Keneth Judd."""
return np.array([[5. * (x[0] ** .2 + x[1] ** .2) ** 4. * .2 * x[0] ** (-.8),
5. * (x[0] ** .2 + x[1] ** .2) ** 4. * .2 * x[1] ** (-.8)],
[4. * (x[0] ** .1 + x[1] ** .4) ** 3. * .1 * x[0] ** (-.9),
4. * (x[0] ** .1 + x[1] ** .4) ** 3. * .4 * x[1] ** (-.6)]])
def function_1(x):
"""Sistema 5.7.2 da pagina 175 do livro
Numerical Methods in Economics, de Keneth Judd."""
return np.array([x[0] ** .2 + x[1] ** .2 - 2.,
x[0] ** .1 + x[1] ** .4 - 2.])
def function_2(x):
"""Sistema 5.7.2 da pagina 175 do livro
Numerical Methods in Economics, de Keneth Judd."""
return np.array([(x[0] ** .2 + x[1] ** .2) ** 5. - 32.,
(x[0] ** .1 + x[1] ** .4) ** 4. - 16.])
def newton(function, jacobian, initial, e=0.00001, d=0.00001, print_iter=True,
plot_iter = False):
"""Funcao que aplica o metodo de Newton-Raphson a funcao 'function' com
jacobiano 'jacobian' inseridos como funcoes do Python. 'initial' e o ponto
inicial para as iteracoes. 'e' e o nivel de precisao desejado e d e a
precisao desejada em relacao ao resultado. 'print_iter' escolhe se quer
que imprima os pontos intermediarios das iteracoes e 'plot_iter' escolhe
se quer imprimir estes pontos num grafico."""
x = initial
iter = 0
while True:
try:
if print_iter:
print x
if plot_iter:
plt.scatter(x[0], x[1])
s = np.dot(-inv(jacobian(x)), function(x)) # Equacao do metodo
x_ante = x
x = x + s
iter += 1
if norm(x_ante - x) <= e * (1 + norm(x)) or isnan(x[0]) or\
isnan(x[1]): # Verifica criterio de convergencia
break
except:
print "Erro"
break
if norm(function(x)) <= d: # Verifica se a solucao atende criterio
print "Convergiu em " + str(iter) + " iterações iniciando em " +\
str(np.round(initial, 2))
return x
else:
print "Método não convergiu iniciando em " + str(np.round(initial, 2))
# O loop abaixo roda o metodo no retangulo [0,2]**2 para function_1
for i in (np.array(range(20))+1.)/10:
for j in (np.array(range(20))+1.)/10:
newton(function_1, jacobian_1, [i, j], print_iter=False)
# O loop abaixo roda o metodo no retangulo [0,2]**2 para function_2
for i in (np.array(range(20))+1.)/10:
for j in (np.array(range(20))+1.)/10:
newton(function_2, jacobian_2, [i, j], print_iter=False)
###############################################################################
# Note que os loops acima mostram que a transformacao da function_2 tornou #
# a convergencia mais rapida, exigindo em media uma iteracao a menos no #
# dominio testado #
###############################################################################
# O comando abaixo imprime um grafico com os pontos da iteracao da function_1
newton(function_1, jacobian_1, [2, 2], plot_iter=True)
# O comando abaixo imprime um grafico com os pontos da iteracao da function_2
newton(function_2, jacobian_2, [2, 2], plot_iter=True)