1.3.2 Raízes de Equações — Raízes de Polinômios — O Método de Muller
1. Introdução
O Método de Muller é uma técnica modificada do Método da Secante, mas que ao contrário dessa, não estima a raiz de uma função prolongando uma reta através de dois pontos — fazendo com que esta reta seja secante à curva da função –, e sim utiliza-se de uma parábola através de três pontos para aproximação da derivada.
A Figura 1 contém uma ilustração de uma aproximação qualquer, utilizando-se pontos em comum entre ambos os métodos.
Através de um processo iterativo, substitui-se os coeficientes para obter o ponto da parábola que intercepta o eixo abscisso.
2. Desenvolvimento do Método
Pelo fato do Método de Muller aproximar uma parábola, é conveniente descrever uma função de aproximação que tenha a forma quadrática, tal como (eq1)
Escolhendo três pontos interceptantes à função: , e para substituí-los na Equação 1 , geramos o conjunto de equações:
Para simplificar a função , tomaremos , o que nos leva a
obter o conjunto de equações:
A partir destas três equações, podemos determinar as três incógnitas: , e . Como dado pela última equação, , que é o valor da função avaliada na terceira aproximação. Ao sobstituirmos este coeficiente nas outras equações, obtemos:
Para simplificar esta demonstração e uma possível codificação, definimos as seguintes variáveis:
Com essas transformações, as Equações 8 e 9 podem ser escritas de uma forma mais simplificada:
isolando os coeficientes a, b e c,
Para encontrarmos a raiz, aplicamos a fórmula quadrática na próxima iteração. Como temos , e , podemos encontrar da seguinte maneira:
isolando ,
Em um programa, os passos acima são repetidos até quando a aproximação estiver com um erro tolerável. Sendo assim, em uma generalização, , , e e o erro relativo é calculado como:
Dos passos acima é possível ver que a equação que estima fornece duas raízes. No método de Muller, a estratégia é escolher o sinal de que majore o denominador, aumentando a velocidade de convergência.
Embora este algoritmo esteja sendo apresentado na busca de raízes de polinômios, ele pode ser utilizado para encontrar raízes de outras classes de funções, da mesma forma que o Método de Newton-Raphson.
</p>
3. Implementação
! -
! muller
! return the [complex] root of a function (using the Muller Method)
!
! F: function pointer
! x: some initial point
! imax: max of iterations allowed
! errto: tolerated error
!
! 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/>
! ago 2009
!
complex function muller(F, x, errto, imax)
implicit none;
complex, external :: F;
real, intent(in) :: errto;
real, intent(in) :: x;
integer, intent(in) :: imax;
! iteration counter !
integer :: iterCount;
! relative error !
real :: dx;
!
complex :: a,b,c;
complex :: x0,x1,x2,x3;
complex h0, h1, d0, d1, delta, denominator;
! initialization block
call random_seed();
call random_number(dx);
x2 = cmplx(x,0.);
x1 = cmplx(x + dx*x,0.);
x0 = cmplx(x - dx*x, 0.);
dx = 100.d0;
iterCount = 0;
do
iterCount = iterCount + 1;
h0 = x1 - x0;
h1 = x2 - x1;
d0 = (f(x1) - f(x2))/h0;
d1 = (f(x2) - f(x1))/h1;
a = (d1 - d0)/(h1+h0);
b = a*h1 + d1;
c = f(x2);
delta = sqrt(b*b - 4.*a*c);
if (abs(b+delta) > abs(b-delta)) then
denominator = b + delta;
else
denominator = b - delta;
end if
dx = -2.*c/denominator;
x3 = x2 + dx;
if( (abs(dx) < errto*real(x) ) .or. iterCount >imax )then
exit;
else
x0 = x1;
x1 = x2;
x2 = x3;
end if
end do
muller = x3;
return;
end function muller
Para inicializar x0 e x1 usamos uma variável aleatória para definir o intervalo em torno do ponto passado x. Isso serve para delegar à função a escolha da estimativa inicial dos pontos utilizados pelo Método de Muller.
Esta abordagem pode gerar um comportamento diferente em cada execução. Ou seja, para um mesmo x, em uma execução podem ser necessárias y iterações, enquanto que ao rodar o programa outra vez, a raiz seja convergida em z iterações.
Podemos evitar este comportamento substituindo as linhas
x1 = cmplx(x + dx*x, 0.);
x0 = cmplx(x – dx*x, 0.);
por
x1 = cmplx(x + s*x, 0.);
x0 = cmplx(x – s*x, 0.);
onde s é o passo que delimitará os outros pontos de avaliação. Veja que também não é necessário utilizar o mesmo s, pois qualquer ponto pode ser utilizado. Poderíamos utilizar sem problemas:
x1 = cmplx(x + s*x, 0.);
x0 = cmplx(x – k*x, 0.);
onde .
Este código pode ser encontrado em http://www.sawp.com.br/code/rootfind/muller.f90 </p>
4. 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/.