［無限遠までの積分］　by ddt³さん：ねこ騙し数学：SSブログ

	ブログをはじめるログイン

［無限遠までの積分］　by ddt³さん　[ddt³さんの部屋] [編集]

［無限遠までの積分］

＞それはそれとして、ネムネコが考えるに、数値計算屋のddt³さんが、きっと、これに関係する記事を投稿してくれると思う・・・

　・・・と言われましても、いったい何を書けばいいんだろう？(^^;)。ネコ先生よりかっちょ良い変数変換でも行えばいいのだろうか？。いや、そうじゃない。「数値計算屋」のご指定があるからには（否定はしないけど(^^;)）、ネコ先生のお望みはきっと台形公式のみを使った強行突破だ！。変数変換なんてもっての外の、とんでもはっぷん！・・・、なのかな？(^^)。

　問題の被積分関数は、

でしたね（図-1）。

　いやぁ～、これはまた性質の良さそうな関数です。

1から滑らかに0を収束する優美な姿(^^)。

でも問題は、(1)を0～∞で積分せよ、というものでした。

　鬱陶しいのは、積分区間の上限が∞である事です。積分区間の上限が∞だと、どこで数値的足し算を打ち切って良いのかを判断しなきゃなりません。自分は有限幅の積分区間の場合、区間幅の100分割を標準にしてます。理由は、100等分して折れ線で結んでやると、大概の関数は十分に滑らかに描けるから。ところで折れ線近似の面積とは、台形公式による結果ですよね？。だから100分割で十分と・・・(^^)。

　図-1をみると0～∞の積分の大勢は、0～10くらいで決まりそうです。なので区間分割幅δは、δ＝10/100＝0.1とします。図-1の白丸は、0.1刻みです。十分ですよね？。

$mseki-fig-001.png$ 　次に積分の打ち切り範囲を考えます。図-1から明らかなように、x＝10から先は、ほぼ傾き0の直線で近似できそうです。そこで10＜xでのyの高さから、その点における接線を引き、x軸と交わる点x＋Lを計算して、図-2に示す三角形の面積を計算します。

なので

から、

三角形の面積をδSとすれば、

　このδSで、xから先の残りの面積を近似してやっても、そんなに間違いではないでしょう。δSが十分小さくなった時、計算を打ち切れば良い訳です。小ささの基準は次のように決めます。

　xの分割幅をδとすると、個々の台形の面積A_jには、A_j×δ²程度の誤差が伴われると考えられます（←ネコ先生がそう言っていた(^^)←それは、（有界な区間の）大域的な誤差で、局所的な誤差はδ³程度ですよ(^^)）。

全体では∑_jA_j×δ²～F(x)･δ²です。ここでF(x)はxまで積分した時の正確な面積。これをF(x)＝∑_jA_jとして台形公式による面積で代用してやると、δ＝0.1だったので、0.01･F(x)。ここにF(x)は、台形公式でxまで積分した結果。

　要するにF(x)の1％くらいの誤差は許容できる、という事になりますので、残りの面積の近似δSについては、F(x)の0.1％くらいになれば十分じゃないでしょうか？。

　つまり、F(x)＝∑_jA_jとδS(x)をモニターしながら台形公式による積算を行い、δS(x)＜F(x)/1000の時点で計算を打ち切る。

　実際に0～10を台形公式で計算すると、F＝1.471なので、この結果をもとに打ち切り工数を概算すると、

から、169.938＜xで、169.938/0.1＝1699～1700ステップくらいになります。

　エッ？、1700ステップは多いって？。そもそも数値計算をやろうって奴が、1000や2000や10000でガタガタ言うんじゃありません。エッ？、計算はExcelでやりたいって？。だからぁ～、現在のExcelは1,048,576行まで使えるんだって・・・。

　エッ？、10,000,000ステップだったら、どうするって？。100万行の列を10列つくりゃぁ～いいだろう！。現在のExcelのフル性能は、1,048,576行×16,384列＝17,179,869,184ステップ分もあるんだ、どうだ参ったか？（← Microsoftの狗(^^;)）。

　なに～っ！、シートを全部ステップで埋めたら、結果を書く場所がないって？。隣のシートを利用せんかい！。もぉ～うるさいなぁ～、100万行も数式をマウスでコピーできないって？。数値計算をやるなら、Excelでマウスは使うな。[Cnt]＋[Shift]＋[上下左右カーソル]のショートカットキーをおぼえろぉ～！。

　・・・もっとも50万行の列が1個あっただけで、重くて嫌になりますけどね(^^;)。で、実際にやってみると1598ステップ，x＝159.8で、F＝1.564539になります。|F－π/2|/(π/2)＝0.004。Fに残りの予想値δS＝0.001564を加えるとF＝1.566103、|F＋δS－π/2|/(π/2)＝0.003です(^^)。収束状況を、図-3に示します。