6 Polinômios

Esse capítulo descreve funções para avaliar e resolver polinômios. Existem rotinas para encontrar raízes reais e complexas de equações quadráticas e cúbicas usando métodos analíticos. Um resolvedor polinomial analítico está também disponível para encontrar raízes de polinômios em geral com coeficientes reais (de qualquer ordem). As funções são declaradas no arquivo de cabeçalho ‘gsl_poly.h’.

6.1 Avaliação de polinômios

As funções descritas aqui avaliam o polinômio P(x) = c[0] + c[1] x + c[2] x^2 + \dots + c[len-1] x^{len-1} usando o método de Horner para estabilidade. Versões modificadas dessa função são usada quando HAVE_INLINE for definida.

Function: double gsl_poly_eval (const double c[], const int len, const double x): Essa função avalia um polinômio com coeficientes reais para a variável real x.

Function: gsl_complex gsl_poly_complex_eval (const double c[], const int len, const gsl_complex z): Essa função avalia um polinômio com coeficientes reais para a variável complexa z.

Function: gsl_complex gsl_complex_poly_complex_eval (const gsl_complex c[], const int len, const gsl_complex z): Essa função avalia um polinômio com coeficientes complexos para a variável complexa z.

Function: int gsl_poly_eval_derivs (const double c[], const size_t lenc, const double x, double res[], const size_t lenres): Essa função avalia um polinômio e suas derivadas armazenando o resultado no vetor estático res de tamanho lenres. O vetor estático de saída contém os valores de d^k P/d x^k para o valor especificado de x iniciando com k = 0.

6.2 Representação da Diferença Dividida de Polinômios

As funções descritas aqui tratam polinômios armazenados na representação de diferenças divididas de Newton. O uso de diferenças divididas é descrito em Abramowitz & Stegun seções 25.1.4 e 25.2.26 e em Burden e Faires, capítulo 3, e discutido brevemente abaixo.

Dado uma função f(x), um polinômio de interpolação de n-ésimo grau P_{n}(x) pode ser construído o qual concorda com f em n+1 pontos distintos x_0,x_1,...,x_{n}. Esse polinômio pode ser escrito de uma forma conhecida como representação de diferença dividida de Newton: onde as diferenças divididas [x_0,x_1,...,x_k] são definidas na seção 25.1.4 de Abramowitz e Stegun. Adicionalmente, é também possível construir um polinômio de interpolação de grau 2n+1 que também coincide com a primeira derivada de f nos pontos x_0,x_1,...,x_n. É chamado de polinômio de interpolação de Hermite e é definido como onde os elementos de z = \{x_0,x_0,x_1,x_1,...,x_n,x_n\} são definidos por z_{2k} = z_{2k+1} = x_k. As diferenças divididas [z_0,z_1,...,z_k] são discutidas em Burden e Faires, seção 3.4.

Function: int gsl_poly_dd_init (double dd[], const double xa[], const double ya[], size_t size): Essa função calcula uma representação de diferenças divididas do polinômio interpolado para os pontos (x, y) armazenados nos vetores estáticos xa e ya de comprimento size. Na saída as diferenças divididas de (xa,ya) são armazenados no vetor dd, também de comprimento size. Usando a notação acima, dd[k] = [x_0,x_1,...,x_k].

Function: double gsl_poly_dd_eval (const double dd[], const double xa[], const size_t size, const double x): Essa função avalia o polinômio armazenado na forma de diferença dividida nos vetores estáticos dd e xa de comprimento size no ponto x. Uma versão modificada dessa função é usada quando HAVE_INLINE for definida.

Function: int gsl_poly_dd_taylor (double c[], double xp, const double dd[], const double xa[], size_t size, double w[]): Essa função converte a representação de diferença dividida de um polinômio para uma expansão de Taylor. A representação de diferença dividida é fornecida nos vetores estáticos dd e xa de comprimento size. Na saída os coeficientes de Taylor do polinômio expandido sobre o ponto xp são armazenados no vetor estático c também de comprimento size. Um espaço de trabalho de comprimento size deve ser fornecido no vetor estático w.

Function: int gsl_poly_dd_hermite_init (double dd[], double za[], const double xa[], const double ya[], const double dya[], const size_t size): This function computes a divided-difference representation of the interpolating Hermite polynomial for the points (x, y) stored in the arrays xa and ya of length size. Hermite interpolation constructs polynomials which also match first derivatives dy/dx which are provided in the array dya also of length size. The first derivatives can be incorported into the usual divided-difference algorithm by forming a new dataset z = \{x_0,x_0,x_1,x_1,...\}, which is stored in the array za of length 2*size on output. On output the divided-differences of the Hermite representation are stored in the array dd, also of length 2*size. Using the notation above, dd[k] = [z_0,z_1,...,z_k].

6.3 Equações Quadráticas

Function: int gsl_poly_solve_quadratic (double a, double b, double c, double * x0, double * x1)

Essa função encontra as raízes reais da equação quadrática,

a x² + b x + c = 0

O número de raízes reais (ou zero, ou uma ou duas) é retornado, e a localização das raízes estão armazenadas em x0 e x1. Se nenhuma raíz real for encontrada então x0 e x1 não são modificadas. Se uma raíz real for encontrada (i.e. se a=0) então ela é armazenada em x0. Quando duas raízes reais são encontradas elas são armazenadas em x0 e x1 em ordem ascendente. O caso de raízes iguais não é considerado especial. Por exemplo (x-1)^2=0 irá ter duas raízes, que possuem valores exatamente iguais.

O número de raízes encontradas depende do sinal do discriminante b^2 - 4 a c. O discriminante estará sujeito a erros de arredondamento e cancelamento quando calculado em precisão dupla, e estará também sujeito a erros se os coeficientes do polinômios forem inexatos. Esses erros podem causar uma discreta modificação no número de raízes. Todavia, para polinômios com pequenos coeficientes inteiro o discriminante pode sempre ser calculado com grande exatidão.

Function: int gsl_poly_complex_solve_quadratic (double a, double b, double c, gsl_complex * z0, gsl_complex * z1)

Essa função encontra as raízes complexas da equação quadrática,

a z² + b z + c = 0

O número de raízes complexas é retornado (ou uma ou duas) e a localização das raízes estão armazenadas em z0 e z1. As raízes são retornadas em ordem ascendente, ordenadas primeiramente por suas componentes reais e a seguir por suas componentes imaginárias. Se somente uma raíz real for encontrada (i.e. se a=0) então essa raíz é armazenada em z0.

6.4 Equações Cúbicas

Function: int gsl_poly_solve_cubic (double a, double b, double c, double * x0, double * x1, double * x2)

Essa função encontra as raízes reais da equação cúbica,

x³ + a x² + b x + c = 0

com um coeficiente lider sendo a unidade. O número de raízes reais (ou uma ou três) é retornado, e suas localizações são armazenadas em x0, x1 e x2. Se uma raíz real for encontrada então somente x0 é modificada. Quando três raízes reais forem encontradas elas são armazenadas em x0, x1 e x2 em ordem ascendente. O caso de raízes iguais não é considerado especial. Por exemplo, a equação (x-1)^3=0 irá ter três raízes com valores exatamente iguais. Com no caso da quadrática, precisão finita pode causar raízes iguais ou muito próximas para ajustar o eixo real no plano complexo, levando a uma discreta modificação no número de raízes reais.

Function: int gsl_poly_complex_solve_cubic (double a, double b, double c, gsl_complex * z0, gsl_complex * z1, gsl_complex * z2)

Essa função encontra as raízes complexas da equação cúbica,

z³ + a z² + b z + c = 0

O número de raízes complexas é retornado (sempre três) e a localização das raízes está armazenada em z0, z1 e z2. As raízes são retornadas na ordem ascendente, ordenada primeiramente por suas componentes reais e a seguir pelas suas componentes imaginárias.

6.5 Equações Polinomiais Gerais

As raízes de equações polinomiais não podem ser encontradas analiticamente excetuando-se os casos especiais da quadrática, da cúbica e da equação de quarto grau. O algoritmo descrito nessa seção usa um método iterativo para encontrar uma localização aproximada de raízes de polinômios de grau mais elevado.

Function: gsl_poly_complex_workspace * gsl_poly_complex_workspace_alloc (size_t n)

Essa função aloca espaço para uma estrutura gsl_poly_complex_workspace e um espaço de trabalho adequado para resolver um polinômio com n coeficientes usando a rotina gsl_poly_complex_solve.

A função retorna um apontador para o mais recentemente alocado gsl_poly_complex_workspace se nenhum erro for encontrado, e um apontador nulo (null) no caso de erro.

Function: void gsl_poly_complex_workspace_free (gsl_poly_complex_workspace * w): Essa função libera toda memória associada ao espaço de trabalho w.

Function: int gsl_poly_complex_solve (const double * a, size_t n, gsl_poly_complex_workspace * w, gsl_complex_packed_ptr z): Essa função calcula as raízes do polinômio genérico P(x) = a_0 + a_1 x + a_2 x^2 + ... + a_{n-1} x^{n-1} usando redução QR balanceada da matriz adjunta . O parâmetro n especifica o comprimento do vetor estático do coeficiente. O coeficiente do termo de mais alta ordem deve ser diferente de zero. A função requer um espaço de trabalho w de tamanho apropriado. As n-1 raízes são retornadas dentro do vetor estático complexo z de comprimento 2(n-1), alternando partes reais e imaginárias. As funções retornam GSL_SUCCESS se todas as raízes forem encontradas. Se a redução QR não vier a convergir, o controlador de erro é chamado com um código de erro de GSL_EFAILED. Note que devido à precisão finita, raízes de alta mutiplicidade são retornadas como um grupo de raízes simples com exatidão reduzida. A solução de polinômios com raízes de alta ordem requerem algoritmos especializados que tomam a estrutura de multiplicidade dentro de uma conta (veja e.g. Z. Zeng, Algorithm 835, ACM Transactions on Mathematical Software, Volume 30, Impressão 2 (2004), pp 218–236).

6.6 Exemplos

Para demonstrar o uso do resolvedor de polinômios gerais iremos tomar o polinômio P(x) = x^5 - 1 que tem as seguintes raízes,

1, e^{2πi /5}, e^{4πi /5}, e^{6πi /5}, e^{8πi /5}

O seguinte programa irá encontrar essas raízes.

#include <stdio.h>
#include <gsl/gsl_poly.h>

int
main (void)
{
  int i;
  /* coefficients of P(x) =  -1 + x^5  */
  double a[6] = { -1, 0, 0, 0, 0, 1 };  
  double z[10];

  gsl_poly_complex_workspace * w 
      = gsl_poly_complex_workspace_alloc (6);
  
  gsl_poly_complex_solve (a, 6, w, z);

  gsl_poly_complex_workspace_free (w);

  for (i = 0; i < 5; i++)
    {
      printf ("z%d = %+.18f %+.18f\n", 
              i, z[2*i], z[2*i+1]);
    }

  return 0;
}

A saída do programa é,

$ ./a.out 
z0 = -0.809016994374947451 +0.587785252292473137
z1 = -0.809016994374947451 -0.587785252292473137
z2 = +0.309016994374947451 +0.951056516295153642
z3 = +0.309016994374947451 -0.951056516295153642
z4 = +1.000000000000000000 +0.000000000000000000

que concorda com o resultado analítico, z_n = \exp(2 \pi n i/5).

6.7 Referências e Leituras Adicionais

O método QR balanceado e sua análise de erro são ambos descritos no seguintes artigos,

R.S. Martin, G. Peters and J.H. Wilkinson, “The QR Algorithm for Real Hessenberg Matrizes”, Numerische Mathematik, 14 (1970), 219–231.
B.N. Parlett and C. Reinsch, “Balancing a Matrix for Calculation of Eigenvalues and Eigenvectors”, Numerische Mathematik, 13 (1969), 293–304.
A. Edelman and H. Murakami, “Polynomial roots from companion matrix eigenvalues”, Mathematics of Computation, Vol. 64, No. 210 (1995), 763–776.

As fórmulas para diferenças divididas são fornecidas nos seguintes textos,

Abramowitz e Stegun, Handbook of Matematical Functions, Seções 25.1.4 e 25.2.26.
R. L. Burden and J. D. Faires, Numerical Analysis, 9th edition, ISBN 0-538-73351-9, 2011.

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