Ｃ言語による数値計算プログラミング演習

Transcription

1 8. 数値積分 任意の区間 [,b] における f() の定積分 b () I = f ( ) d の値は, つぎのように n 点の関数値の和により近似的に与えられる () In = Ak f ( k) n k = このとき, k を分点,A k を重みという 8. ニュートン コーツ ( 複合型 ) 積分公式 積分区間 [,b] を等分割して n 個の分点をとり, 被積分関数 f() を n- 次ラグランジュ補間多項式で近似して得られる積分公式を n 点ニュートン コーツ公式という ただし高次多項式による補間は凹凸の激しい不自然な形状に陥ることがあるので, 実際の適用には積分区間 [,b] をまず 個の小区間に等分し, 各小区間ごとに比較的低次の補間公式を用いることが多い これを複合型積分といい, 代表的なものを以下に記す 複合中点則 : 分点として各小区間の中点をとり, 各小区間ごとに 次式補間を行うもの i + i M = h f i= b h = i = + ih i =,, L, であり, 誤差は O(h ) となる N N N N! N = > 複合台形則 : 分点として小区間の端点をとり, 各小区間ごとに 次式補間 ( 面積の台形近似 ) を行うもの h = ( ( i ) + i ) = i i= + + i= T f f ( ) h f ( ) f ( ) f ( ) b h = i = + ih i =,, L, N N N であり, 誤差は O(h ) となる N! N = > 複合シンプソン則 : 分点として小区間の端点と中点をとり, 各小区間ごとに 次式補間を行うもの h h S = ( f ( i ) + f( i ) + f( i) ) = f( ) f( ) f( i) 4 f( ) i i = i= i= b h = i = + ih i =,, L, であり, 誤差は O(h 4 ) となる また複合シンプソン則は以下のように表現できる S = (/3) M + (/3) T N N N " N $ N = >

2 8. ガウス型積分 () 式の定積分の値を求めるための n 点ガウス ルジャンドル積分公式は n b Gn = wi f( i ) i= + b b i = + ui = = > > > = > = となる ここに u i は分点 ( -< u i < ; n 次ルジャンドル多項式の零点 ) w i は重み ( < w i <= ) であるが, それらの値は通常, 数表に掲載されている 分点 u,u,,u n は原点に対称に配置されており (n が奇数のときは u (n+)/ =), 対称な点に対する重みは等しい n 点ガウス型公式は, 被積分関数 f() が n- 次以下の多項式であれば厳密な積分値を与える B N N K アルゴリズム 8.. n 点ガウス公式の分点 u,u,,u n と重み w,w,,w n を数表からセットする s:= for i= to n i :=(+b)/+(b-)/u i s:=s+w i f( i ) end s:= (b-)/s 問題 つぎの定積分 (8.4 節参照 ) e cos d=.37846l の値を n 点ガウス公式 (n=3,4,5) により求める プログラム /* Gussin qudrture */ #include <stdio.h> #include <th.h> flot gu[]= ,., ,.33998,.86363,., ,.96798,.3869,.6694, ,.,.45845,.7453,.94979, ,.55534, ,.96899,.,.34534,.63374,.8363,.9686, , , , ,

3 ; flot gw[]=., , ,.6545, , , ,.36969, ,.36766,.7345,.47959,.3883,.79754,.9485, ,.33766,.38,.85,.33394,.3347,.667,.8648,.8744,.9554,.69667,.9864,.49453, ; int inde[]=,,,, 3, 5, 8,,5,9,4; flot f(flot ) return( ep()*cos() ); void in(void) int n,,i,k,n,n; flot,b,b,b,s,wi,i,f; flot u[], w[]; printf("order of Gussin qudrture: n? "); scnf("%d", &n); if((n<) (n>)) /* error */ printf("n should be =< n =< n"); eit(); else /* setup Gussin points nd weights */ n= n/; n=(n+)/; for(i=;i<n;i++) u[n+i]= gu[inde[n]+i]; w[n+i]= gw[inde[n]+i]; for(i=;i<n;i++) u[i]=-u[n--i]; w[i]= w[n--i]; printf("intervl of qudrture:, b? "); scnf("%g%g", &, &b); printf(" n t u t t weight t t t f() n"); b=(+b)/; b=(b-)/; s = ; for(i=;i<n;i++) i = b+b*u[i]; f = f(i); 3

4 s += w[i]*f; printf("%5.6e %5.6e %5.6e %5.6e n", u[i], w[i], i, f); s *=b; printf(" nintegrl=%3.7g n", s); 実行例 実行開始 order of Gussin qudrture: n? 3 intervl of qudrture:, b?. 3 点ガウス公式 u weight f() e e-.77e-.97e+.e e- 5.e e e e e e+ integrl= おしまい 実行開始 order of Gussin qudrture: n? 4 intervl of qudrture:, b?. 4 点ガウス公式 u weight f() e e e-.6936e e e e-.3593e e e e-.5377e e e e-.549e+ integrl= おしまい 実行開始 order of Gussin qudrture: n? 5 intervl of qudrture:, b?. 5 点ガウス公式 u weight f() e e- 4.69e e e e e-.675e+.e e- 5.e e e e e-.5548e e e e-.59e+ integrl= おしまい

5 8.3 ロンバーグ積分 区間 [,b] における f() の定積分 () の値を求めるためのロンバーグ積分のアルゴリズムは, 以下のとおり 小区間数を n =( 小区間幅 :(b-)) から始めて倍々にして n k = k- ( 小区間幅 :h k =(/) k- (b-)) としていくとき,k=,, に対し, 収束するまでつぎの Step, Step を繰り返す Step) k 段目 ( 小区間数 :n k, 小区間幅 :h k ) の初期近似として, 複合台形則 T(h k ) により T k, を求める k= のとき h T, = T( h ) = ( f( ) + f( b) ) k のとき n k Tk, = T( hk) = hk f( ) + f( + ihk) + f( b) i= これは,k- 段目の複合台形則 T(h k- ) を用いて効率的に計算することができる n k Tk, = T( hk) = T( hk ) + hk f ( + (i ) hk) Step) k では, 線形補外の反復により =h = における多項式補外値を求める T k, j+ 収束判定条件は下式とする i= j 4 Tk, j Tk, j =, j =,, L, k j 4 T T kk, k, k T kk, < ε アルゴリズム 8.3. n:= h:=b- t :=(h/)(f()+f(b)) k= for k=,3, k=,3, : 収束するまで繰り返す for j= to k- end t old j:=t j end n:=n Step h:=h/ 区間数を倍 ( 区間幅を/) にして複合台形則 T(h) により初期近似を行う t :=T(h) T(h) の計算は, 上記のような効率的な計算に置き換えてもよい r:= for j= to k- Step 線形補外の反復により,h=における多項式補外値を求める r:=4r r=4 j t j+ :=(r t j -t old j)/(r-) t j =T k,j, t old j=t k-,j end if (t k -t old k-)/t k <ε then brek 収束したらループから出る 5

6 問題 8. 節の問題と同じく, つぎの定積分 (8.4 節参照 ) e cos d=.37846l の値をロンバーグ積分により求める プログラム /* Roberg's qudrture */ /* (trpezoid rule + Richrdson etrpoltion) */ #include <stdio.h> #include <th.h> #define KMAX flot,b,h; int n; flot f(flot ) return( ep()*cos() ); flot trpezoid(flot tk) flot su; int i; h /= ; su=; for(i=;i<n;i++) su += f(+h*(*i+)); n *=; return(tk/+h*su); void in(void) flot eps,r,t[kmax],t[kmax]; int k,j; printf("torelence of Roberg's qudrture: eps? "); scnf("%g", &eps); printf("intervl of qudrture:, b? "); scnf("%g%g", &, &b); printf(" n n n=; h=b-; trpezoid t t*** etrpoltion *** n"); 6

7 t[]=h*(f()+f(b))/; printf("%d %3.6e n", n, t[]); for(k=;k<kmax;k++) for(j=;j<k;j++) t[j]=t[j]; t[] = trpezoid(t[]); r =; for(j=;j<k;j++) r *= 4; t[j+] = (r*t[j]-t[j])/(r-); printf("%d", n); for(j=;j<=k;j++) printf(" %3.6e", t[j]); printf(" n"); if( fbs((t[k]-t[k-])/t[k])<eps ) brek; if(k>=kmax) printf("divergent? n"); printf(" nintegrl=%3.7g n", t[k]); 実行例 実行開始 torelence of Roberg's qudrture: eps?.e-6 intervl of qudrture:, b?. n trpezoid *** etrpoltion ***.34347e+.3468e+.3764e e e+.3788e e+.3787e+.3785e+.3785e e+.3784e+.3785e+.3785e+.3785e+ integrl= おしまい 精度の比較検証定積分 I = e cos d を数値積分公式により計算して, その誤差を調べてみよう 被積分関数 f ( ) = e cos の原始関数は F( ) = e (cos+ sin )/ であるから, 定積分の値 ( 最も近い浮動小数点数に丸められた 6 桁の値 ) は, 以下のようになる I = e cos d=

8 いくつかの数値積分公式により上記の定積分の値を算出し, 誤差を比較検討する なお, ここで示す計算結果はエクセル ( 精度 5 桁?) を用いたものである 複合中点則, 複合台形則, 複合シンプソン則により数値積分を行ったときの誤差 E と分点区間数 n との関係を左図に, 誤差 E と分点間隔 h との関係を右図に両対数グラフで示す 主要な誤差項が E=O(h p ) Ch p である場合,log をとると log E log C +p(log h) となり, 誤差 E は勾配 p の直線上にほぼ並ぶはずである 複合中点則と複合台形則の誤差は O(h ), シンプソン則の誤差は O(h 4 ) であり, 右図中で計算値の誤差はそれぞれ勾配 と勾配 4 の直線上にほぼ並んでいることが確認される また複合シンプソン則の誤差は他の 者に比べて小さいことがグラフに現れている 誤差 E.E+.E-.E- 複合中点則 複合台形則 複合シンプソン則 誤差 E.E+.E-.E- 複合中点則 複合台形則 複合シンプソン則.E-3.E-3.E-4.E-4.E-5.E-5.E-6.E-6.E-7.E-7.E-8 分点区間数 n.e-8.. 分点間隔 h n 点ガウス ルジャンドル積分公式の誤差を分点数 n を横軸にとって左図に示す また, ロンバーグ積分 T,T,T 33,T 44,T 55 の誤差 ( それぞれ O(h ),O(h 4 ),O(h 3 6 ),O(h 4 8 ),O(h 5 ), h k =(/) k- ) を分点区間数 n を横軸にとって上述の複合積分則とともに右図に示す T は 区間の台形則,T は 区間 ( 分点区間数 ) のシンプソン則と同じであるので, 誤差も同様であることが確認できる 以上の図より, ガウス積分の精度が一番良いことが確認される 誤差 E.E-.E-3.E-4.E-5.E-6.E-7 ガウス積分 誤差 E.E+.E-.E-.E-3.E-4 複合中点則複合台形則複合シンプソン則ロンバーグ積分.E-8.E-5.E-9.E-6.E-.E-7.E-.E-8.E-.E-3.E-4.E 分点数 n.e-9.e-.e-.e- 参考文献.E-3 分点区間数 n 葉子 : 数値計算の基礎 解法と誤差, コロナ社 (7) 森口繁一, 伊理正夫, 武市正人編 :C による算法痛論, 東京大学出版会 () 森口繁一 : 数値計算工学, 岩波書店 (989) Heth, Michel T.: Scientific Coputing, An Introductory Survey, McGrw-Hill() 8