33 Encontrando Raízes Unidimensionais

Esse capítulo descreve rotinas para encontrar raízes de funções arbitrárias unidimensionais. A bilbioteca fornece componentes de baixo nível para uma variedade de resolvedores de equações e testes de convergência. Esses resolvedores e testes podem ser combinados pelo usuário para alcançar a solução desejada, com acesso completo aos passos intermediários da iteração. Cada classe de métodos usa a mesma estrutura, de forma que você pode alternar entre resolvedores em tempo de execução sem precisar recompilar seu programa. Cada instância de um resolvedor mantém registros de seu próprio estado, permitindo aos resolvedores serem usados em programas com várias linhas de execução.

O arquivo de cabeçalho ‘gsl_roots.h’ conté protótipos para funções que encontram raízes e declarações relacionadas.

33.1 Visão Geral

Algoritmos de busca de raízes unidimensionais podem ser divididos em duas classes, delimitadores de raízes e melhoradores de raízes. Algoritmos que procedem a delimitação de uma raíz são garantidamente convergentes. Algoritmos de delimitação iniciam com uma região delimitada onde sabidamente contém uma raíz. O tamanho dessa região delimitada é reduzido, iterativamente, até que abrace a raíz com uma tolerância desejada. Esses algoritmos fornecem uma rigorosa estimativa de erro para a localização da raíz.

A técnica do melhorador de raíz tenta melhorar uma suposição inicial para a raíz. Esses algoritmos convergem somente se iniciados “perto o suficiente” da raíz, e sacrificam uma rigorosa delimitação de erro em troca de maior velocidade. Aproximando o comportamento de uma função nas vizinhanças de uma raíz os algoritmos de melhoramento de raízes tentam encontrar um aperfeiçoamento de ordem elevada de uma suposição inicial. Quando o comportamento da função é compatível como o algoritmo e uma boa suposição inicial está disponível um algoritmo de melhoramento pode fornecer uma rápida convergência.

Na GSL ambos os tipos de algoritmo estão disponíveis em estruturas semelhantes. O usuário fornece um norteador de alto nível para os algoritmos, e a biblioteca fornece as funções individuais necessárias para cada um dos passos. Existem três fases principais de iteração. Os passos são,

• inicializar o estado do resolvedor, s, para o algoritmo T
• atualizar o estado do resolvedor s usando o iterador T
• testa o estado s com relação à convergência, e repetir a iteração se necessário

O estado para resolvedores por delimitação é mantido em uma estrutura gsl_root_fsolver. O procedimento de atualização usa somente funções de avaliação (não derivadas). O estado de resolvedores por melhoria é mantido em uma estrutura gsl_root_fdfsolver. A atualização requer que ambas a função e sua derivada (daí o nome fdf) sejam fornecidas pelo usuário.

33.2 Ressalvas

Note que funções que encontram raízes podem somente buscar por uma raíz de cada vez. Quando houverem muitas raízes na área de busca, a primeira raíz a ser encontrada irá ser retornada; todavia é difícil prever qual das raízes essa irá ser encontrada. Na maioria dos casos, nenhum erro irá ser informado se você entar encontrar uma raíz e uma área onde existe mais de uma raíz.

Cuidado deve ser tomado quando uma função puder ter uma raíz múltipla (tal como f(x) = (x-x_0)^2 ou f(x) = (x-x_0)^3). Não é possível usar algoritmos de delimitação de raízes sobre raízes de multiplicidade par. Para esses algoritmos o intervalo inicial deve conter um zero cruzado, onde a função é negativa em uma extremidade do intervalo e positiva em outra extremidade. Raízes com multiplicidade par não cruzam zero, mas somente tocam o zero enstantâneamente. Algoritmos baseados em delimitação de raízes irão ainda trabalhar para raízes de multiplicidade ímpar (e.g. cubicas, de quinto grau, …). Algorítmos de melhoramento de raízes geralmente trabalham com alta multiplicidade de raízes, mas em uma razão reduzida de convergência. Nesses casos o algoritmos de Steffenson pode ser usado para acelerar a convergência de multiplas raízes.

Enquanto não é absolutamente requerido que f tenha uma raíz dentro da região de busca, funções que encontram raízes numéricas não devem ser usadas atabalhoadamente para verificar a existência de raízes. Existem melhores caminhos para fazer isso. Pelo fato de ser fácil criar situações onde funções que procuram raízes numéricas podem falhar, é uma péssima idéia iniciar uma rotina que encontra raízes em uma função que você não sabe muito sobre ela. Em geral é melhor examinar a função visualmente montando um gráfico antes de buscar por uma raíz.

33.3 Inicializando o Resolvedor

Function: gsl_root_fsolver * gsl_root_fsolver_alloc (const gsl_root_fsolver_type * T)

Essa função retorna um apontador para uma instância recentemente alocada de um resolvedor do tipo T. Por exemplo, o seguinte código cria uma instancia de resolvedor por bissecção,

const gsl_root_fsolver_type * T 
  = gsl_root_fsolver_bisection;
gsl_root_fsolver * s 
  = gsl_root_fsolver_alloc (T);

Se não houver memória suficiente para criar o resolvedor então a função retorna um apontador nulo e o controlador de erro é chamado com um código de erro de GSL_ENOMEM.

Function: gsl_root_fdfsolver * gsl_root_fdfsolver_alloc (const gsl_root_fdfsolver_type * T)

Essa função retorna um apontador para uma recentemente alocada instância de um resolvedor baseado em derivadas do tipo T. Por exemplo, o seguinte código cria uma instância de um resolvedor de Newton-Raphson,

const gsl_root_fdfsolver_type * T 
  = gsl_root_fdfsolver_newton;
gsl_root_fdfsolver * s 
  = gsl_root_fdfsolver_alloc (T);

Se não houver memória suficiente para criar o resolvedor então a função retorna um apontador nulo e o controlador de erro é chamado com um código de erro de GSL_ENOMEM.

Function: int gsl_root_fsolver_set (gsl_root_fsolver * s, gsl_function * f, double x_lower, double x_upper)

Essa função inicializa, ou reinicializa, um resolvedor existente s para usar a função f e o intervalo inicial de busca [x_lower, x_upper].

Function: int gsl_root_fdfsolver_set (gsl_root_fdfsolver * s, gsl_function_fdf * fdf, double root)

Essa função incializa, ou reinicializa, um resolvedor existente s para usar a função e a derivada fdf e a suposição inicial root.

Function: void gsl_root_fsolver_free (gsl_root_fsolver * s)

Function: void gsl_root_fdfsolver_free (gsl_root_fdfsolver * s)

Essas funções liberam toda memória associada ao resolvedor s.

Function: const char * gsl_root_fsolver_name (const gsl_root_fsolver * s)

Function: const char * gsl_root_fdfsolver_name (const gsl_root_fdfsolver * s)

Essas funções retornam um apontador para o nome do reolvedor. Por exemplo,

printf ("s é um resolvedor '%s'\n",
        gsl_root_fsolver_name (s));

pode imprimir alguma coisa como s é resolvedor 'bisection'.

33.4 Fornecendo a função a ser resolvida

Você deve fornecer uma função contínua de uma variável para a rotina que procura raízes, e, algumas vezes, sua primeira derivada. Com o objetivo de permitir parâmetros gerais as funções são definidas pelos seguintes tipos de dados:

Data Type: gsl_function

Esses tipos de dados define uma função geral com parâmetros.

double (* function) (double x, void * params)

essa função deve retornar o valor f(x,params) para argumentos x e parâmetros params

void * params

um apontador para os parâmetros da função

Aqui está um exemplo para a função quadrática geral,

struct my_f_params { double a; double b; double c; };

double
my_f (double x, void * p) {
   struct my_f_params * params 
     = (struct my_f_params *)p;
   double a = (params->a);
   double b = (params->b);
   double c = (params->c);

   return  (a * x + b) * x + c;
}

gsl_function F;
struct my_f_params params = { 3.0, 2.0, 1.0 };

F.function = &my_f;
F.params = &params;

A função f(x) pode ser avaliada usando a macro GSL_FN_EVAL(&F,x) definida em ‘gsl_math.h’.

Data Type: gsl_function_fdf

Esse tipo de dado define uma função geral com parâmetros e sua primeira derivada.

double (* f) (double x, void * params)

essa função deve retornar o valor de f(x,params) para argumento x e parâmetros params

double (* df) (double x, void * params)

essa função deve retornar o valor da derivada de f em relação a x, f'(x,params), para o argumento x e parâmetros params

void (* fdf) (double x, void * params, double * f, double * df)

essa função deve ajustar os valores da função f para f(x,params) e sua primeira derivada df para f'(x,params) para o argumento x e parâmetros params. Essa função fornece uma otimização das funções eparadas para f(x) e f'(x)—é sempre mais rápido calcular a função e sua derivada ao mesmo tempo.

void * params

um apontador para os parâmetros da função

Aqui está um exemplo onde f(x) = 2\exp(2x):

double
my_f (double x, void * params)
{
   return exp (2 * x);
}

double
my_df (double x, void * params)
{
   return 2 * exp (2 * x);
}

void
my_fdf (double x, void * params, 
        double * f, double * df)
{
   double t = exp (2 * x);

   *f = t;
   *df = 2 * t;   /* uses existing value */
}

gsl_function_fdf FDF;

FDF.f = &my_f;
FDF.df = &my_df;
FDF.fdf = &my_fdf;
FDF.params = 0;

A função f(x) pode ser avaliada usando a macro GSL_FN_FDF_EVAL_F(&FDF,x) e a derivada f'(x) pode ser avaliada usando a macro GSL_FN_FDF_EVAL_DF(&FDF,x). Ambas a função y = f(x) e sua derivada dy = f'(x) podem ser avaliados ao mesmo tempo usando a macro GSL_FN_FDF_EVAL_F_DF(&FDF,x,y,dy). A macro armazena f(x) no seu argumento y e f'(x) no seu argumento dy—ambas essas devem ser apontadores para double.

33.5 Busca de Limites e Suposições

Você fornece ou limites de busca ou uma suposição inicial; essa seção explana como limites de busca e suposisções de trabalho e como argumentos de funções as controlam.

Uma suposição é simplesmente um valor de x que é iterado até que esse x esteja dentro da precisão desejada de uma raíz. Essa raíz dentro de uma precisão toma a forma de um double.

Limites de busca são extremidas de um intervalo que é iterado até que o comprimento do intervalo seja menor que a precisão requisitada. O intervalo é definido por dois valores, o limite inferior e o limite superior. Se as extremidas são pensadas para serem inclídas no intervalo ou não depende do contexto no qual o intervalo é usado.

33.6 Iteração

A seguinte função norteiam a iteração de cada algoritmo. Cada função executa uma iteração para atualizar o estado de qualquer resolvedor do tipo correspondente. A mesma função trabalha para todos os resolvedores de forma que diferentes métodos podem ser substituídos durante a execução do programa sem modificações no código.

Function: int gsl_root_fsolver_iterate (gsl_root_fsolver * s)

Function: int gsl_root_fdfsolver_iterate (gsl_root_fdfsolver * s)

Essas funções executam uma iteração simples do resolvedor s. Se a iteração encontrar um problema inexperado então um código de erro irá ser retornado,

GSL_EBADFUNC

a iteração encontrou um ponto singular onde a função ou sua derivada avaliou para Inf ou NaN.

GSL_EZERODIV

a derivada da função tende a zero no ponto de, evitando que o algoritmo continue sem uma divisão por zero.

O resolvedor mantém uma melhor estimativa atual da raíz o tempo todo. Os resolvedores por delimitação também mantem registros dos melhores limites de intervalo atual da raíz. Essa informação pode ser acessada com as seguintes funções auxiliares,

Function: double gsl_root_fsolver_root (const gsl_root_fsolver * s)

Function: double gsl_root_fdfsolver_root (const gsl_root_fdfsolver * s)

Essas funções retornam a estimativa atual da raíz para o resolvedor s.

Function: double gsl_root_fsolver_x_lower (const gsl_root_fsolver * s)

Function: double gsl_root_fsolver_x_upper (const gsl_root_fsolver * s)

Essas funções retornam o atual intervalo delimitado para o resolvedor s.

33.7 Busca por Parâmetros de Parada

Um procedimento de localização de raízes deve parar quando uma das seguintes condições for verdadeira:

• Uma raíz tiver sido encontrada dentro da precisão especificada pelo usuário.
• Um número número máximo de iterações especificado pelo usuário tiver sido alcançado.
• Um erro ocorrer.

A manipulação dessas condições está sob o controle do usuário. As funções abaixo permitem ao usuário testar a precisão do resultado atual de muitas formas padronizadas.

Function: int gsl_root_test_interval (double x_lower, double x_upper, double epsabs, double epsrel)

Essa função testa a convergência do [x_lower, x_upper] com erro absoluto epsabs e erro relativo epsrel. O teste retorna GSL_SUCCESS se a seguinte condição for alcançada,

|a − b| < epsabs + epsrel   min (|a|,|b|)

quando o intervalo x = [a,b] não inclui a orígem. Se o intervalo inclui a orígem então \min(|a|,|b|) é substituído por zero (que o menor valor de |x| sobre o intervalo). Isso garante que o erro relativo é precisamente estimado para raízes próximas da orígem.

Essa condição sobre o intervalo também implica que qualquer estimativa da raíz r no intervalo satisfaz as mesmas condições com relação à verdadeira raíz r^*,

|r − r*| < epsabs + epsrel  r*

assumindo que a verdadeira raíz r^* esteja contida dentro do intervalo.

Function: int gsl_root_test_delta (double x1, double x0, double epsabs, double epsrel)

Essa função testa a convergência da sequência …, x0, x1 com erro absoluto epsabs e erro relativo epsrel. O teste retorna GSL_SUCCESS se a seguinte condição for alcançada,

|x1 − x0| < epsabs + epsrel  |x1|

e retorna GSL_CONTINUE de outra forma.

Function: int gsl_root_test_residual (double f, double epsabs)

Essa função testa o valor residual f contra o limite de erro absoluto epsabs. O teste retorna GSL_SUCCESS se a seguinte condição for alcançada,

|f| < epsabs

e retorna GSL_CONTINUE de outra forma. Esse critério é adequado para situações onde a localização precisa da raíz, x, é fornecida de forma secundária um valor pode ser encontrado onde o residual, |f(x)|, for pequeno o suficiente.

33.8 Algoritmos Delimitadores de Raízes

O algoritmo delimitador de raízes descrito nessa seção requer um intervalo inicial que garantidamente contém uma raíz—se a e b forem as extremidades do interval então f(a) deve diferir em sinal de f(b). Isso garante que a função cruza a raíz em ao menos uma voz no intervalo. Se um intervalo válido inicial é usado então esse algoritmo não pode falhar, contanto que a função seja bem comportada.

Note que um algoritmo de delimitação não encontra raízes de grau par, uma vez que esses não cruzam o eixo x.

Solver: gsl_root_fsolver_bisection

O algoritmo de bissecção é o método mais simples de delimitação de raízes de uma função. É também o algoritmo mais lento fornecido pela pela biblioteca, com convergência linear.

A cada iteração, o intervalo é bisectado e o valor da função no ponto médio é calculado. O sinal desse valor é usado para determinar qual metade do intervalo não contém uma raíz. Aquela metade é descartada para fornecer um novo intervalo, menor, contendo a raíz. Esse procedimento pode ser continuado indefinidamente até que o intervalo seja suficientemente pequeno.

Em qualquer tempo a estimativa atual da raíz é tomado como o ponto médio do intervalo.

Solver: gsl_root_fsolver_falsepos

O algoritmo de falsa posição é um método de encontrar raízes baseado em interpolação linear. Sua convergência é linear, mas é comumente mais rápido que a bissecção.

A cada iteração uma linha é desenhada entre as extremidades (a,f(a)) e (b,f(b)) e o ponto onde essa linha cruza o eixo x tomado como um “ponto médio”. O valor da função nesse ponto é calculado e seu sinal é usado para determinar qual lado do intervalo não contém uma raíz. Aquele lado é descartado para fornecer um novo intervalo, menor, contendo a raíz. Esse procedimento pode ser continuado indefinidamente até que o intervalo seja suficientemente pequeno.

A melhor estimativa da raíz é tomada a partir da interpolação linear do intervalo sobre a iteração atual.

Solver: gsl_root_fsolver_brent

O método de Brent-Dekker (65) (referido aqui como método de Brent) combina uma estratégia de interpolação com o algoritmo de bissecção. Essa combinação produz um rápido algoritmo que além disso é robusto.

A cada iteração o método de Brent aproxima a função usando uma curva de interpolação. Na primeira iteração a interpolação é linear entre dois pontos extremidades. Nas iterações subsequêntes o algoritmo usa um ajuste quadrado inverso para os últimos três pontos, para maior precisão. O ponto de intercecção da curva de interpolação com o eixo x é tomado como uma suposição para a raíz. Se esse ponto de intersecção encontra-se dentro dos limites do intervalo atual então ponto interpolado é aceito, e usado para gerar um intervalo menor. Se o ponto interpolado não for aceito então o algoritmo retorna para um passo comum de bissecção.

A melhor estimativa da raíz é tomada a partir da mais recente interpolação ou bissecção.

33.9 Algoritmos para Encontrar Raízes Derivativos

Os algoritmos de melhoria de raíz descritos nessa seção requerem uma suposição inicial para a localização da raíz. Não existe garantia absoluta de convergência—a função deve ser adequada para essa técnica e a suposição inicial deve estar suficientemente próxima à raíz para que o procedimento funcione. Quando essas condições forem satisfeitas então a convergência é quadrática.

Esses algoritmos fazem uso de ambos a função e sua derivada.

Derivative Solver: gsl_root_fdfsolver_newton

O Método de Newton é o algoritmo padrão de melhoria de raízes. O algoritmo inicia com uma suposição inicial para a localização da raíz. A cada iteração, uma linha tangente à função f é desenhada naquela posição. O ponto onde essa linha cruza o eixo x torna-se a nova suposição. A iteração é definida pela seguinte sequência,

xi+1 = xi − f(xi) f′(xi)

O método de Newton converge quadráticamente para raízes simples, e linearmente para raízes múltiplas.

Derivative Solver: gsl_root_fdfsolver_secant

O método da secante é uma versão simplificada do método de Newton que não requer o cálculo de derivadas em todos os passos.

Na primeira iteração o algoritmo começa com o método de Newton, usando a derivada para calcular um primeiro passo,

x1 = x0 − f(x0) f′(x0)

As iterações subsequêntes evitam a avaliação da derivada substituindo a derivada por uma estimativa numérica, a inclinação da linha que passa pelos dois pontos anteriores,

xi+1 = xi − f(xi) f′est  onde  f′est = f(xi) − f(xi−1) xi − xi−1

Quando a derivada não muda de forma significativa nas vizinhanças da raíz o método da secante é útil e econômico. Assintóticamente o método da secante é mais rápido que o método de Newton sempre que o custo de avaliação da derivada é mais que 0.44 vezes o custo da avaliação da própria função. Da mesma forma que todos os métodos de cálculo do valor numérico de uma derivada a estimativa pode ser vítima de erros de cancelamento se a separação dos pontos tornar-se muito pequena.

Para raízes simples, o método tem uma convergência de ordem (1 + \sqrt 5)/2 (aproximadamente 1.62). O método da secante converge linearmente para raízes múltiplas.

Derivative Solver: gsl_root_fdfsolver_steffenson

O Método de Steffenson(66) (67) fornece a mais rápida convergência de todas as rotinas. O método de Steffensen combina o algoritmos básico de Newton com uma aceleração “delta-quadrada” de Aitken. Se as iterações de Newton forem x_i então o procedimento de aceleração gera uma nova sequência R_i,

Ri = xi − (xi+1 − xi)2 (xi+2 − 2 xi+1 + xi)

que converge mais rapidamente que a sequência original sob condições razoáveis. A nova sequência requer três termos antes de poder produzir seu primeiro valor de forma que o método retorne valores acelerados a partir da segunda e subsequêntes iterações. Na primeira iteração o método retorna a estimativa comum de Newton. A iteração de Newton é também retornada se o denominador do termo de aceleração em alguma ocasião tornar-ser zero.

Da mesma forma que com todos os procedimentos de aceleração esse método pode tornar-se instável se a função não for bem comportada.

33.10 Exemplos

Para qualquer algoritmo de procura por raízes precisamos preparar a função que vai ser resolvida. Para o exemplo adiante iremos usar a equação quadrática geral descrita anteriormente. Primeiro precisamos de um arquivo de cabeçalho (‘demo_fn.h’) para definir os parâmetros da função,

struct quadratic_params
  {
    double a, b, c;
  };

double quadratic (double x, void *params);
double quadratic_deriv (double x, void *params);
void quadratic_fdf (double x, void *params, 
                    double *y, double *dy);

Colocamos as definições da função em um arquivo separado (‘demo_fn.c’),

double
quadratic (double x, void *params)
{
  struct quadratic_params *p 
    = (struct quadratic_params *) params;

  double a = p->a;
  double b = p->b;
  double c = p->c;

  return (a * x + b) * x + c;
}

double
quadratic_deriv (double x, void *params)
{
  struct quadratic_params *p 
    = (struct quadratic_params *) params;

  double a = p->a;
  double b = p->b;
  double c = p->c;

  return 2.0 * a * x + b;
}

void
quadratic_fdf (double x, void *params, 
               double *y, double *dy)
{
  struct quadratic_params *p 
    = (struct quadratic_params *) params;

  double a = p->a;
  double b = p->b;
  double c = p->c;

  *y = (a * x + b) * x + c;
  *dy = 2.0 * a * x + b;
}

O primeiro programa usa o resolvedor de funções gsl_root_fsolver_brent para o método de Brent e a quadrática geral acima para resolver a seguinte equação,

#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_roots.h>

#include "demo_fn.h"
#include "demo_fn.c"

int
main (void)
{
  int status;
  int iter = 0, max_iter = 100;
  const gsl_root_fsolver_type *T;
  gsl_root_fsolver *s;
  double r = 0, r_expected = sqrt (5.0);
  double x_lo = 0.0, x_hi = 5.0;
  gsl_function F;
  struct quadratic_params params = {1.0, 0.0, -5.0};

  F.function = &quadratic;
  F.params = &params;

  T = gsl_root_fsolver_brent;
  s = gsl_root_fsolver_alloc (T);
  gsl_root_fsolver_set (s, &F, x_lo, x_hi);

  printf ("using %s method\n", 
          gsl_root_fsolver_name (s));

  printf ("%5s [%9s, %9s] %9s %10s %9s\n",
          "iter", "lower", "upper", "root", 
          "err", "err(est)");

  do
    {
      iter++;
      status = gsl_root_fsolver_iterate (s);
      r = gsl_root_fsolver_root (s);
      x_lo = gsl_root_fsolver_x_lower (s);
      x_hi = gsl_root_fsolver_x_upper (s);
      status = gsl_root_test_interval (x_lo, x_hi,
                                       0, 0.001);

      if (status == GSL_SUCCESS)
        printf ("Converged:\n");

      printf ("%5d [%.7f, %.7f] %.7f %+.7f %.7f\n",
              iter, x_lo, x_hi,
              r, r - r_expected, 
              x_hi - x_lo);
    }
  while (status == GSL_CONTINUE && iter < max_iter);

  gsl_root_fsolver_free (s);

  return status;
}

Aqui está os resultados das iterações,

$ ./a.out 
using brent method
 iter [    lower,     upper]      root        err  err(est)
    1 [1.0000000, 5.0000000] 1.0000000 -1.2360680 4.0000000
    2 [1.0000000, 3.0000000] 3.0000000 +0.7639320 2.0000000
    3 [2.0000000, 3.0000000] 2.0000000 -0.2360680 1.0000000
    4 [2.2000000, 3.0000000] 2.2000000 -0.0360680 0.8000000
    5 [2.2000000, 2.2366300] 2.2366300 +0.0005621 0.0366300
Converged:                            
    6 [2.2360634, 2.2366300] 2.2360634 -0.0000046 0.0005666

Se o programa for modificado para usar o resolvedor por bissecção ao invés do método de Brent, mudando gsl_root_fsolver_brent para gsl_root_fsolver_bisection a lenta convergência do método de Bissecção pode ser observada,

$ ./a.out 
using bisection method
 iter [    lower,     upper]      root        err  err(est)
    1 [0.0000000, 2.5000000] 1.2500000 -0.9860680 2.5000000
    2 [1.2500000, 2.5000000] 1.8750000 -0.3610680 1.2500000
    3 [1.8750000, 2.5000000] 2.1875000 -0.0485680 0.6250000
    4 [2.1875000, 2.5000000] 2.3437500 +0.1076820 0.3125000
    5 [2.1875000, 2.3437500] 2.2656250 +0.0295570 0.1562500
    6 [2.1875000, 2.2656250] 2.2265625 -0.0095055 0.0781250
    7 [2.2265625, 2.2656250] 2.2460938 +0.0100258 0.0390625
    8 [2.2265625, 2.2460938] 2.2363281 +0.0002601 0.0195312
    9 [2.2265625, 2.2363281] 2.2314453 -0.0046227 0.0097656
   10 [2.2314453, 2.2363281] 2.2338867 -0.0021813 0.0048828
   11 [2.2338867, 2.2363281] 2.2351074 -0.0009606 0.0024414
Converged:                            
   12 [2.2351074, 2.2363281] 2.2357178 -0.0003502 0.0012207

O programa seguinte resolve a mesma função usando um resolvedor por derivadas no lugar da bissecção.

#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_roots.h>

#include "demo_fn.h"
#include "demo_fn.c"

int
main (void)
{
  int status;
  int iter = 0, max_iter = 100;
  const gsl_root_fdfsolver_type *T;
  gsl_root_fdfsolver *s;
  double x0, x = 5.0, r_expected = sqrt (5.0);
  gsl_function_fdf FDF;
  struct quadratic_params params = {1.0, 0.0, -5.0};

  FDF.f = &quadratic;
  FDF.df = &quadratic_deriv;
  FDF.fdf = &quadratic_fdf;
  FDF.params = &params;

  T = gsl_root_fdfsolver_newton;
  s = gsl_root_fdfsolver_alloc (T);
  gsl_root_fdfsolver_set (s, &FDF, x);

  printf ("using %s method\n", 
          gsl_root_fdfsolver_name (s));

  printf ("%-5s %10s %10s %10s\n",
          "iter", "root", "err", "err(est)");
  do
    {
      iter++;
      status = gsl_root_fdfsolver_iterate (s);
      x0 = x;
      x = gsl_root_fdfsolver_root (s);
      status = gsl_root_test_delta (x, x0, 0, 1e-3);

      if (status == GSL_SUCCESS)
        printf ("Converged:\n");

      printf ("%5d %10.7f %+10.7f %10.7f\n",
              iter, x, x - r_expected, x - x0);
    }
  while (status == GSL_CONTINUE && iter < max_iter);

  gsl_root_fdfsolver_free (s);
  return status;
}

Aqui está os resultados para o método de Newton,

$ ./a.out 
using newton method
iter        root        err   err(est)
    1  3.0000000 +0.7639320 -2.0000000
    2  2.3333333 +0.0972654 -0.6666667
    3  2.2380952 +0.0020273 -0.0952381
Converged:      
    4  2.2360689 +0.0000009 -0.0020263

Note que o erro pode ser estimado mais precisamente tomando a diferença entre a iteração atual e a seguinte em lugar de usar a iteração anterior. Os outros resolvedores derivativos podem ser investigados modificando gsl_root_fdfsolver_newton para gsl_root_fdfsolver_secant ou gsl_root_fdfsolver_steffenson.

33.11 Referências e Leituras Adicionais

Para mais informação sobre o algoritmos de Brent-Dekker veja os seguintes dois artigos,

• R. P. Brent, “An algorithm with guaranteed convergence for finding a zero of a function”, Computer Journal, 14 (1971) 422–425
• J. C. P. Bus and T. J. Dekker, “Two Efficient Algorithms with Guaranteed Convergence for Finding a Zero of a Function”, ACM Transactions of Mathematical Software, Vol. 1 No. 4 (1975) 330–345

Esse documento foi gerado em 23 de Julho de 2013 usando texi2html 5.0.