1.2.7 Raízes de Equações – Métodos Abertos – Interpolação Quadrática Inversa
1. Introdução
Interpolação Quadrática Inversa é uma técnica usada para obtermos zeros de equações não-lineares na forma . Ele é raramente implementado em aplicações como único método de busca de raízes, sendo usado como fator acelerador da convergência.
2. Desenvolvimento do Método
Seja a fórmula obtida a partir do Teorema de Expansão de Lagrange:
$$ x=y+\sum _{k=1}^{\infty }{\dfrac {a}^{k}{\dfrac {d^{d-1}{d{x}^{d-1}} \left( \left( g \left( x \right) \right) ^{k}\right) }{k!} $$
Onde . A partir desta expressão, podemos obter uma fórmula iterativa se tomarmos . Por questão de conveniência, expressamos na Equação 1 , e alteramos o índice para um intervalo finito tal que . Desta forma, a Equação 1 fica na forma:
$$ x_{i+1}=x_{i}+\sum _{j=1}^{m+1}{\dfrac { \left( -1 \right) ^{j}{\dfrac {\partial ^{j-1}{\partial {x}^{j-1}} \left( \left( g \left( x_{i}\right) \right) ^{j} \right) }{j!} $$
onde
$$ {\dfrac {\partial ^{j-1}{\partial {x^{j-1}}g \left( x_{j} \right)} = {\dfrac {\partial ^{j-1}{\partial {x^{j-1}} {\dfrac {1}{f(x_{i})}} $$
Desenvolvendo o somatório até , obtemos:
$$ x_{i+1}=x_{i}-{\dfrac {f \left( x_{i} \right) }{\dfrac {d}{dx}f \left( x_{i} \right) } – \dfrac {1}{2}\,{\dfrac { \left( f \left( x_{i}\right) \right) ^{2}{\dfrac {d^{2}{d{x}^{2}}f \left( x_{i}\right) }{ \left( {\dfrac {d}{dx}f \left( x_{i} \right) \right) ^{3}} $$
Aproximando a função por ela mesmo através da Interpolação de Lagrange por dois pontos e , temos:
$$ f \left( x \right) ={\frac { \left( x-x_{i} \right) f \left( x_{i-1} \right) }{x_{i-1}-x_{i}}+{\frac { \left( x-x_{i-1} \right) f \left( x_{i} \right) }{x_{i}-x_{i-1}} $$
Igualmente para a primeira derivada e segunda derivada:
$$ {\dfrac {d}{dx}f \left( x_{i} \right) ={\dfrac {f \left( x_{i-1} \right) -f \left( x_{i} \right) }{x_{i-1}-x_{i}} $$
$$ {\dfrac {d^{2}{d{x}^{2}}f \left( x_{i} \right) =-6\,{\dfrac {f \left( x_{i} \right) -f \left( x_{i-1} \right) }{ \left( x_{i} – x_{i-1} \right) ^{2}}+2\,{\dfrac {2\,{\dfrac {d}{dx}f \left( x_{i} \right) +{\dfrac {d}{dx}f \left( x_{i-1} \right) }{x_{i}-x_{i-1} }} $$
Levando as Expressões 5, 6 e 7 na Equação 2, obteremos:
$$ x_{i+1}={\frac {f \left( x_{i-1} \right) f \left( x_{i} \right) x_{i-2}}{ \left( f \left( x_{i-2} \right) -f \left( x_{i-1} \right) \right) \left( f \left( x_{i-2} \right) -f \left( x_{i} \right) \right) }+{\frac {f \left( x_{i-2} \right) f \left( x_{i} \right) x_{i-1}}{ \left( f \left( x_{i-1} \right) -f \left( x_{i-2} \right) \right) \left( f \left( x_{i-1} \right) -f \left( x _{i} \right) \right) }+{\frac {f \left( x_{i-2} \right) f \left( x_{i-1} \right) x_{i}}{ \left( f \left( x_{i} \right) -f \left( x_{i-2} \right) \right) \left( f \left( x_{i} \right) -f \left( x_{i-1} \right) \right) } $$
Para melhor melhorar a notação, tomemos . Desta forma, a fórmula 8 fica expressa como:
$$ x_{i+1}={\frac {y_{i-1}y_{i}x_{i-2}}{ \left( y_{i-2}-y_{i-1} \right) \left( y_{i-2}-y_{i} \right) }+{\frac {y_{i-2}y_{i}x_{i-1}}{ \left( y_{i-1}-y_{i-2} \right) \left( y_{i-1}-y_{i} \right) }+{\frac {y_{i-2}y_{i-1}x_{i}}{ \left( y_{i}-y_{i-2} \right) \left( y_{i}-y_{i-1} \right) }$$
3. Critérios usados na implementação
Nas implementações apresentadas a seguir, utilizamos o Método de Steffensen para obtermos a primeira aproximação e o Método da Secante para a segunda, para que de posse dessas, pudéssemos refinar a aproximação utilizando a Interpolação Quadrática Inversa.
Desta forma, na derivada:
$$ {\frac {d}{dx}f \left( x \right) ={\frac {f \left( x+h \right) -f \left( x \right) }{h}$$
aplicamos o critério de Steffensen — — na primeira aproximação. Ou seja,
$$ h=f \left( x_{i-2} \right) $$
o que nos permite obter a aproximação seguinte:
$$ x_{i-1}=x_{i-2}-{\frac {h}^{2}{f \left( x_{i-2}+h \right) -h} $$
Para próxima aproximação, aplicamos o Método da Secante:
$$x_{i}=x_{i-1}-{\frac {f \left( x_{i-1} \right) \left( x_{i-2} -x_{i-1} \right) }{f \left( x_{i-2} \right) -f \left( x_{i-1}\right) }$$
Por questão de conveniência, renomeamos as variáveis,
$$x_{0}=x_{i-2} $$
$$x_{1}=x_{i-1} $$
$$x_{2}=x_{i} $$
$$x_{3}=x_{i+1} $$
$$y_{0}=y_{i-2} $$
$$y_{1}=y_{i-1} $$
$$y_{2}=y_{i} $$
o que gera o conjunto de expressões utilizados no programa:
$$ x_{1}=x_{0}-{\frac {y_{0}}^{2}{f \left( x_{0}+y_{0}\right) -y_{0}} $$
$$ x_{2}=x_{1}-{\frac {y_{1} \left( x_{0}-x_{1} \right) }{y_{0}-y_{1}} $$
$$ x_{3}={\frac {y_{1}y_{2}x_{0}}{ \left( y_{0}-y_{1} \right) \left( y_{0}-y_{2} \right) }+{\frac {y_{0}y_{2}x_{1}}{ \left( y_{1}-y_{0} \right) \left( y_{1}-y_{2} \right) }+{\frac {y_{0}y_{1}x_{2}}{ \left( y_{2}-y_{0} \right) \left( y_{2}-y_{1} \right) } $$
</p>
4. Implementação
A primeira implementação apresentada a seguir consiste em aplicar os critérios apresentados na seção anterior a cada iteração, enquanto que a segunda utiliza tais critérios apenas para inicialização das variáveis, aplicando apenas a Fórmula 8 ao longo da execução.
Após algumas simulações, a primeira abordagem pareceu convergir mais rápido, desde que fosse dada uma boa aproximação inicial. Por outro lado, a implementação que utiliza somente a Interpolação Quadrática Inversa convergiu, independendo da aproximação inicial dada, ainda que sendo mais lenta.
Segue as implementações destas funções abaixo:
def IQI(f, x0=1.0, errto=0.001, imax=100):
"""
Return the root of a function (Inverse quadratic interpolation)
IQI(f, x0=1.0, errto=0.001, imax=100)
* f: Function where find roots
* x0: next point (estimative)
* errto: tolerated error ( for aprox. )
* imax: max of iterations allowed
return: root next x0
Author: Pedro Garcia [sawp@sawp.com.br]
see: http://www.sawp.com.br
License: Creative Commons
<http ://creativecommons.org/licenses/by-nc-nd/2.5/br/>
Jan 2009
"""
errno = errto + 1
iterCount = 0
while errno > errto and iterCount < imax:
y0 = f(x0)
if fabs(f(x0 + y0) - y0) != 0:
x1 = x0 - y0*y0/(f(x0 + y0) - y0)
else:
x3 = x0
break
y1 = f(x1)
if y0 != y1:
x2 = x1 - y1*(x0 - x1)/(y0 - y1)
else:
x3 = x1
break
# try Aitkens optimization
if (x2 - 2*x1 + x0) != 0:
alpha = (x0*x2 - x1*x1)/(x2 - 2*x1 + x0)
C1 = fabs((alpha - x1)/(alpha - x0))
C2 = fabs((alpha - x2)/(alpha - x1))
if (C1 < 1) and (C2 < 1):
x2 = alpha
y2 = f(x2)
delta = 7*[0]
delta[1] = y0 - y1
delta[2] = y0 - y2
delta[3] = y1 - y0
delta[4] = y1 - y2
delta[5] = y2 - y0
delta[6] = y2 - y1
def disturb(x):
if(x == 0):
return errto
else:
return x
delta = map(disturb, delta)
x3 = y1*y2*x0/((delta[1])*(delta[2])) + \
y0*y2*x1/((delta[3])*(delta[4])) + \
y0*y1*x2/((delta[5])*(delta[6]))
iterCount += 1
if x3 != 0:
errno = fabs((x3-x0)/x3)
x0 = x3
return x3
A segunda forma:
def IQI2(f, x0=1.0, errto=0.001, imax=100):
"""
Return the root of a function (Inverse quadratic interpolation)
IQI2(f, x0=1.0, errto=0.001, imax=100)
* f: Function where find roots
* x0: next point (estimative)
* errto: tolerate error ( for aprox. )
* imax: max of iterations allowed
return: root next x0
Author: Pedro Garcia [sawp@sawp.com.br]
see: http://www.sawp.com.br
License: Creative Commons
<http ://creativecommons.org/licenses/by-nc-nd/2.5/br/>
Jan 2009
"""
errno = errto + 1
iterCount = 0
y0 = f(x0)
if fabs(f(x0 + y0) - y0) != 0:
x1 = x0 - y0*y0/(f(x0 + y0) - y0)
else:
x1 = errto
y1 = f(x1)
if y0 != y1:
x2 = x1 - y1*(x0 - x1)/(y0 - y1)
else:
x2 = errto
while errno > errto and iterCount < imax:
y0 = f(x0)
y1 = f(x1)
y2 = f(x2)
delta = 7*[0]
delta[1] = y0 - y1
delta[2] = y0 - y2
delta[3] = y1 - y0
delta[4] = y1 - y2
delta[5] = y2 - y0
delta[6] = y2 - y1
def disturb(x):
if(x == 0):
return errto
else:
return x
delta = map(disturb, delta)
x3 = y1*y2*x0/((delta[1])*(delta[2])) + \
y0*y2*x1/((delta[3])*(delta[4])) + \
y0*y1*x2/((delta[5])*(delta[6]))
iterCount += 1
if x3 != 0:
errno = fabs((x3-x0)/x3)
x0 = x1
x1 = x2
x2 = x3
return x3
Nos arquivos http://www.sawp.com.br/code/rootfind/IQI.py e http://www.sawp.com.br/code/rootfind/IQI2.py apresento exemplos onde o uso do primeiro caso não converge, enquanto que o segundo
converge rapidamente. Nos mesmos arquivos, podemos ver outro exemplo onde a raiz é encontrada para ambas implementações. Nesse caso, podemos ver que o primeiro encontra o zero da função utilizando menos iterações. </p> </p>
5. Copyright
Este documento é disponível sob a licença Creative Commons. As regras dos direitos de cópia deste conteúdo estão acessíveis em http://creativecommons.org/licenses/by-nc-nd/2.5/br/.