61.1.1, Iniciando com orthopoly

load ("orthopoly") torna o pacote orthopoly disponível para uso.

Para encontrar o polinómio de Legendre de terceira ordem,

(%i1) legendre_p (3, x);
                      3             2
             5 (1 - x)    15 (1 - x)
(%o1)      - ---------- + ----------- - 6 (1 - x) + 1
                 2             2

Para expressar esse polinómio como uma soma de potências de x, aplique ratsimp ou rat para o resultado anterior.

(%i2) [ratsimp (%), rat (%)];
                        3           3
                     5 x  - 3 x  5 x  - 3 x
(%o2)/R/            [----------, ----------]
                         2           2

Alternativamente, faça o segundo argumento para legendre_p (sua variável “principal”) uma expressão racional canónica (CRE) usando rat(x) em lugar de somente x.

(%i1) legendre_p (3, rat (x));
                              3
                           5 x  - 3 x
(%o1)/R/                   ----------
                               2

Para avaliação em ponto flutuante, orthopoly usa uma análise de erro durante a execução para estimar uma associação superior para o erro. Por exemplo,

(%i1) jacobi_p (150, 2, 3, 0.2);
(%o1) interval(- 0.062017037936715, 1.533267919277521E-11)

intervalos possuem a forma interval (c, r), onde c é o centro e r é o raio do intervalo. Uma vez que Maxima não suporta aritmética sobre intervalos, em algumas situações, tais como em gráficos, vai querer suprimir o erro e sair somente com o centro do intervalo. Para fazer isso, escolha a variável de opção orthopoly_returns_intervals para false.

(%i1) orthopoly_returns_intervals : false;
(%o1)                         false
(%i2) jacobi_p (150, 2, 3, 0.2);
(%o2)                  - 0.062017037936715

Veja a secção veja Avaliação em Ponto Flutuante para maiores informações.

Muitas funções em orthopoly possuem uma propriedade gradef; dessa forma

(%i1) diff (hermite (n, x), x);
(%o1)                     2 n H     (x)
                               n - 1
(%i2) diff (gen_laguerre (n, a, x), x);
              (a)               (a)
           n L   (x) - (n + a) L     (x) unit_step(n)
              n                 n - 1
(%o2)      ------------------------------------------
                               x

A função de um único passo no segundo exemplo previne um erro que poderia de outra forma surgir através da avaliação de n para 0.

(%i3) ev (%, n = 0);
(%o3)                           0

A propriedade gradef somente aplica para a variável “principal”; dderivadas com relação a outros argumentos usualmente resultam em uma mensagem de erro; por exemplo

(%i1) diff (hermite (n, x), x);
(%o1)                     2 n H     (x)
                               n - 1
(%i2) diff (hermite (n, x), n);

Maxima doesn't know the derivative of hermite with respect the first argument
 -- an error.  Quitting.  To debug this try debugmode(true);

Geralmente, funções em orthopoly mapeiam sobre listas e matrizes. Para o mapeamento para avaliação total, as variáveis de opção doallmxops e listarith devem ambas serem true (o valor padrão). Para ilustrar o mapeamento sobre matrizes, considere

(%i1) hermite (2, x);
                                     2
(%o1)                    - 2 (1 - 2 x )
(%i2) m : matrix ([0, x], [y, 0]);
                            [ 0  x ]
(%o2)                       [      ]
                            [ y  0 ]
(%i3) hermite (2, m);
               [                             2  ]
               [      - 2        - 2 (1 - 2 x ) ]
(%o3)          [                                ]
               [             2                  ]
               [ - 2 (1 - 2 y )       - 2       ]

No segundo exemplo, o elemento i, j do valor é hermite (2, m[i,j]); isso não é o mesmo que calcular -2 + 4 m . m, como visto no próximo exemplo.

(%i4) -2 * matrix ([1, 0], [0, 1]) + 4 * m . m;
                    [ 4 x y - 2      0     ]
(%o4)               [                      ]
                    [     0      4 x y - 2 ]

Se avaliar uma função em um ponto fora do seu domínio, geralmente orthopoly retorna uma função não avaliada. Por exemplo,

(%i1) legendre_p (2/3, x);
(%o1)                        P   (x)
                              2/3

orthopoly suporta tradução em TeX; orthopoly também faz saídas bidimensionais em um terminal.

(%i1) spherical_harmonic (l, m, theta, phi);
                          m
(%o1)                    Y (theta, phi)
                          l
(%i2) tex (%);
$$Y_{l}^{m}\left(\vartheta,\varphi\right)$$
(%o2)                         false
(%i3) jacobi_p (n, a, a - b, x/2);
                          (a, a - b) x
(%o3)                    P          (-)
                          n          2
(%i4) tex (%);
$$P_{n}^{\left(a,a-b\right)}\left({{x}\over{2}}\right)$$
(%o4)                         false

61.1.2, Limitations

Quando uma expressão envolve muitos polinómios ortogonais com ordens simbólicas, é possível que a expressão actualmente tenda para zero, e ainda ocorre também que Maxima estar incapacitado de simplificar essa expressão para zero. Se fizer uma divisão por tal quantidade que tende a zero, poderá ficar em apuros. Por exemplo, a seguinte expressão tende para zero para inteiros n maiores que 1, e ainda ocorre também que Maxima está incapacitado de simplificar essa expressão para zero.

(%i1) (2*n - 1) * legendre_p (n - 1, x) * x - n * legendre_p (n, x) + (1 - n) * legendre_p (n - 2, x);
(%o1)  (2 n - 1) P     (x) x - n P (x) + (1 - n) P     (x)
                  n - 1           n               n - 2

Para um n específico, podemos reduzir a expressão a zero.

(%i2) ev (% ,n = 10, ratsimp);
(%o2)                           0

Geralmente, a forma polinomial de um polinómio ortogonal esteja adequada de forma hostil para avaliaçao em ponto flutuante. Aqui está um exemplo.

(%i1) p : jacobi_p (100, 2, 3, x)$

(%i2) subst (0.2, x, p);
(%o2)                3.4442767023833592E+35
(%i3) jacobi_p (100, 2, 3, 0.2);
(%o3)  interval(0.18413609135169, 6.8990300925815987E-12)
(%i4) float(jacobi_p (100, 2, 3, 2/10));
(%o4)                   0.18413609135169

O verdadeiro valor está em torno de 0.184; ess calculo suporta erro de cancelamento por extremo subtrativo.Expandindo o polinómio e então avaliando, fornecendo um melhor resultado.

(%i5) p : expand(p)$
(%i6) subst (0.2, x, p);
(%o6) 0.18413609766122982

Essa não é uma regra geral; expandindo o polinómio não resulta sempre em expressões que são melhores adaptadas a avaliação numérica. Com grande folga, o melhor caminho para fazer avaliação numérica é fazer um ou mais argumentos da função serem números em ponto flutuante. Em função disso, algoritmos especializados em ponto flutuante são usados para avaliação.

A função float do Maxima é até certo ponto indiscriminada; se aplicar float a uma expressão envolvendo um polinómio ortogonal com um grau simbólico ou um parâmetro de ordem, esses parâmetos (inteiros) podem ser convertido em números em ponto flutuante; após o que, a expressão não irá avaliar completamente. Considere

(%i1) assoc_legendre_p (n, 1, x);
                               1
(%o1)                         P (x)
                               n
(%i2) float (%);
                              1.0
(%o2)                        P   (x)
                              n
(%i3) ev (%, n=2, x=0.9);
                             1.0
(%o3)                       P   (0.9)
                             2

A expressão em (%o3) não irá avaliar para um número em ponto flutuante; orthopoly não reconhece valores em ponto flutuante em lugares onde deve haver valores inteiros. Similarmente, avaliação numérica da função pochhammer para ordens que excedam pochhammer_max_index pode ser perturbador; considere

(%i1) x :  pochhammer (1, 10), pochhammer_max_index : 5;
(%o1)                         (1)
                                 10

Aplicando float não avalia x para um número em ponto flutuante

(%i2) float (x);
(%o2)                       (1.0)
                                 10.0

Para avaliar x para um número em ponto flutuante, irá precisar associar pochhammer_max_index a 11 ou mais e aplicar float a x.

(%i3) float (x), pochhammer_max_index : 11;
(%o3)                       3628800.0

O valor padrão de pochhammer_max_index é 100; modifique esse valor após chama orthopoly.

Finalmente, tenha consciência que os livros citados nas referências adotam diferentes definições de polinómios ortogonais; geralmente adotamos as convenções citadas nas convenções de Abramowitz e Stegun.

Antes de suspeitar de um erro no pacote orthopoly, verifique alguns casos especiais para determinar se suas definições coincidem com aquelas usadas por orthopoly. Definitions muitas vezes diferem por uma normalização; ocasionalmente, autores utilizam versões “modificadas” das funções que fazem a família ortogonal sobre um intervalo diferente do intervalo (-1, 1). Para definir, por exemplo, um polinómio de Legendre que é ortogonal a (0, 1), defina

(%i1) shifted_legendre_p (n, x) := legendre_p (n, 2*x - 1)$

(%i2) shifted_legendre_p (2, rat (x));
                            2
(%o2)/R/                 6 x  - 6 x + 1
(%i3) legendre_p (2, rat (x));
                               2
                            3 x  - 1
(%o3)/R/                    --------
                               2

61.1.3, Avaliação em Ponto Flutuante

Muitas funções em orthopoly utilizam análise de erro durante a execução para estimar o erro em avaliações em ponto flutuante; as exceções são funções de Bessel esféricas e os polinómios associados de Legendre do segundo tipo. Para avaliações numéricas, as funções de Bessel esféricas chamam funções da colecção de programas SLATEC. Nenhum método especializado é usado para avaliação numérica dos polinómios associados de Legendre do segundo tipo.

A análise de erro durante a execução ignora erros que são de segunda ordem ou maior na máquina (também conhecida como perda de algarismos). A análise de erro durante a execução também ignora alguns poucos outros tipos de erro. É possível (embora não provável) que o erro actual exceda o estimado.

Intervalos possuem a forma interval (c, r), onde c é o centro do intervalo e r é seu raio. O centro de um intervalo pode sr um número complexo, e o raio é sempre um número real positivo.

Aqui está um exemplo.

(%i1) fpprec : 50$

(%i2) y0 : jacobi_p (100, 2, 3, 0.2);
(%o2) interval(0.1841360913516871, 6.8990300925815987E-12)
(%i3) y1 : bfloat (jacobi_p (100, 2, 3, 1/5));
(%o3) 1.8413609135168563091370224958913493690868904463668b-1

Vamos testar o quanto o erro actual é é menor que o erro estimado

(%i4) is (abs (part (y0, 1) - y1) < part (y0, 2));
(%o4)                         true

Realmente, por esse exemplo o erro estimado é um maior que o erro verdadeiro.

Maxima não suporta aritmética sobre intervalos.

(%i1) legendre_p (7, 0.1) + legendre_p (8, 0.1);
(%o1) interval(0.18032072148437508, 3.1477135311021797E-15)
        + interval(- 0.19949294375000004, 3.3769353084291579E-15)

Um utilizador pode definir operadores aritméticos que fazem matemática de intervalos. Para definir adição de intervalos, podemos definir

(%i1) infix ("@+")$

(%i2) "@+"(x,y) := interval (part (x, 1) + part (y, 1), part (x, 2) + part (y, 2))$

(%i3) legendre_p (7, 0.1) @+ legendre_p (8, 0.1);
(%o3) interval(- 0.019172222265624955, 6.5246488395313372E-15)

As rotinas eseciais em ponto flutuante são chamadas quando os argumentos forem complexos. Por exemplo,

(%i1) legendre_p (10, 2 + 3.0*%i);
(%o1) interval(- 3.876378825E+7 %i - 6.0787748E+7, 
                                           1.2089173052721777E-6)

Let’s compare this to the true value.

(%i1) float (expand (legendre_p (10, 2 + 3*%i)));
(%o1)          - 3.876378825E+7 %i - 6.0787748E+7

Adicionalmente, quando os argumentos forem grandes números em ponto flutuante, as rotinas especiais de ponto flutuante são chamadas; todavia, tos grandes números em ponto flutuante são convertidos para números em ponto flutuante de dupla precisão e o resultado final é número em ponto flutuante de precisão dupla.

(%i1) ultraspherical (150, 0.5b0, 0.9b0);
(%o1) interval(- 0.043009481257265, 3.3750051301228864E-14)

61.1.4, Gráficos e orthopoly

Para desenhar gráficos de expressões que envolvem polinómios ortogonais, deverá fazer duas coisas:

1. Escolher a variável de opção orthopoly_returns_intervals para false,
2. Colocar apóstrofo em qualquer chamada a funções do pacote orthopoly.

Se chamadas a funções não receberem apóstrofo, Maxima irá avaliá-las para polinómios antes de montar o gráfico; consequêntemente, as rotinas especializadas em ponto flutuante não serão chamadas. Aqui está um exemplo de como montar o gráfico de uma expressão que envolve um polinómio de Legendre.

(%i1) plot2d ('(legendre_p (5, x)), [x, 0, 1]), orthopoly_returns_intervals : false;
(%o1)

A expressão completa legendre_p (5, x) recebe apóstrofo; isso é diferente de apenas colocar apóstrofo no nome da função usando 'legendre_p (5, x).

61.1.5, Funções Diversas

O pacote orthopoly define o síbolo de Pochhammer e uma função de passo de unidade. orthopoly utiliza a função delta de Kronecker e a função de passo de unidade em declarações gradef.

Para converter os símbolos Pochhammer em quocientes da funções gama, use makegamma.

(%i1) makegamma (pochhammer (x, n));
                          gamma(x + n)
(%o1)                     ------------
                            gamma(x)
(%i2) makegamma (pochhammer (1/2, 1/2));
                                1
(%o2)                       ---------
                            sqrt(%pi)

Derivadas de símbolos de Pochhammer são fornecidas em termos de psi function.

(%i1) diff (pochhammer (x, n), x);
(%o1)             (x)  (psi (x + n) - psi (x))
                     n     0             0
(%i2) diff (pochhammer (x, n), n);
(%o2)                   (x)  psi (x + n)
                           n    0

É preciso ser cuidadoso com expressões como (%o1); a diferença das funções psi possuem polinómios quando x = -1, -2, .., -n. Esses polinómios cacelam-se com factores em pochhammer (x, n) fazendo da derivada um polinómio de grau n - 1 quando n for um inteiro positivo.

O símbolo de Pochhammer é definido de ordens negativas até sua representação como um quociente de funções gama. Considere

(%i1) q : makegamma (pochhammer (x, n));
                          gamma(x + n)
(%o1)                     ------------
                            gamma(x)
(%i2) sublis ([x=11/3, n= -6], q);
                               729
(%o2)                        - ----
                               2240

Alternativamente, podemos tomar ese resultado directamente.

(%i1) pochhammer (11/3, -6);
                               729
(%o1)                        - ----
                               2240

A função passo de unidade é contínua à esquerda; dessa forma

(%i1) [unit_step (-1/10), unit_step (0), unit_step (1/10)];
(%o1)                       [0, 0, 1]

Se precisar de uma função de degrau unitário que seja ou contínua à esquerda ou contínua à direita do zero, defina a sua própria função usando signum; por exemplo,

(%i1) xunit_step (x) := (1 + signum (x))/2$

(%i2) [xunit_step (-1/10), xunit_step (0), xunit_step (1/10)];
                                1
(%o2)                       [0, -, 1]
                                2

Não redefina a própria unit_step; alguns código em orthopoly requerem que a função de passo de unidade seja contínua à esquerda.

61.1.6, Algorítmos

Geralmente, orthopoly faz avaliações simbólicas pelo uso de uma representação hipergeométrica de polinómios ortogonais. As funções hipegeométricas são avaliadas usando as funções (não documetadas) hypergeo11 e hypergeo21. As excessões são as funções de Bessel metade inteiras e a função de Legendre associada de segundo tipo. As funções de Bessel metade inteiras são avaliadas usando uma representação explícita, e a função de Legendre associada de segundo tipo é avaliada usando recursividade.

Para avaliação em ponto flutuante, nós novamente convertemos muitas fuções em uma forma hipergeométrica; nós avaliamos as funções hipergeométricas usando recursividade para frente. Novamente, as excessões são as funções de Bessel metade inteiras e a função de Legendre associada de segundo tipo. Numericamente, as funções de Bessel meio inteiras são avaliadas usando o código SLATEC.