アトムの物理ノート

積分は物理や工学などでは頻繁に現われる計算です。しかし多くの場合計算が解析的になされることはありません。現実的な場面で出会う積分は複雑すぎるのです。しかしどうしても計算しなくてはならない値がある、そういった場合は数値積分を行ないます。

数値計算のテクニックはコンピュータと共に発展してきたと思われますが、最近は一から全てを自分でプログラムすることは殆どありません。なぜなら様々なライブラリーが既に存在しているからです。又はヌメリカル・レシピなどを参考すればそのアルゴリズムと具体的なプログラミングまで直ぐに学ぶ事ができます（ここでヌメリカル・レシピに関して専門家からの批判的な意見がることを付け加えておきます。それでも私自身はこの本は非常に参考になったと思っています。：http://nakano.webmasters.gr.jp/nr.html。私は数値計算を道具として使った経験があるだけで専門家ではありません。ヌメリカルレシピから多くのことを学びましたがそれがベストな方法か、またはプログラムにバグがなかったかどうかは知りません。ただし大抵２つか３つのプログラムでチェックしてましたから計算にバグはなかったと思いたいです。）。


今回は数値積分について書きたいと思います。例えば

A =∫_0〜∞ sin(x)/x dx

を考えます。この積分に関してFortran(またはCなど）で数値計算を行なう場合の問題点とその解決方法を調べます。この積分は答えが知られていて、mathematicaのような高度な数式処理ソフトなら

A = π/2 = 1.5707963・・・・

という結果を与えてくれます。ところがFortranで積分を数値的に行なう際に良く聞くのが「被積分関数がｘ＝０で発散するからできない」というものです。実際には発散ではないのですが、sin(0)/0 のような不定形はコンピューターでは扱えません。こういった場合に良く使われるのは変数変換で見かけの特異性を除去してしまうもの、または見かけ発散を差っ引いて被積分関数を完全に正則なものに変えてしまう方法です。今回は二番目の方法を説明します。つまりsin(x)/x　の x=0での振る舞いを巧く差っ引いてよい精度の答えを得る方法です。便宜上積分を０〜π　と　π　〜∞の二つの区間に分けて考えます

A1 = ∫_0〜π sin(x)/x dx
A2 = ∫_π〜∞ sin(x)/x dx

---

ところでsin(x)のテイラー展開の公式をしっていれば被積分関数は

sin(x)/x = (x-x³/3!+...)/x　=　1 - x²/3! + ....

なのでx=0でも全く発散はないことに気がつきます。もちろんコンピュータはここまで理解していませんから

sin(0)/0 = [sin(0)の値]/[xのゼロでの値]

という計算をする際に分母に0が入っているので即計算がストップしてしまいます。そこで実際にはx=0ではなく小さな数x₀=10^(-6)を積分の下限にして計算します。

A1'　=　∫_x0〜π sin(x)/x dx

これなら分母にゼロがきませんから計算できます。この時のA1の値の精度はというと

A1' - A1 = ∫_0〜x0sin(x)/x dx　= ∫_0〜x0[ x - x³/3!+...]/x dx　= x₀ - x₀³/(3*3!)+...　〜 x₀

で与えられます。x₀=10^(-6)なら誤差は10^(-6)です。大して大きな誤差ではありませんが、もっと精度が必要な場合もあります。そういった時には原理的にはx₀=10^(-12)などと小さくすれば良いのですが、計算機であまり小さな数を使うのは勧められません。計算時間の増大や丸め誤差の原因になるためです。こんな場合には被積分関数を操作してやって

A1 = ∫_0〜πsin(x)/x dx = ∫_0〜π [sin(x)-x]/x dx + ∫_0〜π x/x dx
= ∫_0〜π [sin(x)-x]/x dx + π

と書きます。二項目は誰でも積分できますからワザワザ数値計算することもありません。被積分関数のxが小さい場合の寄与を差っ引いて、その分は手で計算するわけです。こうすれば被積分関数のx=0での振舞いは元のものよりマイルドになっています。この表式で積分の下限をx₀で置き換えると

A1'' = π　+ ∫_x0〜π [sin(x)-x]/x dx

ですが、その誤差は

A1''- A1 = ∫_0〜x0 [sin(x)-x]/x dx = ∫_0〜x0 [-x³/3!]/x dx = - x³/(3*3!) +・・・・〜　(10^-6)³

となります。

---

それではA2の評価ですがここでもまたもや発散が現われます。数学的な発散ではなく、計算機の問題です。なぜなら

A2 = ∫_π〜∞ sin(x)/x dx

ですが、積分の上限の∞という表式は現実には存在しません。少なくとも計算機で∞を扱う事はできません。よって積分の上限を有限な値で置き換えてやります。

A2' = ∫_π〜nπ sin(x)/x dx

A2積分は被積分関数の無限遠点での振る舞いが1/xですから一見定義されていないように見えますが、分子のsin(x)が正の値と負の値に振動しますから無限遠点での1/xの振る舞い幾分か相殺して収束性のよいものになり、結果として積分は定義されています。よって積分の上限の置き換えは許されますが、誤差を押さえるためには上限nπをどのくらいに大きく取れば良いかが問題になってきます。こういった場合にはnπを色々変えて試して見るのが一番手っ取り早い方法です。nを大きくしていった場合のA2'の値に先に求めたA1''を足してグラフ描きました。横軸がnで縦はA1''+ A2'の値です。赤い実線は比較のための正確な値A = 1.5707963・・・です。ｎを変えた場合の値が上下に分裂している理由は、nが偶数か奇数かによってsin(x)/xの相殺の残りの効果が振動するからです。この図だけからはn=100のオーダーではA=1.57±0.01程度の精度しかでない事が分ります。例えば６桁の精度を出そうと思ったらA２の評価に何か別な方法を用いなければならないことが分ります。