1.1.1 Raízes de Equações — Métodos Intervalares — Método da Bissecção

1. Introdução (sec1)

Métodos intervalares que buscam raízes de funções consistem em checar se uma função muda de sinal dentro de um intervalo específico. Conhecidas também como “bracketing methods”, elas exigem duas estimativas iniciais para a raiz. Estas “estimativas” são os limites do intervalo inicial, contenedor da raiz. A partir do intervalo inicial, delimitado pelas estimativas, utilizam-se diferentes métodos para diminuir a largura deste intervalo a fim de aproximar a solução correta.

2. Métodos de buscas incrementais (sec2)

Para entendermos os métodos numéricos que encontram zeros de funções devemos compreender antes o conceito de busca incremental.

Métodos de buscas incrementais são métodos intervalares que detectam subintervalos no qual a função muda de sinal. Como definido na introdução como “bracking methods”, devemos alterar o intervalo de forma refinar a busca pela raiz.

A busca incremental divide o intervalo formado inicialmente em diversos subintervalos, depois busca nestes onde ocorre a mudança de sinal. O processo é repetido no subintervalo onde a raiz foi detectada, dividindo o subintervalo em pedaços menores, a fim de refinar a aproximação da raiz.

Um problema da potencial para a escolha da busca incremental é a escolha do comprimento do incremento. Se o comprimento é muito pequeno, a busca pode demorar muito tempo. Se o comprimento for muito grande, há possibilidade de o próprio incremento incorporar várias raízes, se elas estiverem muito próximas, não sendo possível detectá-las.

3. Método da bissecção (Método de Bolzano) (sec3)

O método da bissecção pode ser chamado de “método da divisão do intervalo na metade” e tem este nome porque é uma busca incremental no qual sempre se divide o intervalo ao meio, gerando dois novos subintervalos, contidos no anterior.

Após cada subdivisão, aplicam-se os critérios que caracterizam este método: se a função muda de sinal em um intervalo, calcula-se o seu ponto médio. A posição da raiz é determinada como sendo este ponto médio dentro do subintervalo no qual a mudança de sinal ocorreu, repetindo este processo até que a aproximação esteja próxima o suficiente.

O algoritmo, basicamente, pode ser descrito da seguinte forma:

1. Escolhemos os limites do intervalo que incorpora a raiz, chamando estes limites de $x\_{l}$ e $x\_{r}$ . Como uma raiz significa um ponto onde a função muda de sinal, o intervalo pode ser verificado como sendo contenedor da raiz se $F(x\_{r})F(x\_{l}) < 0$.
2. Uma estimativa da raiz pode ser aproximada como $x = \dfrac{ x_{l} + x{r} }{2}$ . Com isso, divide-se o intervalo inicial em dois subintervalos idênticos.
3. Verificamos para saber em qual subintervalo a raiz está:
   1. se $F(x)F(x\_{l}) < 0$ , a raiz está no intervalo da esquerda. Faça $x\_{r} = x$ e volte ao passo $2$ .
   2. se $F(x)F(x\_{l}) > 0$ , a raiz não esta no intervalo da esquerda, logo, está no da direita, então faça $x\_{l} = x$ e volte ao passo $2$ .
   3. se $F(x)F(x_{l}) = 0$ , a raiz foi encontrada e é $x$ .

4. Implementação (sec4)

def Bissecc(F, xl, xr, imax=500, errto=0.01):
    """
    Return the root of a function using Bissection Method

    Bissecc(fun, xl, xr, imax, errto)

    * fun: Function where find the roots
    * xl: left bound of interval
    * xr: right bound of interval
    * imax: max of iterations allowed
    * errto: tolerated error 

    Author: Pedro Garcia [sawp@sawp.com.br]
    License: Creative Commons
         <http ://creativecommons.org/licenses/by-nc-nd/2.5/br/>
    """

    x = 0
    iterCount = 0
    errno = 1
    while errno > errto and iterCount < imax:
        xold = x
        x = (xl + xr) / 2
        iterCount += 1
        if x != 0:
            errno = fabs((x - xold) / x) * 100
        test = F(xl) * F(x)
        if test < 0:
            xr = x
        elif test > 0:
            xl = x
        else:
            # test == 0 is when the F(x) is the root
           errno = 0
    return x

Esta implementação funciona bem para qualquer função F que retorne um valor real, sendo que outras variantes destes códigos podem ser implementadas, podendo usar tipos de tamanho maior, recursos de paralelização, etc. Apesar do código acima funcionar bem, ainda podemos melhorá-lo com algumas otimizações.

4.1. Otimizando o Método da Bissecção

Ainda que o código anterior permita bons resultados, em um programa real pode acontecer da função F precisar ser chamada milhares de vezes. Caso isso ocorra, é extremamente interessante buscar otimizações no método a fim de fornecer um algoritmo menos caro para ser processado.

No caso de um método numérico, como o da Bissecção, uma possível melhoria de se implementação consiste em minimizar o overhead das chamada de funções, como é feito no código abaixo. Ao comparamos ele com a implementação anterior, é possível observar que algumas variáveis foram inseridas a fim de evitar que uma mesma função seja chamada mais de uma vez para calcular o mesmo valor. Por exemplo:

def Bissecc_O1(F, xl, xr, imax=500, errto=0.01):
    """
    Return the root of a function using Bissection Method

    Bissecc_O1(fun, xl, xr, imax, errto)

    * fun: Function where find the roots
    * xl: left bound of interval
    * xr: right bound of interval
    * imax: max of iterations allowed
    * errto: tolerated error ( for aprox. ), in percents

    Author: Pedro Garcia [sawp@sawp.com.br]
    License: Creative Commons
         <http ://creativecommons.org/licenses/by-nc-nd/2.5/br/>
    """

    x = 0
    iterCount = 0
    errno = 1
    fl = F(xl) # used in O1
    while errno > errto and iterCount < imax:
        xold = x
        x = (xl + xr) / 2
        iterCount += 1
        fx = F(x) # O1
        if x != 0:
            errno = fabs((x - xold) / x) * 100
        test = fl * fx # O1 reduce the funcion call:test = F(xl)*F(x)
        if test < 0:
            xr = x
        elif test > 0:
            xl = x
        else:
            # test == 0 is when the F(x) is the root
           errno = 0
    return x

As linhas inseridas em relação ao código original possuem o comentário “#O1” , indicando a primeira otimização feita.

Sem as novas linhas inseridas, o código não minimiza a chamada de funções. Enquanto que no caso não otimizado, para cada iteração há duas chamadas de função que calculam o mesmo valor, no caso otimizado esta repetição de cálculo é desfeita quando acumulamos o valor da função anterior na variável “fl” . Assim, reduzimos a quantidade de chamadas da função de $2n$ vezes para apenas $(n+1)$, onde $n$ é o número de iterações necessárias para convergir à raiz.

Em casos onde $n$ seja muito grande, reduzir à metade o número de chamadas de função pode significar uma melhoria considerável para eficiência de todo sistema.

Outro passo para uma otimização interessante consiste em aplicar alguma espécie de otimização matemática — reduzindo-se a complexidade do algoritmo, aplicando-se otimização numérica, estocástica etc. Sempre que for aplicável uma otimização numérica, ela deve ser feita.

Uma propriedade útil do Método da Bissecção é que ele permite que o número de iterações necessárias para se chegar ao erro absoluto — o erro de tolerância que é passado na implementação, usado como critério de parada — sejam calculadas antecipadamente, ou seja, antes de começar as iterações.

Desta forma, antes de começar o algoritmo, o erro absoluto $E$ pode ser pré-determinado da seguinte maneira:

(eq1)

$$ E^0 = x_{r}^0 – x_{l}^0 $$

chamando o intervalo inicial delimitado por $\left[ x\_{l}, x\_{r} \right]$ de

(eq2)

$$ \Delta x^0 = x_{r}^0 – x_{l}^0 $$

temos que, o erro absoluto no início do programa será

(eq3)

$$ E^0 = \Delta x^0 $$

Portanto, este erro será dado como o tamanho do intervalo, que será sempre reduzido à metade a cada iteração. Sendo assim, logo na primeira iteração o erro terá a forma:

(eq4)

$$ E^0 = \dfrac{1}{2} \Delta x^0 $$

como a cada iteração, este processo é repetido, o n-ésimo erro absoluto será sempre a metade do erro anterior, sendo possível relacioná-lo com o primeiro intervalo dado:

(eq5)

$$ E^n = \dfrac{ \Delta x^0 }{ 2^n } $$

ou seja, se o erro de tolerância necessita de $n$ iterações, então ele será o n-ésimo erro dado pela fórmula acima.

Chamaremos de $E_{wanted}$ o erro máximo desejado, que é definido de acordo com as necessidades do programa em que a implementação está envolvida. A partir dele podemos prever o número de iterações máximas que serão necessárias para o método convergir satisfatoriamente à raiz.

(eq6)

$$ E_{wanted} = \dfrac{ \Delta x^0 }{ 2^n } $$
$$\Downarrow $$

(eq7)

$$ 2^n = \dfrac{ \Delta x^0 }{ E_{wanted} } $$
$$\Downarrow $$

(eq8)

$$ ln(2^n) = ln \left( \dfrac{ \Delta x^0 }{ E_{wanted} } \right) $$
$$\Downarrow $$

(eq9)

$$ n~ln(2) = ln \left( \dfrac{ \Delta x^0 }{ E_{wanted} } \right) = ln(\Delta x^0) – ln(E_{wanted}) $$

logo, o número de iterações máximas necessárias será:

(eq10)

$$ n = \dfrac{ ln \left( \dfrac{\Delta x^0}{ E_{wanted} } \right) }{ ln(2) } $$

Em um código real, ao aplicarmos este recurso, podemos reduzir o número de parâmetros passados, fornecendo à função apenas o erro de tolerância desejado. Um exemplo desta abordagem implementada pode ser vista abaixo:

def Bissecc_O2(F, xl, xr, errto=0.01):
    """
    Return the root of a function using Bissection Method

    Bissecc_O2(fun, xl, xr, imax, errto)

    * fun: Function where find the roots
    * xl: left bound of interval
    * xr: right bound of interval
    * imax: max of iterations allowed
    * errto: tolerated error ( for aprox. ), in percents

    Author: Pedro Garcia [sawp@sawp.com.br]
    License: Creative Commons
         <http ://creativecommons.org/licenses/by-nc-nd/2.5/br/>
    """

    x = 0
    iterCount = 0
    errno = 1
    fl = F(xl) # used in O1
    imax = int(ceil(fabs(log(fabs(xr-xl)/errto)/log(2)))) # used in O2
    while errno > errto and iterCount < imax:
        xold = x
        x = (xl+xr)/2
        iterCount += 1
        fx = F(x) # O1
        if x != 0:
            errno = fabs((x - xold) / x) * 100
        test = fl * fx # O1 reduce the funcion call:test = F(xl)*F(x)
        if test < 0:
            xr = x
        elif test > 0:
            xl = x
        else:
            # test == 0 implies x is the root
            errno = 0
    return x

Como o objetivo destas modificações é otimizar o programa, retiramos os testes baseados na convergência, visto que eles são desnecessários, pois o Método da Bissecção sempre converge e o número de iterações máximo é previsto. Desta forma, a variável errno se torna desnecessária, e sua retirada reduz o número de desvios condicionais que o programa deve fazer. Desta forma, o código mais otimizado possível é:

def Bissecc_O3(F, xl, xr, errto=0.01):
    """
    Return the root of a function using Bissection Method

    Bissecc_O3(fun, xl, xr, imax, errto)

    * fun: Function where find the roots
    * xl: left bound of interval
    * xr: right bound of interval
    * imax: max of iterations allowed
    * errto: tolerated error

    Author: Pedro Garcia [sawp@sawp.com.br]
    License: Creative Commons
         <http ://creativecommons.org/licenses/by-nc-nd/2.5/br/>
    """

    x = 0
    iterCount = 0
    fl = F(xl) # used in O1
    imax = int(ceil(fabs(log(fabs(xr-xl)/errto)/log(2)))) # used in O2
    while iterCount < imax:
        xold = x
        x = (xl + xr) / 2
        iterCount += 1
        fx = F(x) # O1
        if x != 0:
            errno = fabs((x - xold) / x)
        test = fl * fx # O1 reduce the funcion call:test = F(xl)*F(x)
        if test < 0:
            xr = x
        elif test > 0:
            xl = x
        else:
            # test == 0 is when the F(x) is the root
           errno = 0
    return x

Estes códigos estão disponíveis para download em:

• http://www.sawp.com.br/code/rootfind/bisseccao.py
• http://www.sawp.com.br/code/rootfind/bisseccao_optimized1.py
• http://www.sawp.com.br/code/rootfind/bisseccao_optimized2.py
• http://www.sawp.com.br/code/rootfind/bisseccao_optimized3.py

5. Copyright

Este documento está 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] Ralston, Anthony and Rabinowitz, Philip, A First Course in Numerical Analysis.,(2nd ed.) McGraw-Hill and Dover, (2001), 40-42.</p>

</a>
[2] N.B Franco, Cálculo Numérico, Pearson Prentice Hall, (2006), 66.

</fieldset>