風の迷路

移転先:http://fermiumbay13.hatenablog.com/ ここはヤフーブログ終了に伴い12/15閉鎖します

全体表示

[ リスト ]

ルジャンドル多項式

物理でよくでてくるというルジャンドルの微分方程式というものがあるようですが、
その特別な場合の解である多項式をルジャンドル多項式と呼びます。
このルジャンドル多項式を応用すれば、定積分値がいとも簡単に求まってしまうわけです。
まあ……人間的な方法ぢゃないですけどねぇ。

ルジャンドル多項式とは、以下のような漸化式で求まる多項式です。
P0(x)=1
P1(x)=x
P_(n+1)(x)=(2n+1)/(n+1) xPn(x)−n/(n+1) P_(n-1)(x)
例えばP2を求めたければ、↑の式のn=1を代入して、
P2(x)=(2×1+1)/(1+1) xP1(x)−1/2 P0(x)
=3/2 x^2−1/2
1/2 (3x^2−1)
となりました。
このようにしていけば、すべてのPnを求めていくことが可能です。
しかしこれでP50(x)とか求めようとするとえらい大変なので困っていたのですが、
どうやら一般式があったそうで……ロドリグの公式というみたいです。

Pn(x)=1/(2^n n!)(d/dx)^n {(x^2−1)^n}

なんて便利な……良くできてるなあ。
例えばP2(x)を求めたかったらn=2を代入すれば良いだけです。
P2(x)=1/(2^2×2!)(d/dx)^2 {(x^2−1)^2}
=1/8(d/dx)^2(x^4−2x^2+1)
=1/8(d/dx)(4x^3−4x)
=1/8 (12x^2−4)
1/2 (3x^2−1)
出ました。


このようなルジャンドル多項式を、定積分値を求めることにも応用できます。
ガウス求積と呼ばれる手法です。

∫[-1→1]f(x)dx≒Σ[i=1→n]Aif(xi)

関数f(x)を範囲−1〜1で定積分したときの値が、単なる足し算で求まるわけです。
ここで、Aiとは重みと呼ばれる係数で、次の式で求められます。
Ai=2/{(1−xi^2)(P'n(xi))^2}
P'n(x)とは、ルジャンドル多項式のPn(x)をxで1階微分したものです。
xiとは、ルジャンドル多項式の零点になります。
零点とはPn(xi)=0としたときの解xiのことです。
ルジャンドル多項式Pn(x)はn次多項式になりますので、解xiはn個あります。
それぞれに対する重みAiが求まりますので、全部でn個の重みになります。
「積分値を求めたい関数に零点を代入したもの」と「その零点に対する重み」との積の総和が、
求めたい積分値になるわけです。

実際に何かやってみましょうね。

∫[0→3]x^2 dx

これは手計算でも簡単に求めることが出来ます。
∫[0→3]x^2 dx=[x^3/3](範囲:0→3)=9
になります。

これをガウス求積で求めることにします。
ガウス求積で求めるには積分範囲を−1〜1に改めなければなりませんので、
変数変換をして積分範囲を変えることにします。
一般に、∫[a→b]f(x)dxを求める際には、
t=2x/(b−a)+(a+b)/(a−b)と置換してやれば、範囲が−1〜1になります。
なんか激しいですが、連立方程式を解けば求まります。
今はa=0、b=3ですから、
t=2/3 x−1 と置換すればよいことになりますね。すると
x=3/2 (t+1)
dx=3/2 dt
∫[0→3]x^2 dx=3/2 ∫[-1→1]{3/2 (t+1)}^2 dt
と変形されます。範囲が−1〜1になったので、ガウス求積が可能です。
ややこしいので、積分の中身の{3/2 (t+1)}^2=f(t)とします。

精度を初めに決めます。ルジャンドル多項式Pn(x)の、
nを大きくすればするほど高精度になり、その代わり計算量が増えます。
n=3を選択してみましょうね。

まずは零点を求めます。
P3(x)=1/2 (5x^3−3x) ですから、零点は
1/2 (5x^3−3x)=0 の解xとなります。3つあるはずです。
5x^3−3x=0
x(x^2−3/5)=0
x{x−√(3/5)}{x+√(3/5)}=0
となりますので、x=0,±√(3/5)となります。x1=0,x2=√(3/5),x3=−√(3/5)とします。
重みを求めるため、ルジャンドル多項式の微分をします。
P'3(x)=1/2 (15x^2−3)ですね。
これをもとに、上の方の式を使って、それぞれに対する重みを求めます。
A1=2/〔(1−0^2){1/2 (15×0−3)}^2〕=8/9
A2=2/〔(1−3/5){1/2 (15×3/5−3)}^2〕=5/9
A3=2/〔(1−3/5){1/2 (15×3/5−3)}^2〕=5/9
となりました。
あとは、もとの関数に零点を代入したものと、それに対する重みを掛け算したものをすべて足すだけです。

3/2 ∫[-1→1]f(t) dt
≒3/2 {A1f(x1)+A2f(x2)+A3f(x3)}
=3/2 {8/9 f(0)+5/9 f(√(3/5))+5/9 f(−√(3/5))}
=3/2 {8/9 f(0)+5/9 f(√(3/5))+5/9 f(−√(3/5))}
=3/2 {8/9 {3/2 (0+1)}^2+5/9 {3/2 (√(3/5)+1)}^2+5/9 {3/2 (−√(3/5)+1)}^2}
あとは掛け算とか足し算するだけです……
≒3/2 {2+3.9365+0.0635}=9
確かに9って求まりましたねぇ。n=3ぐらいになると結構精度が良くなります。


こんなこと実に長ったらしくて全然役に立たないだろうと思われますが、
この手法の良いところは、どんな関数でも零点や重みは不変で、
実際それを使い回しすれば、後は足し算掛け算とかで積分が求まっちゃう所にあります。

例えば、∫[-1→1]e^x dx を求めるとします。
範囲が−1〜1になっている場合は変数変換の必要はありません。
ここでn=3におけるガウス求積を求めようと思えば、零点と重みを使い回しして
∫[-1→1]e^x dx≒A1e^(x1)+A2e^(x2)+A3e^(x3)
=8/9 e^(0)+5/9 e^(√(3/5))+5/9 e^(−√(3/5))
を計算すればよいことになります。関数電卓があればすぐでしょう。
これは、2.350337ぐらいになります。よって、∫[-1→1]e^x dx≒2.350337
真値は2.350402なので、かなり真値に近いことが分かります。

同じく定積分値を求める手法として台形和やシンプソンの公式というものもありますが、
足し算をしなければならない項がとても多くなってしまうので、電卓で打とうにも大変になります。
一方、ガウス求積は必要な精度に相当する重みや係数が分かっていれば、
それらの値を使ってすぐに求められるのが特徴です。

高精度にするなら、高次のルジャンドル多項式を求めて、零点と重みを求めて
表にでもしておけば良いですね^^;

閉じる コメント(4)

顔アイコン

なんと物理でよく出てくるのですか@@ 自分の知っている物理には登場しませんでした。
数値が大きくなればなるほど求めづらくなるのは仕方がないですが、これは大変です(*_*;

一度求めてしまえば使い回しができるとはかなり便利そうです!
関数電卓を使用して求めるにはうってつけなんですね。

何やら見たことがない式はもちろん、記号もあるのですが全て暗記なさっているのでしょうか@@
これらは一体どんな物理に使われているのでしょう(^_^;)

2013/1/4(金) 午前 2:21 [ Dik ] 返信する

顔アイコン

よく使われているは、多分結構嘘です。何かで出てくるんじゃないですか;
2階の漸化式なので、やたら面倒です。一般式があって助かりました。

面倒ですが、初めの計算さえしておけば後は掛け算と足し算でいいので結構便利です。
関数の積分を求めるには関数の計算が必要なので……手計算はめんどくさそうです。

そんなばかな、確認のために調べてます; ウィキペディアは案外便利ですよ^^;
量子力学でルジャンドル関係の何かが出てきましたが、
はて、この微分方程式は使うのかなあ。よく知らないのです……

2013/1/5(土) 午前 0:01 フェルミウム湾 返信する

顔アイコン

課題でルジャンドルの式たちが出てきてかなりナウな記事ですよこれ…

2013/1/6(日) 午前 2:03 [ _ ] 返信する

顔アイコン

貴殿もルジャンドルの使い手であったか……@@
そうすると既にこの便利さに気付いていらっしゃることでしょう。

2013/1/6(日) 午後 11:46 フェルミウム湾 返信する

コメント投稿

顔アイコン

顔アイコン・表示画像の選択

名前パスワードブログ
絵文字
×
  • オリジナル
  • SoftBank1
  • SoftBank2
  • SoftBank3
  • SoftBank4
  • docomo1
  • docomo2
  • au1
  • au2
  • au3
  • au4
投稿

.


プライバシー -  利用規約 -  メディアステートメント -  ガイドライン -  順守事項 -  ご意見・ご要望 -  ヘルプ・お問い合わせ

Copyright (C) 2019 Yahoo Japan Corporation. All Rights Reserved.

みんなの更新記事