プログラム：数値積分
近ごろは数値積分の計算なんてライブラリーにあるサブルーチンを使えばすむことなのですが、やはり原理を知っておくのは意味があると思います。いったいどうやって数値積分が行なわれているのか、それを知っておくと計算の精度を上げる事や計算時間の短縮のためには何をすれば良いかとかという事も自ずと見えてくるでしょう。そんなわけで今日は自分で積分プログラムを書いてみようという記事です。とはいってもプログラミングについて書くつもりはありません。あくまで原理の話です。
数値積分でもっとも単純な方法は、積分領域をΔｘの短冊に切って、各短冊の面積を長方形近似して足し上げる。つまり
を使う方法です。この式を使ってｎ個に分割した各短冊の面積を足し上げると
∫dx[x0,x1] f(x)
≒ Δx Σi=0, n-1 f(x0 + iΔx) ..........(2)
= Δx ( f(x0) + f(x0+Δx) + f(x0+2Δx) + f(x0+3Δx) + f(x0+4Δx) + .... + f(x0+(n-1)Δx) )
です。積分が面積を求めることだと知っていれば、この計算は最も分りやすい方法です。それは下の図を見れば一目瞭然でしょう。
実際には短冊を無限に小さく刻んでゆく事など不可能ですからこの計算には誤差がついてきます。図をみれば分るとおり、それは短冊とグラフの間にできた隙間の面積です。この隙間の面積は底辺がΔxで高さがおよそ f(x0+Δx)-f(x0) = Δx [f(x0+Δx)-f(x0)]/Δx 〜 Δx f ' (x0)です。つまり誤差は
一つの隙間から来る誤差 〜 (1/2)底辺×高さ〜Δx2
となります（１/2やf ' (x0)は無視しました）。よって短冊切りによる積分の誤差はΔx2*n=(x1-x0)2/nです。nを大きくすれば1/nで誤差は小さくなってゆきます。六ケタの精度を出したいとするなら大体 n=106個の短冊切りをする必要があります。しかしこんな細かい短冊に切って計算するというのは現実的ではありません。仮に一つの短冊の面積を求める際に10-6程度の数値的な丸め誤差があったとして、それが１０6個の短冊和の計算で単純に積ったら全ての努力が水の泡。そもそも積分は塵を集めて山にするプロセスですから、塵の塵を集めて小山くらいにはなる可能性があります。そんなわけで短冊きりをあまり細かくするのは現実的ではありません。誤差を小さくするにすこし工夫が必要です。 取り合えずこの積分ルーチンを foolintegrate[ f[x] , {x,x0,x1} , n]としておきます。このプログラムで計算した具体例を見てみましょう。
foolintegrate[x2, {x, 0, 1}, 20] = 0.33375 (正確な値は1/3=0.33333333... )
foolintegrate[sin(x), {x, 0 , π}, 20] = 1.99589 (正確な値は 2)
20個の短冊に分割すればこの程度の値は出るということですが、逆に言えばこの程度しか精度がでません。
シンプソンによれば一個の短冊の面積を次の公式で求める方が精度が良いことが知られています。
何ですかこのデタラメな公式は？ と思っちゃいますが、言われたとおりにやってみると
simpsonintegrate[ f(x), {x, x0, x1}, n]
= Δx/6 Σi=0, n-1 [ f(x0) + 4 f(x0 + Δx/2) + f(x0+Δx) ] ..........(3)
simpsonintegrate[x2, {x, 0, 1}, 20] = 0.333333
simpsonintegrate[sin(x), {x, 0 , π}, 20] = 2.00001
となって確かに正確な値に近いです。実はシンプソン公式で計算したx2の積分には誤差がありません。0.3333の３はずっと続いてゆきます。（実際には数値を丸めるために当然、どこかで桁がきれてしまいますが。）シンプソンの公式を調べてみましょう。当然ですがデタラメな公式ではないのです(笑）。その方法はテイラー展開です。つまり
f(x) = f(x0+(x-x0))
= f(x0) + f '(x0)(x-x0) + f ''(x0) (x-x0)2/2! + f '''(x0)(x-x0)3/3! + ......
を使えば
∫dx[x0 ,x0+Δx] f(x)
= f(x0)Δx + f '(x0)Δx2/2! + f ''(x0) Δx3/3! + f '''(x0)Δx4/4! + ......
第一項は短冊きりで公式です、それ以降が短冊で開いた隙間を埋めてゆく補正です。現実問題として短冊の幅Δxは十分小さくとりますから補正項は適当なところまで取り込めばほぼ正確な値に近くなります。 さてシンプソンの公式でΔｘに関するテイラー展開をやってみましょう。
Δx/6 [ f(x0)+ 4 f(x0 + Δx/2) + f(x0+Δx) ]
= Δx/6 [ f(x0)]
+ 4Δx/6 [ f(x0) + f ' (x0) Δx/2 + f '' (x0) (Δx/2)2/2! + ..... ]
+ Δx/6 [ f(x0) + f ' (x0) Δx + f '' (x0) Δx2/2! + .....]
= Δx [ f(x0) + f ' (x0) Δx/2! + f '' (x0) Δx2/3! + f ''' (x0) Δx3/4! + .....]
なんとこれはテイラー展開で求めた公式と一致してしまいます。違いはΔx5 の項からです。つまりシンプソンの公式では誤差が小さくなるように作られたうまい公式なわけです。そしてそれはテイラー展開から得られた面積の公式の最初の３項を正確に再現するようになっています。というのもテイラー展開公式の３項は、被積分関数に関する３つの情報 f(x0) , f '(x0) , f ''(x0) を使いますがシンプソンの公式では [x0 , x0+Δx] の領域での f(x0) , f(x0＋Δx/2) , f(x0＋Δx)の三つの情報から微係数３っつを作り出すわけです。つまり何も不思議でも、ましてやデタラメでもなく、ちゃんとした理屈があるわけです。この事からシンプソン公式では正確な値との誤差はΔx5となるため、領域をあまり細かくする必要がありません。これは数値計算上非常に大きなメリットです。
同じような考えの下に作られた積分公式としてガウシアン公式というものがあります。ガウシアン公式は２点での関数値から３次までのテイラー展開を再現するように作られたものです。シンプソン公式の考えさえ理解できればこれも難しくはありません。詳しくはヌメリカルレシピなどを参考にしてください。今回言いたかったことは、美味しいおでんを作りたければ時間をかけて煮込むだけではなく、予め各素材に少し工夫を加えるべしという事です。
数値積分でもっとも単純な方法は、積分領域をΔｘの短冊に切って、各短冊の面積を長方形近似して足し上げる。つまり
| 短冊切り積分公式 ∫dx[x0,x1] f(x) ≒ Δx f(x0) .......(1) |
を使う方法です。この式を使ってｎ個に分割した各短冊の面積を足し上げると
∫dx[x0,x1] f(x)
≒ Δx Σi=0, n-1 f(x0 + iΔx) ..........(2)
= Δx ( f(x0) + f(x0+Δx) + f(x0+2Δx) + f(x0+3Δx) + f(x0+4Δx) + .... + f(x0+(n-1)Δx) )
です。積分が面積を求めることだと知っていれば、この計算は最も分りやすい方法です。それは下の図を見れば一目瞭然でしょう。
実際には短冊を無限に小さく刻んでゆく事など不可能ですからこの計算には誤差がついてきます。図をみれば分るとおり、それは短冊とグラフの間にできた隙間の面積です。この隙間の面積は底辺がΔxで高さがおよそ f(x0+Δx)-f(x0) = Δx [f(x0+Δx)-f(x0)]/Δx 〜 Δx f ' (x0)です。つまり誤差は
一つの隙間から来る誤差 〜 (1/2)底辺×高さ〜Δx2
となります（１/2やf ' (x0)は無視しました）。よって短冊切りによる積分の誤差はΔx2*n=(x1-x0)2/nです。nを大きくすれば1/nで誤差は小さくなってゆきます。六ケタの精度を出したいとするなら大体 n=106個の短冊切りをする必要があります。しかしこんな細かい短冊に切って計算するというのは現実的ではありません。仮に一つの短冊の面積を求める際に10-6程度の数値的な丸め誤差があったとして、それが１０6個の短冊和の計算で単純に積ったら全ての努力が水の泡。そもそも積分は塵を集めて山にするプロセスですから、塵の塵を集めて小山くらいにはなる可能性があります。そんなわけで短冊きりをあまり細かくするのは現実的ではありません。誤差を小さくするにすこし工夫が必要です。 取り合えずこの積分ルーチンを foolintegrate[ f[x] , {x,x0,x1} , n]としておきます。このプログラムで計算した具体例を見てみましょう。
foolintegrate[x2, {x, 0, 1}, 20] = 0.33375 (正確な値は1/3=0.33333333... )
foolintegrate[sin(x), {x, 0 , π}, 20] = 1.99589 (正確な値は 2)
20個の短冊に分割すればこの程度の値は出るということですが、逆に言えばこの程度しか精度がでません。
シンプソンによれば一個の短冊の面積を次の公式で求める方が精度が良いことが知られています。
| シンプソンの積分公式 ∫dx[x0,x1] f(x) ≒ Δx/6 [ f(x0) + 4 f(x0 + Δx/2) + f(x0+Δx) ] |
何ですかこのデタラメな公式は？ と思っちゃいますが、言われたとおりにやってみると
simpsonintegrate[ f(x), {x, x0, x1}, n]
= Δx/6 Σi=0, n-1 [ f(x0) + 4 f(x0 + Δx/2) + f(x0+Δx) ] ..........(3)
simpsonintegrate[x2, {x, 0, 1}, 20] = 0.333333
simpsonintegrate[sin(x), {x, 0 , π}, 20] = 2.00001
となって確かに正確な値に近いです。実はシンプソン公式で計算したx2の積分には誤差がありません。0.3333の３はずっと続いてゆきます。（実際には数値を丸めるために当然、どこかで桁がきれてしまいますが。）シンプソンの公式を調べてみましょう。当然ですがデタラメな公式ではないのです(笑）。その方法はテイラー展開です。つまり
f(x) = f(x0+(x-x0))
= f(x0) + f '(x0)(x-x0) + f ''(x0) (x-x0)2/2! + f '''(x0)(x-x0)3/3! + ......
を使えば
∫dx[x0 ,x0+Δx] f(x)
= f(x0)Δx + f '(x0)Δx2/2! + f ''(x0) Δx3/3! + f '''(x0)Δx4/4! + ......
第一項は短冊きりで公式です、それ以降が短冊で開いた隙間を埋めてゆく補正です。現実問題として短冊の幅Δxは十分小さくとりますから補正項は適当なところまで取り込めばほぼ正確な値に近くなります。 さてシンプソンの公式でΔｘに関するテイラー展開をやってみましょう。
Δx/6 [ f(x0)+ 4 f(x0 + Δx/2) + f(x0+Δx) ]
= Δx/6 [ f(x0)]
+ 4Δx/6 [ f(x0) + f ' (x0) Δx/2 + f '' (x0) (Δx/2)2/2! + ..... ]
+ Δx/6 [ f(x0) + f ' (x0) Δx + f '' (x0) Δx2/2! + .....]
= Δx [ f(x0) + f ' (x0) Δx/2! + f '' (x0) Δx2/3! + f ''' (x0) Δx3/4! + .....]
なんとこれはテイラー展開で求めた公式と一致してしまいます。違いはΔx5 の項からです。つまりシンプソンの公式では誤差が小さくなるように作られたうまい公式なわけです。そしてそれはテイラー展開から得られた面積の公式の最初の３項を正確に再現するようになっています。というのもテイラー展開公式の３項は、被積分関数に関する３つの情報 f(x0) , f '(x0) , f ''(x0) を使いますがシンプソンの公式では [x0 , x0+Δx] の領域での f(x0) , f(x0＋Δx/2) , f(x0＋Δx)の三つの情報から微係数３っつを作り出すわけです。つまり何も不思議でも、ましてやデタラメでもなく、ちゃんとした理屈があるわけです。この事からシンプソン公式では正確な値との誤差はΔx5となるため、領域をあまり細かくする必要がありません。これは数値計算上非常に大きなメリットです。
同じような考えの下に作られた積分公式としてガウシアン公式というものがあります。ガウシアン公式は２点での関数値から３次までのテイラー展開を再現するように作られたものです。シンプソン公式の考えさえ理解できればこれも難しくはありません。詳しくはヌメリカルレシピなどを参考にしてください。今回言いたかったことは、美味しいおでんを作りたければ時間をかけて煮込むだけではなく、予め各素材に少し工夫を加えるべしという事です。
コメント
コメントの投稿
トラックバック
http://letsphysics.blog17.fc2.com/tb.php/241-d64e0387
