読者です 読者をやめる 読者になる 読者になる

naoya_t@hatenablog

いわゆるチラシノウラであります

Pythonで数値積分 〜フーリエ級数展開を例に、てかscipy.integrate.quad()かわいいよscipy.integrate.quad()

はいはいまたグラフ描きたいだけのエントリですよ。てか数値積分が簡単に出来ちゃうSciPy(・∀・)イイネ!!

\displaystyle\int_a^bf(x) dx

from scipy.integrate import quad

y, abserr = quad(f, a, b)

で計算できちゃうのです。(y が積分値、abserrは推定誤差の絶対値)

基礎編

\displaystyle\int_0^2 x^3 dx
\displaystyle[\frac14{x^4}]_0^2=\frac{2^4}4-\frac{0^4}4=4-0=4 を期待。

# -*- coding: utf-8 -*-
from scipy.integrate import quad

y, abserr = quad(lambda x:x*x*x, 0, 2)
print "∫_0^2 x^3 dx = %g ± %g" % (y, abserr)

→ ∫_0^2 x^3 dx = 4 ± 4.44089e-14
正解。推定誤差も小さい。
相手はSciPy様なので人間に正解とか言われてもおそらくあれですね。

\displaystyle\int_0^{\frac\pi3} \sin x dx
\displaystyle[-\cos x]_0^{\frac\pi3}=(-\cos\frac\pi3)-(-\cos0)=(-\frac12)-(-1)=\frac12 を期待。

# -*- coding: utf-8 -*-
from numpy import sin
from scipy.integrate import quad

y, abserr = quad(sin, 0, pi/3)
print "∫_0^{¥pi/3} sin x dx = %g ± %g" % (y, abserr)

→ ∫_0^{\pi/3} sin x dx = 0.5 ± 5.55112e-15
これも正解。

応用編

積分を使ってて簡単そうなところで(かつグラフが描けそうなところで←ここ重要)、フーリエ級数展開とかしてみる。
\displaystyle f(x)=\frac{a_0}2+\sum_{n=1}^{\infty}a_n\cos nx dx+b_n\sin nx dx
の係数の計算に定積分を使うあれです。

  • \displaystyle a_n=\frac1\pi\int_{-\pi}^{\pi}f(x)\cos nx dx (n=0,1,2,...)
  • \displaystyle b_n=\frac1\pi\int_{-\pi}^{\pi}f(x)\sin nx dx (n=1,2,...)
# -*- coding: utf-8 -*-
from pylab import *
from scipy.integrate import quad

def fourier(fun, n_max):
  a = []
  b = []
  for n in xrange(0, 1+n_max):
    res, err = quad(lambda x:fun(x)*cos(n*x), -pi, pi)
    a.append(res/pi)
    res, err = quad(lambda x:fun(x)*sin(n*x), pi, pi)
    b.append(res/pi)
  def fn(x):
    sum = a[0] / 2
    for n in xrange(1, 1+n_max):
      sum += a[n]*cos(n*x) + b[n]*sin(n*x)
    return sum
  return fn

def f(x):
  x = (x + pi) % (pi * 2) - pi
  if x >= 0:
    return 1
  else:
    return 0

x_min = -pi - 0.5
x_max = pi*2 + 0.5
y_min = -0.25
y_max = 1.25
axis([x_min, x_max, y_min, y_max])

f_fn = fourier(f, 10)
xs = linspace(x_min, x_max, 256)
plot(xs, amap(f, xs), 'b:', lw=1)
plot(xs, amap(f_fn, xs), 'r-', lw=1)
show()

f:id:n4_t:20120303211816p:plain

いくつか遊んでみた結果をアニメGIFにしてみた。ImageMagick

$ convert -delay 42 f1_*.png -loop 0 anime_f1.gif

みたいに。





Gistに貼ったソースをここにembedしてみるテスト。(2/27から出来るようになったらしいので)
ライセンスはNYSLでどぞ。