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 \(f(x) = 0 \) . 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:
(eq1)


\( 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 \(x = x(y) \) . A partir desta expressão, podemos obter uma fórmula iterativa se tomarmos \(y=f \left( x \right) \) . Por questão de conveniência, expressamos na Equação 1 \(y_{{i}} = f_{i} = f \left( x_{{i}} \right) \), \(a = -1 \) e alteramos o índice \(k\) para um intervalo finito tal que \( k = 1 \ldots m+1 \) . Desta forma, a Equação 1 fica na forma:

(eq2)


\( 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
(eq3)

\( {\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é \(m = 1 \) , obtemos:
(eq4)


\( 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 \(x_{i} \) e \(x_{i-1} \), temos:
(eq5)


\( 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:
(eq6)

\( {\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}}}} \)

(eq7)


\( {\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:
(eq8)


\( 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 \(y_{k} = f(x_{k}) \) . Desta forma, a fórmula 8 fica expressa como:

(eq9)


\( 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:
(eq10)


\( {\frac {d}{dx}}f \left( x \right) ={\frac {f \left( x+h \right) -f \left( x \right) }{h}}\)

aplicamos o critério de Steffensen — \(h=f \left( x \right) \) — na primeira aproximação. Ou seja,
(eq11)

\( 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:
(eq12)


\(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:
(eq13)

\( x_{{1}}=x_{{0}}-{\frac {{y_{{0}}}^{2}}{f \left( x_{{0}}+y_{{0}}\right) -y_{{0}}}} \)

(eq14)


\( x_{{2}}=x_{{1}}-{\frac {y_{{1}} \left( x_{{0}}-x_{{1}} \right) }{y_{{0}}-y_{{1}}}} \)

(eq15)


\( 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) }} \)

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.

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/.

References


[1] http://en.wikipedia.org/wiki/Inverse_quadratic_interpolation
[2] Ralston, Anthony and Rabinowitz, Philip, A First Course in Numerical Analysis.,(2nd ed.) McGraw-Hill and Dover, (2001), 358–359.