積分をPythonで理解する

積分もやはり、Sympy というライブラリが非常に便利です。

Sympy を使った方法

import sympy as sym
from sympy.plotting import plot
sym.init_printing(use_unicode=True)
%matplotlib inline

指定した文字をシンボル（変数を表す文字）として扱います。

a, b, c, x, y = sym.symbols("a b c x y")

主要な数学関数を sympy からインポートしておきます。

from sympy import sin, cos, tan, log, exp

積分

# 積分の式
expr  = x ** a
integ = sym.Integral(expr, x)
print(integ)
integ

Integral(x**a, x)

$$\int x^{a}, dx$$

# 積分の実行
integ.doit()

$$\begin{cases} \frac{x^{a + 1}}{a + 1} & \text{for}: a \neq -1 \\log{\left (x \right )} & \text{otherwise} \end{cases}$$

# このように書いても同じ
sym.integrate(expr, x)

$$\begin{cases} \frac{x^{a + 1}}{a + 1} & \text{for}: a \neq -1 \\log{\left (x \right )} & \text{otherwise} \end{cases}$$

# sym.Eq を使って等式として表記
eq = sym.Eq(sym.Integral(expr, x), sym.integrate(expr, x))
print(eq)
eq

Eq(Integral(x**a, x), Piecewise((x**(a + 1)/(a + 1), Ne(a, -1)), (log(x), True)))

$$\int x^{a}, dx = \begin{cases} \frac{x^{a + 1}}{a + 1} & \text{for}: a \neq -1 \\log{\left (x \right )} & \text{otherwise} \end{cases}$$

いろんな積分の公式

expr = 1/x
eq = sym.Eq(sym.Integral(expr, x), sym.integrate(expr, x))
print(eq)
eq

Eq(Integral(1/x, x), log(x))

$$\int \frac{1}{x}, dx = \log{\left (x \right )}$$

expr = a ** x
eq = sym.Eq(sym.Integral(expr, x), sym.integrate(expr, x))
print(eq)
eq

Eq(Integral(a**x, x), Piecewise((a**x/log(a), Ne(log(a), 0)), (x, True)))

$$\int a^{x}, dx = \begin{cases} \frac{a^{x}}{\log{\left (a \right )}} & \text{for}: \log{\left (a \right )} \neq 0 \x & \text{otherwise} \end{cases}$$

expr  = exp(x)
eq = sym.Eq(sym.Integral(expr, x), sym.integrate(expr, x))
print(eq)
eq

Eq(Integral(exp(x), x), exp(x))

$$\int e^{x}, dx = e^{x}$$

expr = sin(x)
eq = sym.Eq(sym.Integral(expr, x), sym.integrate(expr, x))
print(eq)
eq

Eq(Integral(sin(x), x), -cos(x))

$$\int \sin{\left (x \right )}, dx = - \cos{\left (x \right )}$$

expr = cos(x)
eq = sym.Eq(sym.Integral(expr, x), sym.integrate(expr, x))
print(eq)
eq

Eq(Integral(cos(x), x), sin(x))

$$\int \cos{\left (x \right )}, dx = \sin{\left (x \right )}$$

expr = tan(x)
eq = sym.Eq(sym.Integral(expr, x), sym.integrate(expr, x))
print(eq)
eq

Eq(Integral(tan(x), x), -log(cos(x)))

$$\int \tan{\left (x \right )}, dx = - \log{\left (\cos{\left (x \right )} \right )}$$

expr = log(x)
eq = sym.Eq(sym.Integral(expr, x), sym.integrate(expr, x))
print(eq)
eq

Eq(Integral(log(x), x), x*log(x) - x)

$$\int \log{\left (x \right )}, dx = x \log{\left (x \right )} - x$$

expr = 1/cos(x)**2
eq = sym.Eq(sym.Integral(expr, x), sym.integrate(expr, x))
print(eq)
eq

Eq(Integral(cos(x)**(-2), x), sin(x)/cos(x))

$$\int \frac{1}{\cos^{2}{\left (x \right )}}, dx = \frac{\sin{\left (x \right )}}{\cos{\left (x \right )}}$$

expr = 1/sin(x)**2
eq = sym.Eq(sym.Integral(expr, x), sym.integrate(expr, x))
print(eq)
eq

Eq(Integral(sin(x)**(-2), x), -cos(x)/sin(x))

$$\int \frac{1}{\sin^{2}{\left (x \right )}}, dx = - \frac{\cos{\left (x \right )}}{\sin{\left (x \right )}}$$

expr = (x-a)**b
eq = sym.Eq(sym.Integral(expr, x), sym.integrate(expr, x))
print(eq)
eq

Eq(Integral((-a + x)**b, x), Piecewise(((-a + x)**(b + 1)/(b + 1), Ne(b, -1)), (log(-a + x), True)))

$$\int \left(- a + x\right)^{b}, dx = \begin{cases} \frac{\left(- a + x\right)^{b + 1}}{b + 1} & \text{for}: b \neq -1 \\log{\left (- a + x \right )} & \text{otherwise} \end{cases}$$

expr = 1/sin(x)
eq = sym.Eq(sym.Integral(expr, x), sym.integrate(expr, x))
print(eq)
eq

Eq(Integral(1/sin(x), x), log(cos(x) - 1)/2 - log(cos(x) + 1)/2)

$$\int \frac{1}{\sin{\left (x \right )}}, dx = \frac{\log{\left (\cos{\left (x \right )} - 1 \right )}}{2} - \frac{\log{\left (\cos{\left (x \right )} + 1 \right )}}{2}$$

定積分

expr = (x-a)*(b-x)
eq = sym.Eq(sym.Integral(expr, (x, a, b)), sym.integrate(expr, (x, a, b))).factor()
print(eq)
eq

Eq(-Integral((-a + x)*(-b + x), (x, a, b)), -(a - b)**3/6)

$$- \int_{a}^{b} \left(- a + x\right) \left(- b + x\right), dx = - \frac{\left(a - b\right)^{3}}{6}$$

expr = x/(x**2 + 1)
eq = sym.Eq(sym.Integral(expr, (x, 1, 2)), sym.integrate(expr, (x, 1, 2)))
print(eq)
eq

Eq(Integral(x/(x**2 + 1), (x, 1, 2)), -log(2)/2 + log(5)/2)

$$\int_{1}^{2} \frac{x}{x^{2} + 1}, dx = - \frac{\log{\left (2 \right )}}{2} + \frac{\log{\left (5 \right )}}{2}$$

Sympy ではこれ以上展開してくれる（小数表示にしてくれる）方法が見つからなかったので、とりあえず numpyなどから関数を借りてきて次のように計算してみました。

from numpy import log
-log(2)/2 + log(5)/2

$$0.4581453659370775$$

円周率

円周率は次のような積分で表されます。

expr = 4/(x**2 + 1)
eq = sym.Eq(sym.Integral(expr, (x, 0, 1)), sym.integrate(expr, (x, 0, 1)))
print(eq)
eq

Eq(Integral(4/(x**2 + 1), (x, 0, 1)), pi)

$$\int_{0}^{1} \frac{4}{x^{2} + 1}, dx = \pi$$

Scipy.integrate を使った方法

from scipy import integrate

積分したい式を関数で用意します。

func = lambda x: 4 / (x ** 2 + 1)

integrate.quad で、積分と、推定誤差が得られます。

integrate.quad(func, 0, 1)

$$\left ( 3.1415926535897936, \quad 3.4878684980086326e-14\right )$$

台形則とシンプソン則

台形則やシンプソン則により積分するには、あらかじめ刻み幅を指定しておく必要があります。

import numpy as np
x=np.linspace(0,1,5)

integrate.trapz(func(x), x) # 台形則

$$3.131176470588236$$

integrate.simps(func(x), x) # シンプソン則

$$3.1415686274509804$$

課題

• 積分を数値的に求める区分求積法について、長方形近似、台形則、シンプソン則の違いを説明してください。

• 上記の Integral(x/(x**2 + 1), (x, 1, 2)) の式が解析的に解けるのは上に示した通りですが、これを数値的に解く以下の３つのプログラムを Python で作成してください。ただし、Sympy も Scipy.integrate も使わないこと。インターネットで検索して見つけたプログラムを適宜改変して使っても良いですが、その場合は出典のURLを明記してください。

  • 長方形近似で解くプログラム
  • 台形則で解くプログラム
  • シンプソン則で解くプログラム