Legendre関数

 

P(x,n):=という書き方だと、subsを使うことになるのでなにかと不便です。

P(x,n):=1/2^n/n!*diff((x^2-1)^n,x$n);

print(P(x,0));
for i from 0 to 3 do
  print(P(x,i)=expand(subs(P(x,n),n=i)));
end_for

diff((x^2 - 1)^n, x $ n)/(2^n*n!)
1
1 = diff(1)
x = diff(x^2 - 1, x)/2
(3*x^2)/2 - 1/2 = (3*x^2)/2 - 1/2
(3*x*(x^2 - 1))/2 + x^3 = (5*x^3)/2 - (3*x)/2

0と1の時は何故か途中で計算が終わってしまいました。

 

P:=(x,n)->という書き方が楽ちんですが、P(x,0)といった表示ができなくなる欠点があります。

そこでhold関数を使います。

P:=(x,n)->1/2^n/n!*diff((x^2-1)^n,x$n);
for i from 0 to 3 do
  print(hold(P)(x,i)=expand(P(x,i)));
end_for

(x, n) -> ((1/2^n)/n!)*diff((x^2 - 1)^n, x $ n)
P(x, 0) = 1
P(x, 1) = x
P(x, 2) = (3*x^2)/2 - 1/2
P(x, 3) = (5*x^3)/2 - (3*x)/2

0から9まで増やしてみます。

for i from 0 to 9 do
  print(hold(P)(x,i)=expand(P(x,i)));
end_for

P(x, 0) = 1
P(x, 1) = x
P(x, 2) = (3*x^2)/2 - 1/2
P(x, 3) = (5*x^3)/2 - (3*x)/2
P(x, 4) = (35*x^4)/8 - (15*x^2)/4 + 3/8
P(x, 5) = (63*x^5)/8 - (35*x^3)/4 + (15*x)/8
P(x, 6) = (231*x^6)/16 - (315*x^4)/16 + (105*x^2)/16 - 5/16
P(x, 7) = (429*x^7)/16 - (693*x^5)/16 + (315*x^3)/16 - (35*x)/16
P(x, 8) = (6435*x^8)/128 - (3003*x^6)/32 + (3465*x^4)/64 - (315*x^2)/32 + 35/128
P(x, 9) = (12155*x^9)/128 - (6435*x^7)/32 + (9009*x^5)/64 - (1155*x^3)/32 + (315*x)/128

 

グラフにしてみました。

for i from 0 to 7 do
  f:=subs(P(x,n),n=i);
  print(hold(P)(x,i)=f);
  plot(plot::Function2d(f,x=-1..1, Color=RGB::Blue));
end_for

P(x, 0) = 1
MuPAD graphics
P(x, 1) = x
MuPAD graphics
P(x, 2) = (3*x^2)/2 - 1/2
MuPAD graphics
P(x, 3) = (3*x*(x^2 - 1))/2 + x^3
MuPAD graphics
P(x, 4) = 3*x^2*(x^2 - 1) + (3*(x^2 - 1)^2)/8 + x^4
MuPAD graphics
P(x, 5) = (15*x*(x^2 - 1)^2)/8 + 5*x^3*(x^2 - 1) + x^5
MuPAD graphics
P(x, 6) = (15*x^4*(x^2 - 1))/2 + (5*(x^2 - 1)^3)/16 + (45*x^2*(x^2 - 1)^2)/8 + x^6
MuPAD graphics
P(x, 7) = (35*x*(x^2 - 1)^3)/16 + (21*x^5*(x^2 - 1))/2 + (105*x^3*(x^2 - 1)^2)/8 + x^7
MuPAD graphics

Legendre関数の直交性の検討をします。P(x,○)とP(x,△)を掛け算したものを-1～1で積分しました。

hold(int(P(x,1)*P(x,1),x=-1..1))=int(P(x,1)*P(x,1),x=-1..1);
hold(int(P(x,0)*P(x,4),x=-1..1))=int(P(x,0)*P(x,4),x=-1..1);
hold(int(P(x,7)*P(x,9),x=-1..1))=int(P(x,7)*P(x,9),x=-1..1);
hold(int(P(x,13)*P(x,17),x=-1..1))=int(P(x,13)*P(x,17),x=-1..1);
hold(int(P(x,15)*P(x,15),x=-1..1))=int(P(x,15)*P(x,15),x=-1..1);

int(P(x, 1)*P(x, 1), x = -1..1) = 2/3
int(P(x, 0)*P(x, 4), x = -1..1) = 0
int(P(x, 7)*P(x, 9), x = -1..1) = 0
int(P(x, 13)*P(x, 17), x = -1..1) = 0
int(P(x, 15)*P(x, 15), x = -1..1) = 2/31

xはcosθとすることが多いです。

gth:=Symbol::theta;
P:=(gth,n)->expand(Simplify(1/2^n/n!*diff((cos(gth)^2-1)^n,gth$n)));
for n from 0 to 5 do
  print(hold(P)(gth,n)=P(gth,n));
end_for

`θ`
(gth, n) -> expand(Simplify(((1/2^n)/n!)*diff((cos(gth)^2 - 1)^n, gth $ n)))
P(`θ`, 0) = 1
P(`θ`, 1) = -cos(`θ`)*sin(`θ`)
P(`θ`, 2) = (3*sin(`θ`)^2)/2 - 2*sin(`θ`)^4
P(`θ`, 3) = 2*cos(`θ`)*sin(`θ`)^3 - (9*cos(`θ`)^3*sin(`θ`)^3)/2
P(`θ`, 4) = (32*sin(`θ`)^8)/3 - (175*sin(`θ`)^6)/12 + (35*sin(`θ`)^4)/8
P(`θ`, 5) = (625*cos(`θ`)^3*sin(`θ`)^7)/24 + (113*cos(`θ`)*sin(`θ`)^7)/24 - (63*cos(`θ`)*sin(`θ`)^5)/8

グラフにしました。

for n from 0 to 5 do
  f:=P(gth,n):
  print(hold(P)(gth,n)=f);
  plot(plot::Function2d(f,gth=-PI/2..PI/2,Color=RGB::Blue));
end_for

P(`θ`, 0) = 1
MuPAD graphics
P(`θ`, 1) = -cos(`θ`)*sin(`θ`)
MuPAD graphics
P(`θ`, 2) = (3*sin(`θ`)^2)/2 - 2*sin(`θ`)^4
MuPAD graphics
P(`θ`, 3) = 2*cos(`θ`)*sin(`θ`)^3 - (9*cos(`θ`)^3*sin(`θ`)^3)/2
MuPAD graphics
P(`θ`, 4) = (32*sin(`θ`)^8)/3 - (175*sin(`θ`)^6)/12 + (35*sin(`θ`)^4)/8
MuPAD graphics
P(`θ`, 5) = (625*cos(`θ`)^3*sin(`θ`)^7)/24 + (113*cos(`θ`)*sin(`θ`)^7)/24 - (63*cos(`θ`)*sin(`θ`)^5)/8
MuPAD graphics