Newton-Raphson em uma dimensão

November 27, 2011
By

This post was kindly contributed by Ridículas - go there to comment and to read the full post.

Trajetória obtida com os valores consecutivos nas iterações do método Newton-Raphson para estimação da raíz de uma função.

O método Newton-Raphson é empregado para encontrar o zero de funções. Em problemas de otimização, é empregado para encontrar o zero da função derivada f'(x) de nossa função objetivo f(x). O procedimento se baseia na aproximação em série de Taylor de primeira ordem ao redor de um valor arbitrário x_0. Considerando que desejamos obter o zero de uma função f(x), temos a expansão ao redor de x_0

f(x) \approx f(x_0) + f'(x_0)(x-x_0)

como queremos saber qual o valor de x para o qual f(x)=0, nós isolamos x na expressão acima, assim

x = x_0 + f(x_0)/f'(x_0)

e como esse valor é baseado em uma aproximação linear ao redor de x_0, agora atualizamos x_0 com o valor encontrado e repetimos esse procedimento até alcançar convergência por algum critério. Normalmente esses critérios são sobre as diferenças entre os valores obtidos em passos consecutivos.

No CMR abaixo eu apresento a obtenção raízes de uma função composta por \sin(x) e x^2. Nos escrevemos as funções, obtemos a derivada e aplicamos o método Newton-Raphson observando graficamente o histórico de iterações. No final eu coloquei um exemplo de como obter o ponto de máximo de uma função aplicando o método Newton-Raphson. Essas funções são didáticas e o método Newton-Raphson está disponível na função micEcon::maxNR() para maximização baseada em Newton-Raphson. No pacote animation a função newton.method() ilustra o funcionamento do método. Mais sobre métodos de otimização estão disponíveis na Task View Optimization do R. Em um post futuro vou ilustrar o uso do método Newton-Raphson em estimação por verossimilhança. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# newton-raphson em uma dimensão
# http://condicaoinicial.com/2011/08/sistema-nao-lineares-newton-raphson.html
# http://www.mat.ufrgs.br/~guidi/grad/MAT01169/laminas.calculo_numerico.5.pdf
# http://www.editora.ufla.br/site/_adm/upload/revista/23-3-1999_25.pdf
# http://www.dpi.inpe.br/~miguel/cursoemilioribas/MaterialRef/Curso1Estatistica_Corina.pdf

f <- quote(sin(x)-x^2/16)        # expressão da função, queremos obter a raíz
fx0 <- function(x){ eval(f) }    # função
curve(fx0, -10, 10); abline(h=0) # gráfico da função

f1 <- D(f,"x")                   # expressão da derivada
fx1 <- function(x){ eval(f1) }   # função
par(new=TRUE); curve(fx1, -10, 10, col=2, axes=FALSE) # gráfico

#-----------------------------------------------------------------------------
# método iterativo de Newton-Raphson

i <- 1       # valor inicial para o passo
dif <- 100   # valor inical para a diferença entre duas avaliações
tol <- 10^-9 # diferência mínima entre duas avaliações (tolerância)
dif <- 100   # número máximo de passos
x <- -7      # valor inicial para a raiz

while(i<100 & dif>tol){
  x[i+1] <- x[i]-fx0(x[i])/fx1(x[i])
  dif <- abs(x[i+1]-x[i])
  i <- i+1
  print(x[i])
}

#-----------------------------------------------------------------------------
# gráfico que ilustra a trajetória até a obtenção da raíz

curve(fx0, -10, 10)
for(j in 2:i){
  arrows(x[j-1], fx0(x[j-1]), x[j], fx0(x[j]), length=0.1, col=3)
  Sys.sleep(0.5)
}
abline(v=x[i], h=fx0(x[i]), col=2)

# experimente outros valores iniciais
# veja que essa função tem 2 raízes, cada valor inicial vai levar
# à apenas uma delas!

#-----------------------------------------------------------------------------
# fazendo animação gif

dir.create("seqpngs")
setwd("seqpngs")

png(file="p%02d.png", width=400, height=300)
for(j in 2:i){
  curve(fx0, -10, 10, main=paste("passo ", j-1, ", (x = ",
                        round(x[j],2), ")", sep=""))
  arrows(x[j-1], fx0(x[j-1]), x[j], fx0(x[j]), length=0.1, col=3, lwd=2)
  abline(v=x[j], h=fx0(x[j]), col="gray50")
}
abline(v=x[i], h=fx0(x[i]), col=2, lwd=2)
text(0, -3, label="CONVERGIU!", cex=2)
dev.off()

system("convert -delay 80 *.png example_1.gif")

#-----------------------------------------------------------------------------
# encontrando o ponto de máximo de uma função com newton-raphson

f <- quote(x*(2+0.5*x)^(-4))     # expressão da função, queremos obter MÁXIMO
fx0 <- function(x){ eval(f) }    # função
curve(fx0, 0, 10)

f1 <- D(f,"x")                   # expressão da derivada, obteremos o zero
fx1 <- function(x){ eval(f1) }   # função
curve(fx1, 0, 10, col=2)         # gráfico

f2 <- D(f1, "x")
fx2 <- function(x){ eval(f2) }   # função

#-----------------------------------------------------------------------------
# aplicando o newton-raphson

i <- 1       # valor inicial para o passo
dif <- 100   # valor inical para a diferença entre duas avaliações
tol <- 10^-9 # diferência mínima entre duas avaliações (tolerância)
dif <- 100   # número máximo de passos
x <- 0       # valor inicial para a raiz

while(i<100 & dif>tol){
  x[i+1] <- x[i]-fx1(x[i])/fx2(x[i])
  dif <- abs(x[i+1]-x[i])
  i <- i+1
  print(x[i])
}

#-----------------------------------------------------------------------------
# trajetória

curve(fx0, 0, 10)
for(j in 2:i){
  arrows(x[j-1], fx0(x[j-1]), x[j], fx0(x[j]), length=0.1, col=3)
  Sys.sleep(0.5)
}
abline(v=x[i], h=fx0(x[i]), col=2)

# experimente outros valores iniciais, como por exemplo 3 e 8
# para essa última função podemos obter o ponto de máximo analiticamente
#-----------------------------------------------------------------------------


Tags: , , ,

Comments are closed.

Contributors