連続関数に、数値積分 Newton-Cotes 法、Legendre 法、Gauss-Chebyshev 法適用し、解析値との比較で、手法の優劣比較を行った。これを、具体的な数値積分計算 1 に示す。計算精度の良い順は、Legendre 法、Gauss-Chebyshev 法、Newton-Cotes法の順である。
数値積分には、積分域 [a.b] を等間隔で行う方法と不等間隔で積分を行う方法がある。
f(x) の積分域 [a.b] を等間隔で積分する方法は、Newton-Cotes 法と呼ばれ、台形則、Simpson1/3 則、
Simpson3/8 則、Boole 則がある。これらは、区分点(分点)の数が異なる。
不等間隔には、Legendre 法と Gauss-Ckebyshev 法がある。Legendre 法は、その次数毎に非等間隔の区分点
(分点) xi の f( xj) に異なる重み係数を設けて行う。 Gauss-Chebyshev 法はその次数毎に分点は非等間隔である
が、 f(xj ) の重み係数は同じである。
なお、Legendre 法は、前章の数値計算 5 直交関数に記載した Legendre 多項式から導きだされるが、Gauss-
Chebyshev 法は前章の Chebyshev 多項式とは関係がない。
その外、無限大領域までを扱う計算方法がある。ここでは触れない。
これらの数値計算法により、Chebyshev多項式の3次、12次、25次について、数値積分を行い優劣を比較した。
これを、具体的な数値積分計算 1 に示す。
Newton-Cotes 法は、積分域の区間を等間隔に分けて、積分する。これらの有名な方法を近似の次数毎に分け、
被積分関数を f(x) とすると、以下となる。
積分域は [x0,x1] で、その幅は h=(x1-x0) である。
1)1次式 :台形則
h/2・(f(x0)+f(x1)) 8-N-1 式
2)2次式 :Simpson 1/3 則
h/2・(f(x0)+4f(x0+h/2)+f(x1) )・1/3 8-N-2 式
3)3次式 :Simpson 3/8 則
h/3・(f(x0)+3f(x0+h/3)+3f(x0+2h/3)+f(x1) )・3/8 8-N-3 式
4)4次式 :Boole 則
h/4・(7f(x0)+32f(x0+h/4)+12f(x0+h/2)+32f(x0+3h/4)+7f(x1))・2/45
8-N-4 式
このように、区分点での f(x) に定数を掛けて、合算すると積分値が算出できる手法である。
また、不等間隔の離散データに対する積分方法は、次章の数値計算 7 に示す。
参考図書;スミルノフ高等数学教程 T巻 第一分冊 P.252
(110.節にシンプソンの公式の説明があり、等間隔の場合について記述されている。)
1、計算の概要
この方法は、積分域 [a,b] を [-1,1] に変換し、Legendre 多項式の解を分点
xk として、分点
xk に対応し
た重み係数を用いる。
8-L-1 式
上記の f(
xj ) の算出は、[-1,1] での Legendre 多項式 P
n の解
xk に領域 [a,b]から変換した
xj を 関数 f(
xj )
に代入し、行う。
8-L-2 式
重み係数 w
j は Legendre 多項式の解である分点 x
k に対応した値で、次の式から算出される。この値は、
積分域 [a,b] には関係がなく、x
k のみによる。従い、Legendre 多項式 P
n の次数毎に、同じ解と同じ重み係
数を持つ事になる。一旦、それぞれの値を知れば、被積分対象の関数に関係なく使用できる。
8-L-3 式
あるいは、
8-L-4 式
上記の計算が可能となる理論的な考えは、多くのweb上でも公開されている。その概要は、f(x) を 2n-1 次
の多項式とし、P
n と Q
n-1 の積と剰余式 R(x) で展開できるとする。P
n(x) の全ての解で成り立つラグランジ
ェ多項式で、f(x) を近似できるので、直交性をもつ Legendre 多項式の解から分点と重み係数を得る事が出
来るという考えである。Legendre 多項式の直交性は前章の
数値計算 5 直交関数に示す。
これらの考えでラグランジェ多項式から Legendre 多項式への転換とその逆が跳躍点である。
2、分点 xk と重み係数 wj 算出
分点 x
k はLegendre多項式の解である。Legendre 多項式 P
n の具体的な式は、
数値計算 5 直交関数に
34 次まで示している。ご参照していただきたい。
6次までのP
nは、以下である。
2次
8-L-5 式
3次
8-L-6 式
4次
8-L-7 式
5次
8-L-8 式
6次
8-L-9 式
上記の式の解、すなわち分点は下表の通りである。
なお、上記の6次 P
6 は解に以下を持つ。
8-L-10 式
しかし、実用的でないので、0.932469514203152とした。
上述の解に対応した重み係数を求めるには、導関数が必要となる。具体的には、4次の場合、P4 の導関数は次
式となる。
8-L-11 式
従い、8-L-3 式は、以下の式で示される。
8-L-12 式
4次の解 x
k=sqrt((15.+2*sqrt(30.))/35) の時は、以下で算出できる。
8-L-13 式
上記のようにして得た P
n の6次までの重み係数を示す。
また、[-1,1] の積分域の中央部の分点の重み係数は他の分点の場合より大きい。Newton-Cotes法の場合も
同じである。
なお、x
k の領域 [a,b] への変換値 x
j の具体的な算出については、
具体的な数値積分計算1に示す。
1、計算の概要
このChebyshev法は、非等間隔の分点 x
j を持つが、重み係数 w
j がその次数で同一である特長を持つ。
f(x
j) と w
j の積の回数が少ない利点を持つ。また、Chebyshev多項式(前章に記載)とは関係がない。
この計算の考えは、
・f(x) を [-1,1] で積分した値が、f(x
1),f(x
2),... f(x
n) の和に定数 k を掛けた値にほぼ等しいと仮定する。
8-C-1式
ただし、上式で R(x) は残余とする。
・この f(x) を (n+1) 階まで導関数を持つとし、マクローリン展開を行う。
8-C-2式
・8-C-1式の右辺のΣの部分を求める。上式の左辺 f(x) の x に x
1,x
2,.. x
n を代入すると以下となる。
8-C-3式
8-C-4式
・・・
8-C-5式
・上式の和 f(x1)+f(x2)+・・・+f(xn) に k を掛けると、8-C-1 式の右辺が求まり、 以下の式となる。
8-C-6式
・8-C-1 式の左辺の [-1,1] での f(x) の積分を求める。 8-C-2 式のマクローリン展開部を積分する。
8-C-7式
次に定数である f(0),f'(0) 等を積分記号の前に出すと
8-C-8式
x の奇数次の積分は、x の偶数次になり、積分域が [-1,1] 故に (-1)
2p - (1)
2p = 0 となり、以下となる。
8-C-9a式
8-C-9b式
8-C-9式
また、xの偶数次の積分は、以下となる。
8-C-10式
結果、8-C-1 式の f(x) の [-1,1] 域の積分は以下となる。
8-C-11式
2、Chebyshev法の数値積分の条件算出
・8-C-1 式の両辺、8-C-6 式と8-C-11式の比較により、次の式(条件)を得る事が出来る。
1)
重み係数となる
(ただし、nは次数あるいは分点数) 8-C-12式
2) 奇数次の xj について
8-C-13式
3)偶数次の xj について
8-C-14式
3、Chebyshev法の重み係数と分点の算出方法
上記 3 条件を満足する x
j を求める事が出来れば、f(x) の積分を、下式のように、近似できる事になる。
8-C-15式
すなわち、x
j を含んだ係数を持つ式として
2 次式 ( 2分点 )は以下の式となる。
8-C-16式
8-C-17式
3次式(3分点)は以下となる。
8-C-18式
8-C-19式
n 次式(n分点)は以下で、それぞれの係数に x
j を含んでいる式となる。
8-C-20式
8-C-21式
上式のそれぞれの係数の値を8-C-12 式、8-C-13 式、8-C-14 式から求める方法を次々節に示す。
4、Chebyshev法の重み係数と分点
・重み係数:
ただし、nは次数 8-C-12式
・2次(
2分点): 2分点を得る式
8-C-22式
・3次
(3分点): 3分点を得る式
8-C-23式
・4次(
4分点): 4分点を得る式
8-C-24式
・5次(
5分点) : 5分点を得る式
8-C-25式
・6次(
6分点):6分点を得る式
8-C-26式
上記のまとめとして、
Chebyshev積分の重み係数と分点を下表に示す。
なお、G
6(x)の解に、以下を得たが、実用的でないので、0.866246818107821とした。
8-C-27式
5、分点を求める式の算出( 8-C-22 式から 8-C-26 式 の根拠)
・
2次式( 8-C-22 式 の算出根拠)
x の 2 次での分点を求める式は、8-C-20 式の 2 次形式で、以下である。
8-C-28式
上式を展開すると、
8-C-29式
2次式では、n=2であるので、重み係数は
k=2/2=1 8-C-30式
第 1 項の係数は、以下となる。
(x
1+x
2)=0 8-C-31式
∵8-C-13 式より (x
1+x
2)=0 8-C-32式
第 2 項の係数は、以下となる。
2・x1・x2 = (x
1+x
2)
2- (x
12+x
22) = 0 - 2/3 = - 2/3 8-C-33式
∵8-C-14式より x
12+x
22 = 2/(2*1+1) 8-C-34式
従い、
8-C-22式
・
3次式( 8-C-23 式の算出根拠)
x の 3 次での分点を求める式は、8-C-20 式の 3 次形式で、以下である。
8-C-35式
上式を展開すると、
8-C-36式
3次式では、n=3 であるので、重み係数は
k=2/3 8-C-37式
第 1 項の係数は、
(x
1+x
2+x
3)=0 8-C-38式
∵8-C-13式より (x
1+x
2+x
3)=0
第 3 項の係数は、
2・(x
1・x
2+x
2・x
3+x
3・x
1)= (x
1+x
2+x
3)
2- (x
12+x
22+x
32) = 0 -1 = -1
8-C-39式
第 4 項の係数は、
-3・x
1・x
2・x
3=
(x
1+x
2+x
3)
3- (x
13+x
23+x
33) - 3 (x
1+x
2+x
3)・(x
1・x
2+x
2・x
3+x
3・x
1)
=0 - 0- 3・0・(-1/2)= 0 8-C-40式
∵8-C-13式より x
13+x
23+x
33=0 8-C-41式
従い、
8-C-23式
・
4次式( 8-C-24 式の算出根拠)
x の 4 次での分点を求める式は、8-C-20 式の 4 次形式で、以下である。
8-C-42式
上式を展開すると
8-C-43式
4次式では、n=4 であるので、重み係数は
k=1/2
第 2 項の係数は、
(x
1+x
2+x
3+x
4)=0 (∵8-C-13式より 0 ) 8-C-44式
第 3 項の係数は、
2・'(x
1・x
2+x
2・x
3+x
3・x
4+x
4・x
1)=(x
1+x
2+x
3+x
4)
2- (x
12+x
22+x
32+x
42)
= 0 - 4/3 = - 4/3 8-C-45式
∵8-C-14式より x
12+x
22+x
32+x
42 = 2・2/1/(2*1+1) = 4/3
8-C-46式
第 4 項の係数は、
-3(x
1・x
2・x
3+x
2・x
3・x
4+x
3・x
4・x
1+x
4・x
1・x
2)
=(x
1+x
2+x
3+x
4)
3
- (x
13+x
23+x
33+x
43)
-3 (x
1+x
2+x
3+x
4)・(x
1・x
2+x
2・x
3+x
3・x
4+x
4・x
1+x
4・x
2+x
1・x
3)
= 0
- 0
- 0
=0 8-C-47 式
∵8-C-13 式より x
1+x
2+x
3+x
4=0
x
13+x
23+x
33+x
43=0
第 5 項の係数は、
-12・x
1・x
2・x
3・x
4
=(x
1+x
2+x
3+x
4)
4
+ 3(x
14+x
24+x
34+x
44)
- 6(x
1・x
2+x
2・x
3+x
3・x
4+x
4・x
1+x
4・x
2+x
1・x
3)
2
- 4(x
13+x
23+x
33+x
43)(x
1+x
2+x
3+x
4)
= 0
+3・(4/5)
-6・(-2/3)・(-2/3)
-4・0
=-4/15
∵8-C-14式より x
14+x
24+x
34+x
44=2・2/1/(2*2+1)=4/5
また、第3項の計算結果を上式の三番目の項に代入する。
∴ x
1・x
2・x
3・x
4=1/45
従い、
8-C-48 式
・
5次式( 8-C-25 式の算出根拠)
この式の各項の係数の算出過程と値は、改訂前は示したが、今は示さない。後日に示す。
現在、数式の貼り付けがいつのまにか出来なくなっているので、調査中 2014/2/7
・
6次式( 8-C-26 式の算出根拠)
計算続行中
6、その他
積分域 [a,b] への [-1,1] 域の分点の変換は具体的な数値積分計算の個所で示す。
参考図書:数値計算ハンドブック
新数値計算法 小国力著、サイエンス社出版
1、数値積分計算 (コード1)
数値積分の手法には、前述のように、各種ある。どの方法が一般的な関数に対して効果的であるかを知るために、具体的に計算を行った。
連続関数である Chebyshev多項式に、数値積分
Newton-Cotes 法、Legendre 法、Gauss-Chebyshev 法適用し、解析値との比較で、手法の優劣比較を行う。
結論を言うと、計算精度の良い順は、Legendre 法、Gauss-Chebyshev 法、Newton-Cotes法の順である。
2、数値積分計算の説明
1) [-1,1] 内の点から積分域[a,b] への点への変換
この積分域の変換は、8-L-2 式の x
j = (b+a)/2 + (b-a)x
k/2 を用いて行う事になる。
b-a = h の時、b = a+h で、(b+a)/2 = (h+2a)2 = a+h/2 となる。積分の下限点 a に h/2 をプラスす
ればよい事になる。従い、以下となる。
2)数値積分の対象
・数値積分の手法
Gauss-Legendre法 :2、3、4、5、6 分点(次数)
Chebyshev法 :2、3、4、5、6 分点(次数)
Newton-Cotes法 :台形則、Simpson 1/3 則、Simpson 3/8 則、Boole則
・被計算式 Chebyshev多項式
(1) 3 次式、 積分域 [2,6]、 解析 1232
T
3(x)=4x
3-3x
(2) 12 次式、 積分域 [-0.15,1.85]、 解析 143632.865175008
T
12(x)=2048x
12-6144x
10+6912x
8-3584x
6+840x
4-72x
2+1
(3) 25 次式、 積分域 [-0.95,1.05]、 解析 14.6430725603351
T
25(x)=16777216x
25-104857600x
23+288358400x
21-458752000x
19+466944000x
17-
-317521920x
15+146227200x
13-45260800x
11+9152000x
9-1144000x
7+80080x
5-
-2600x
3+25x
3、数値積分の手法の評価
上記 (1)、(2)、(3)の積分域 [a,b] を均等に分割し、その区分の中で分点と重み係数を使用し、少ない分割
数で解析値に到達した順より評価した。
分割数の少ない順 (精度高く値を得られる分割数の少ない順で推定した)
順位1、Gauss-Legendre 6分点
順位2、Gauss-Legendre 5分点
順位3、Chebyshev 6分点
順位4、Gauss-Legendre 4分点
順位5、Boole則
順位6、Chebyshev 5分点
順位7、Chebyshev 4分点
順位8、Gauss-Legendre 3分点
順位9、Simpson 3/8則
順位10、Simpson 1/3則
順位11、Chebyshev 3分点
順位12、Chebyshev 2分点、Gauss-Legendre 2分点
順位14、台形則
なお、Chebyshevの2分点と Gauss-Legengreの2分点は、同じ重み、分点であるので、同じ結果を得る。
4、計算結果の一覧
1)Chebyshev多項式 T
3次
2)Chebyshev多項式 T
12次
3)Chebyshev多項式 T
25次
5、使用したSTL
vector<double>fx、inner_product()
具体的な数値計算2(予備)
Information
- 2012年
- この章の統合履歴
2012年1月 2011年の 数値積分と積分方程式の章の一部を統合
- 2012&13年
- この章の変更追加履歴
10)タワー、御苑の画像を追加 2013/9/28
9)防災を環状2号の秋に差替え 2013/9/28
8)トマト、花、Skytree、母子画、市民の画像の追加 2013/9/27
7)画像の追加 2013/05/12
6)Chebyshev積分の4次式の第5項係数の算出式根拠の記載 2012/11/20
5)Chebyshev積分の4次式の第4項係数の算出式根拠の記載 2012/10/26
4)具体的な数値積分計算1の記載 2012/3/25
3)Gauss-Chebyahev法の記載 2012/3/23
2)Legendre法の記載 2012/3/22
1)Newton-Cotes法の記載 2012/3/21