|
|
|
|
O Cálculo de polinômios ortogonais
Os polinômios de Legendre são definidos por:
p0(x) = 1
p1(x) = x
n*pn(x) = (2*n -1)*x*pn-1(x) - (n - 1)*pn-2(x)
O bojetivo dessa seção é desenvolver programas que calculem o polinômio de Legendre para um dado valor de n.
Para uma primeira tentativa porcuramos seuir fielmente a definição da forma mais fiel possível:
Legendre1(n, x) :=
block ( [],
if n = 0
then 1
else
if n = 1
then x
else ((2*n - 1)*x*Legendre1 (n - 1, x)
- (n - 1) *Legendre1 (n - 2, x)) / n
)
Um primeiro teste fornece:
Legendre1(3, z);
2 5 z (3 z - 1) -------------- - 2 z 2 (%o2) -------------------- 3 Legendre1 (5, t);
2 5 t (3 t - 1) 7 t (-------------- - 2 t) 2 2 2 3 (3 t - 1) 5 t (3 t - 1) 9 t (-------------------------- - ------------) 4 (-------------- - 2 t) 3 2 2 ----------------------------------------------- - ------------------------ 4 3 (%o3) -------------------------------------------------------------------------- 5
Perfeito para os amantes de expressões imensas, mas longe da perfeição para aqueles que preferem resultados simplificados. O que fizemos errado?
Nada. Simplificação é uma coisa difícil de fazer e simplesmente ignoramos o problema. Multiplicamos e adicionamos polinômios. Para obter uma representação canônica de polinômios temos que expandir todos os produtos:
Legendre2(n, x) :=
block ( [],
if n = 0
then 1
else
if n = 1
then x
else expand(((2*n - 1)*x*Legendre2 (n - 1, x)
- (n - 1) *Legendre2 (n - 2, x)) / n) )
Legendre2(3, z);
3
5 z 3 z
(%o5) ---- - ---
2 2
Legendre2(5, z);
5 3
63 z 35 z 15 z
(%o6) ----- - ----- + ----
8 4 8
Expandimos a expressão inteira - Essa é a aproximação segura, mas pode ter um custo muito alto. Podemos querer perguntar sobre simplificações com o menor esforço. Isso é suficiente expandir os dois produtos separadamente:
Aqui expandimos somente a expressão que contém uma multiplicação com a variável:
Legendre2(n, x) :=
block ( [],
if n = 0
then 1
else
if n = 1
then x
else (expand((2*n - 1)*x*Legendre2 (n - 1, x)) - (n - 1) *Legendre2 (n - 2, x)) / n )
Isso é melhor:
Legendre2(n, x) :=
block ( [],
if n = 0
then 1
else
if n = 1
then x
else expand((2*n - 1)/n*x*Legendre2 (n - 1, x))
- expand((n - 1)/n *Legendre2 (n - 2, x)) )
Uma olhada restrita no código reela que ele é ineficiente e nem mesmo é elegante como se pode pensar inicialmente. A falha principal é que um monte de cálculos são realizados repetidamente . Para ver isso, você pode usar a facilidade trace do Maxima:
trace(Legendre2);
[Legendre2] Legendre2(7, z)
Isso escreverá uma mensagem todas as vezes que a função Legendre2 for iniciada e concluída.
***
O caminho fácil para evitar isso é evitar ciclos aninhados. Isso é o que iremos fazer adiante.
Vamos agora escrever uma função não cíclica que executa o mesmo cálculo:
LegendreN(n, x) :=
block ( [p0, p1, pn, cnt],
if n = 0
then return (1)
else if n = 1
then return (x),
p0 : 1,
p1 : x,
cnt: 2,
while cnt <= n do
( pn: expand(((2*cnt - 1)*x*p1
- (cnt - 1) *p0) / cnt),
p0: p1,
p1: pn,
cnt: cnt + 1
),
pn
)
Aqui está uma variante que usa uma declaração composta - que é, um bloco - dentro de uma cláusula while.
LegendreN(n, x) :=
block ( [p0, p1, pn, cnt],
if n = 0
then return (1)
else if n = 1
then return (x),
p0 : 1,
p1 : x,
cnt: 2,
while block(pn: expand(((2*cnt - 1)*x*p1
- (cnt - 1) *p0) / cnt),
cnt: cnt + 1,
cnt <= n)
do
(
p0: p1,
p1: pn
),
pn
)
Isso é possível também programar uma repetição com uma instrução de salto. existe uma boa razão para preferir uma variante adequada para a declaração do, mas um programador deve estar apto a entender programas Maxima com instruções de salto.
LegendreN(n, x) :=
block ( [p0, p1, pn, cnt],
if n = 0
then return (1)
else if n = 1
then return (x),
p0 : 1,
p1 : x,
cnt: 2,
beginOfLoop,
pn: expand(((2*cnt - 1)*x*p1
- (cnt - 1) *p0) / cnt),
p0: p1,
p1: pn,
cnt: cnt + 1,
if cnt <= n then go(beginOfLoop),
pn
);
Saltando agora, escrevemos definições que constóem seqüências de
polinômios de grau crescente. S ào todos esses polinômios necessários?
Não, eles não são.
O que é preciso para os resultados intermediários são os
coeficientes de todos esses polinômios - nada mais.
O seguinte exemplo usa três listas para armazenar os coeficientes dos
três polinômios pn, pn-1, pn-2:
LegendreNN(n, x) :=
block ( [cnt, pn, coeffsP0, coeffsP1, coeffsPN, oldList ],
if n = 0
then return (1)
else if n = 1
then return (x),
coeffsP0: makelist (0, x, 0, n),
coeffsP1: makelist (0, x, 0, n),
coeffsPN: makelist (0, x, 0, n),
coeffsP0[1]: 1,
coeffsP1[2]: 1,
cnt: 2,
while cnt <= n do
(coeffsPN[1]: -coeffsP0[1]*(cnt - 1)/cnt,
for idx : 2 thru cnt + 1 do
coeffsPN[idx] : ((2*cnt - 1)*coeffsP1[idx - 1]
-(cnt - 1)*coeffsP0[idx])/ cnt,
oldList: coeffsP0,
coeffsP0: coeffsP1,
coeffsP1: coeffsPN,
coeffsPN: oldList,
cnt: cnt + 1
),
pn: 0,
for idx:1 thru n + 1 do
pn: pn + coeffsP1[idx]*x^(idx - 1),
pn
);
Para comparar os tempos de computação solicitemos ao Maxima que mostre todos os tempos:
showtime:all;
Agora podemos executar as definições que queremos comparar. Para suprimir a saída desse grande polinômio, encerramos nossa entrada com um sinal de dólar:
(%i4) LegendreN(200, x)$
Evaluation took 36.58 seconds (36.58 elapsed)
(%i5) LegendreNN(200, x)$
Evaluation took 13.62 seconds (13.62 elapsed)
Concluímos que aritmética de polinômios e expansões freqüêntes de polinômios consome muito tempo em computações.
|
|
|
|