とらりもん - python入門: 17. 数値積分 Diff

!数値積分

次に積分を考えよう。

積分とは、関数に微小量をかけて足し合わせることである。すなわち, 関数~~f~~(~~x~~)を~~a~~から~~b~~まで積分するとは,

~~S~~ = ~~f~~(~~t~~__1__)(~~x~~__1__-~~x~~__0__) + ~~f~~(~~t~~__2__)(~~x~~__2__-~~x~~__1__) +  ・・・ + ~~f~~(~~t~~__~~n~~-1__)(~~x~~__~~n~~__-~~x~~__~~n~~-1__)

を求めることである。ここで, ~~x~~__0__=~~a~~, ~~x~~__n__=~~b~~であり, ~~t~~__~~k~~__は~~x~~__~~k~~-1__以上~~x~~__~~k~~__以下の数である。数学的には、~~x~~__1__ - ~~x~~__0__等を限りなく0に近づけるという概念操作を行うが、数値的にはそんなことは気にしない。かわりに~~x~~の間隔をじゅうぶんに狭くとっておけばよろしい。その間隔をどこでも等しく設定して, それを~~dx~~と書くと,

~~S~~ = ~~f~~(~~t~~__1__)~~dx~~ + ~~f~~(~~t~~__2__)~~dx~~ +  ・・・ + ~~f~~(~~t~~__~~n~~-1__)~~dx~~

すなわち,

~~S~~ = {~~f~~(~~t~~__1__) + ~~f~~(~~t~~__2__) +  ・・・ + ~~f~~(~~t~~__~~n~~-1__)}~~dx~~

となる。さらに, ~~t~~__~~k~~__として~~x~~__~~k~~-1__を採用することにすれば,

~~S~~ = {~~f~~(~~x~~__0__) + ~~f~~(~~x~~__1__) +  ・・・ + ~~f~~(~~x~~__~~n~~-1__)}~~dx~~

となる。つまり, xの各値での関数値を合計して, ~~dx~~をかければよい。要するに数列の和にdxをかけるだけだ（本当はもっと巧妙で高精度なやり方があるが, この教材ではこの程度で十分である）。

例として, ~~f~~(~~x~~)=exp(-~~x~~^2)~~f~~(~~x~~)=exp(-~~x~~""2"")を, -7から7まで積分してみよう。~~x~~の微小量, つまり~~dx~~は0.1としよう。なお、この関数はガウス関数といって、統計学の正規分布など、広い範囲で応用される。この関数の不定積分は解析的に行うと非常に難しいのだが, -∞から∞までの定積分は解析的に実行できて、パイの平方根（1.7724...）となることがわかっている（標準正規分布の式で係数にパイの平方根が現れるのはそのためである）。ここでも我々は-∞から+∞まで積分したいのだが, 計算機は∞を扱えないので, かわりに-7から7までにしておく。それでも十分な精度が出るのである。

まず, いつものように, 数学計算ライブラリのnumpyと, グラフを描くためのモジュールmatplotlib.pyplotを用意する:
import numpy as np                                                      
import matplotlib.pyplot as plt

次にxを配列（numpy.ndarrayのインスタンス）として用意する:
x=np.arange(-7, 7, 0.1)
これではxは-7以上6.9以下となり, 7にはわずかに届かないが, まあ気にしないでおこう。

そして, xの各値に対応する関数値を求めるには,
np.exp(-x*x)

とやればよい。積分は, これを全部足してdxつまり0.1をかければよい:
np.exp(-x*x).sum()*0.1

すると,
1.7724538509055225
という値が出てくる。これは√πにほぼ等しいはずである。確認してみよう:
print(np.sqrt(np.pi))
すると,
1.7724538509055159
となる。末尾のあたりがやや違うが、かなりの高精度で計算できたことがわかるだろう。

課題17-1: 以上の操作を再現せよ。

ここまでは, 関数を決められたひとつの区間内で積分すること（定積分）を扱った。

こんどは, 原始関数を求める積分（不定積分）もやってみよう。すなわち, この積分を, -7から7までではなく, -7から~~x~~までとし, ~~x~~を徐々に変えていって, それぞれに応じた積分結果を求めるのである。それにはリスト内包表記を使えばよい:

S=[np.exp(-x[:k]*x[:k]).sum()*0.1 for k in range(x.size-1)]

このSというリストに, 原始関数（~~x~~の各値までの積分結果）が入っている。ここでSの長さを見てみよう:

len(S)

すると139になる。一方で, xの長さは,

x.size

と打てば140と出てくる。つまりSの方がひとつ少ない。これは, xの最後の項に対応する値が無いためである（理由は考えてみよう）。

では, 原始関数をもとの関数と一緒にグラフに描いてみよう:

plt.plot(x, np.exp(-x*x))
plt.plot(x[:-1], S)
plt.show()

課題17-2: 以上の操作を再現せよ。

課題17-3: 関数cos ~~x~~を, 0から7までの~~x~~について積分し（要するに原始関数を~~x~~が0以上7未満の各値について求める）, グラフに描け。解析的に求まる原始関数sin ~~x~~と重ねて描くこと。刻みは0.1とすればよい。