風と翼 Windandwing は、c++の数値計算から風と翼の計算へ を目的とした数値計算を行っています。

Addres http://www5.synapse.ne.jp/windandwing/
c++ 数値計算 6  数値積分(連続データ) : Newton-Cotes、                         Legendre、Gauss-Chebyshev

第8章 数値積分(連続データ): Newton-CotesLegendreGauss-Chebyshev

 連続関数に、数値積分 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( x) に異なる重み係数を設けて行う。 Gauss-Chebyshev 法はその次数毎に分点は非等間隔である
が、 f(xj ) の重み係数は同じである。

 なお、Legendre 法は、前章の数値計算 5 直交関数に記載した Legendre 多項式から導きだされるが、Gauss-
Chebyshev 法は前章の Chebyshev 多項式とは関係がない。
   
 その外、無限大領域までを扱う計算方法がある。ここでは触れない。

 これらの数値計算法により、Chebyshev多項式の3次、12次、25次について、数値積分を行い優劣を比較した。
これを、具体的な数値積分計算 1 に示す。

        

Newton-Cotes法

 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.節にシンプソンの公式の説明があり、等間隔の場合について記述されている。)



Legendre法

1、計算の概要
  この方法は、積分域 [a,b] を [-1,1] に変換し、Legendre 多項式の解を分点 xk として、分点 xk に対応し
 た重み係数を用いる。
   
                              8-L-1 式

 上記の f(xj ) の算出は、[-1,1] での Legendre 多項式 Pn の解 xk に領域 [a,b]から変換した xj を 関数 f(xj )
 に代入し、行う。

                                  8-L-2 式

  重み係数 w は Legendre 多項式の解である分点 xk に対応した値で、次の式から算出される。この値は、
 積分域 [a,b] には関係がなく、xk のみによる。従い、Legendre 多項式 Pn の次数毎に、同じ解と同じ重み係
 数を持つ事になる。一旦、それぞれの値を知れば、被積分対象の関数に関係なく使用できる。

                              8-L-3 式

 あるいは、

                                    8-L-4 式

  上記の計算が可能となる理論的な考えは、多くのweb上でも公開されている。その概要は、f(x) を 2n-1 次
 の多項式とし、Pn と Qn-1 の積と剰余式 R(x) で展開できるとする。Pn(x) の全ての解で成り立つラグランジ
 ェ多項式で、f(x) を近似できるので、直交性をもつ Legendre 多項式の解から分点と重み係数を得る事が出
 来るという考えである。Legendre 多項式の直交性は前章の数値計算 5 直交関数に示す。

  これらの考えでラグランジェ多項式から Legendre 多項式への転換とその逆が跳躍点である。


2、分点 xk と重み係数 wj 算出
  分点 xk はLegendre多項式の解である。Legendre 多項式 Pn の具体的な式は、数値計算 5 直交関数
 34 次まで示している。ご参照していただきたい。
  
  6次までのPnは、以下である。
    
    2次
                                    8-L-5 式

    3次
                                     8-L-6 式

    4次
                               8-L-7 式

    5次
                               8-L-8 式

    6次
                          8-L-9 式


  上記の式の解、すなわち分点は下表の通りである。

           

  なお、上記の6次 P6 は解に以下を持つ。

                          8-L-10 式

    しかし、実用的でないので、0.932469514203152とした。

  上述の解に対応した重み係数を求めるには、導関数が必要となる。具体的には、4次の場合、P4 の導関数は次
 式となる。

                                   8-L-11 式

  従い、8-L-3 式は、以下の式で示される。

                             8-L-12 式

  4次の解 xk=sqrt((15.+2*sqrt(30.))/35) の時は、以下で算出できる。

            8-L-13 式

  上記のようにして得た Pn の6次までの重み係数を示す。

       


   また、[-1,1] の積分域の中央部の分点の重み係数は他の分点の場合より大きい。Newton-Cotes法の場合も
  同じである。

   なお、xk の領域 [a,b] への変換値 xj の具体的な算出については、具体的な数値積分計算1に示す。





Gauss-Chebyshev法

1、計算の概要
   このChebyshev法は、非等間隔の分点 xj を持つが、重み係数 wj がその次数で同一である特長を持つ。
  f(xj) と wj の積の回数が少ない利点を持つ。また、Chebyshev多項式(前章に記載)とは関係がない。

  この計算の考えは、
  ・f(x) を [-1,1] で積分した値が、f(x1),f(x2),... f(xn) の和に定数 k を掛けた値にほぼ等しいと仮定する。
    
                              8-C-1式
     
    ただし、上式で R(x) は残余とする。


  ・この f(x) を (n+1) 階まで導関数を持つとし、マクローリン展開を行う。

                  8-C-2式


  ・8-C-1式の右辺のΣの部分を求める。上式の左辺 f(x) の x に x1,x2,.. xn を代入すると以下となる。

                  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 条件を満足する xj を求める事が出来れば、f(x) の積分を、下式のように、近似できる事になる。

                            8-C-15式

  すなわち、xj を含んだ係数を持つ式として

   2 次式 ( 2分点 )は以下の式となる。

                                   8-C-16式

                                8-C-17式


   3次式(3分点)は以下となる。

                               8-C-18式

              8-C-19式


   n 次式(n分点)は以下で、それぞれの係数に xj を含んでいる式となる。

                                   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積分の重み係数と分点を下表に示す。

       


   なお、G6(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 項の係数は、以下となる。

       (x1+x2)=0                                  8-C-31式

              ∵8-C-13 式より (x1+x2)=0                 8-C-32式

    第 2 項の係数は、以下となる。

       2・x1・x2 = (x1+x2)2- (x12+x22) = 0 - 2/3 = - 2/3              8-C-33式

              ∵8-C-14式より  x12+x22 = 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 項の係数は、

       (x1+x2+x3)=0                                8-C-38式

              ∵8-C-13式より (x1+x2+x3)=0

    第 3 項の係数は、

       2・(x1・x2+x2・x3+x3・x1)= (x1+x2+x3)2- (x12+x22+x32) = 0 -1 = -1
                                              8-C-39式

    第 4 項の係数は、

       -3・x1・x2・x3=

            (x1+x2+x3)3- (x13+x23+x33) - 3 (x1+x2+x3)・(x1・x2+x2・x3+x3・x1)

                =0 - 0- 3・0・(-1/2)= 0                   8-C-40式

              ∵8-C-13式より x13+x23+x33=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 項の係数は、

       (x1+x2+x3+x4)=0   (∵8-C-13式より 0 )                  8-C-44式

    第 3 項の係数は、

       2・'(x1・x2+x2・x3+x3・x4+x4・x1)=(x1+x2+x3+x4)2- (x12+x22+x32+x42)

                           = 0 - 4/3 = - 4/3            8-C-45式

              ∵8-C-14式より x12+x22+x32+x42 = 2・2/1/(2*1+1) = 4/3
                                              8-C-46式

    第 4 項の係数は、

       -3(x1・x2・x3+x2・x3・x4+x3・x4・x1+x4・x1・x2)

        =(x1+x2+x3+x4)3

          - (x13+x23+x33+x43)

          -3 (x1+x2+x3+x4)・(x1・x2+x2・x3+x3・x4+x4・x1+x4・x2+x1・x3)

         = 0
          - 0
         - 0

         =0                               8-C-47 式

              ∵8-C-13 式より x1+x2+x3+x4=0

                       x13+x23+x33+x43=0


    第 5 項の係数は、

        -12・x1・x2・x3・x4

         =(x1+x2+x3+x4)4

           + 3(x14+x24+x34+x44)

           - 6(x1・x2+x2・x3+x3・x4+x4・x1+x4・x2+x1・x3)2

           - 4(x13+x23+x33+x43)(x1+x2+x3+x4)

          = 0
          +3・(4/5)
          -6・(-2/3)・(-2/3)
          -4・0

         =-4/15

            ∵8-C-14式より  x14+x24+x34+x44=2・2/1/(2*2+1)=4/5

             また、第3項の計算結果を上式の三番目の項に代入する。

         ∴ x1・x2・x3・x4=1/45


    従い、
                             8-C-48 式



  ・5次式( 8-C-25 式の算出根拠)
    この式の各項の係数の算出過程と値は、改訂前は示したが、今は示さない。後日に示す。
    現在、数式の貼り付けがいつのまにか出来なくなっているので、調査中 2014/2/7

  ・6次式( 8-C-26 式の算出根拠)
    計算続行中


6、その他
  積分域 [a,b] への [-1,1] 域の分点の変換は具体的な数値積分計算の個所で示す。


参考図書:数値計算ハンドブック
     新数値計算法 小国力著、サイエンス社出版





具体的な数値積分計算1(被対象 Chebyshev 多項式)

1、数値積分計算 (コード1)
 数値積分の手法には、前述のように、各種ある。どの方法が一般的な関数に対して効果的であるかを知るために、具体的に計算を行った。
 
 連続関数である Chebyshev多項式に、数値積分 Newton-Cotes 法Legendre 法Gauss-Chebyshev 法適用し、解析値との比較で、手法の優劣比較を行う。

 結論を言うと、計算精度の良い順は、Legendre 法、Gauss-Chebyshev 法、Newton-Cotes法の順である。 

2、数値積分計算の説明
 1) [-1,1] 内の点から積分域[a,b] への点への変換
   この積分域の変換は、8-L-2 式の xj = (b+a)/2 + (b-a)xk/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

      T3(x)=4x3-3x
     
    (2) 12 次式、 積分域 [-0.15,1.85]、  解析 143632.865175008

      T12(x)=2048x12-6144x10+6912x8-3584x6+840x4-72x2+1

    (3) 25 次式、 積分域 [-0.95,1.05]、 解析 14.6430725603351

      T25(x)=16777216x25-104857600x23+288358400x21-458752000x19+466944000x17-
          -317521920x15+146227200x13-45260800x11+9152000x9-1144000x7+80080x5-
          -2600x3+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多項式 T3

        


   2)Chebyshev多項式 T12

        


   3)Chebyshev多項式 T25

        


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

ナビゲーション


環状2号線 シャネルの女性

 横浜 上大岡 夜のmiokaビル-0 

 横浜 上大岡 夜のmiokaビル-1 

 百合の花 

 シンガポール空港 

 横浜 みなとみらい 日本丸 

 横浜 みなとみらい 日本丸 

 横浜 みなとみらい ぷかり桟橋 

 花と花瓶 

 いけばな 

 花 カサブランカ 

 横浜 秋の環状2号線 

 横浜 上大岡イルミネーション 2題 

 横浜 上大岡イルミネーション 2題 

 2013年 元旦 初日の出 

 2013年 元旦の朝の月 

 2013年 正月の夕刻 上大岡駅前 

 2013年 正月の生け花とお飾り 

 2013年正月 鎌倉 小町通り

 2013年正月 鎌倉八幡宮 参道石の階段

 2013年正月 鎌倉八幡宮 大菩薩様
 

2013年正月 鎌倉八幡宮参道でのShot
 

2013年正月 鎌倉八幡宮 参道の人々
 

2013年正月 八幡宮 牡丹園 三花 島錦 
 

2013年正月 八幡宮 牡丹園 三花 明皇宝 
 

2013年正月 牡丹園 三花 八千代椿 
 

Rosetta Stone の文字
 

戸塚駅ビル内 TOtsukana
 

戸塚駅 駅前の階段と交差点
 

横浜駅西口 高島屋を過ぎFedEx通り
 

横浜駅 西口側付近の空
 

橙色のみかん

正月の富士山

さざんかの花

ドイツ フェスタ in 赤レンガ

ドイツ フェスタ in 赤レンガ

サン プリンセス 客船 at 大桟橋 

汽車道 終着駅 (由来説明 下図)  

汽車道 終着駅 由来説明  

赤レンガ倉庫位置(娘さんの左の2本線)
 

ワールドポーターズ (水上バスから)

汽車道 鉄橋(水上バスから)

観覧車(水上バスから)

航海へのサン プリンセス(水上バスから)

象の鼻公園 (水上バス終着桟橋)

象の鼻のソフトクリーム (Rest houseで)

ベイブリッジ (山下公園より)

客船 三隻 大桟橋 (山下公園より)

みなとみらい地区 (山下公園より)

横浜 山下公園前の通り

トマトの実

淡い紫の花

スカイツリー 台東区 西洋美術館より

ラファエル展 西洋美術館

カレーの市民 ロダン作 西洋美術館

東京タワー

新宿御苑 5月