【数値情報処理a】第3回の答え

$x$ の区間 $[a,b]$¹⁾ で関数 $y=f(x)$ と $x$ 軸ではさまれた領域の面積 $S$ は積分（定積分）で求められる。 \[ S=\int_{a}^{b}f(x)dx=[F(x)]_{a}^{b}=F(b)-F(a) \] $x$ を積分変数、 $f(x)$ を被積分関数、 $F(x)$ を $f(x)$ の原始関数という。

$F(x)$ は $x$ で微分すると $f(x)$ になるという関数である。 \[ \frac{dF(x)}{dx}=F^{\prime}(x)=f(x) \] 与えられた $f(x)$ の積分を行うには、上式を満たす $F(x)$ を見つければよい。 $f(x)$ から $F(x)$ をうまく見つけられればよいが、それが難しい場合は数値的に積分を行う方法を用いる。

次の関数を区間 $[a,b]$ で積分を行う。 \[ f(x)=x^{n}\quad(n\neq -1) \]

$f(x)$ の原始関数は \[ F(x)=\frac{1}{n+1}x^{n+1} \] である。

よって定積分は次のようになる。 \[ S=\int_{a}^{b}x^{n}dx=\left[\frac{1}{n+1}x^{n+1}\right]_{a}^{b}=\frac{b^{n+1}-a^{n+1}}{n+1} \]

例えば、$n=1$ の場合 \[ f(x)=x \] 定積分は次のようになる。 \[ S=\int_{a}^{b}x\ dx=\frac{b^{2}-a^{2}}{2}=\frac{(b-a)(b+a)}{2} \] この $n=1$ の場合は、下図のように台形の面積を求めていることになる。 もしくは、大きい直角二等辺三角形から小さい直角二等辺三角形を引いた面積と考えることもできる。

次の関数を区間 $[0,1]$ で積分を行う。 \[ f(x)=4\sqrt{1-x^2} \] この場合、$f(x)$ の原始関数を求めるのは難しいが、積分変数の変数変換 \[ x=\sin\frac{\theta}{2} \] を行うことで積分できる。 \[ S=\int_{0}^{1}4\sqrt{1-x^{2}}dx=\int_{0}^{\pi}(1+\cos\theta)d\theta=[\theta+\sin\theta]_{0}^{\pi}=\pi \] ここで以下の関係を使っている。 \[ \sqrt{1-x^{2}}=\sqrt{1-\sin^{2}\frac{\theta}{2}}=\sqrt{\cos^{2}\frac{\theta}{2}}=\cos\frac{\theta}{2}\quad(0\le\theta\le\pi) \] \[ dx=\frac{d\left(\sin\frac{\theta}{2}\right)}{d\theta}d\theta=\frac{1}{2}\left(\cos\frac{\theta}{2}\right)d\theta \] \[ \cos^{2}\frac{\theta}{2}=\frac{1+\cos\theta}{2} \] \[ \frac{d(\sin\theta)}{d\theta}=\cos\theta \] \[ \sin0=0,\quad\sin\pi=0 \]

この積分は、図形からも求めることができる。

$x$ の区間 $[0,1]$ で $y=\sqrt{1-x^{2}}$ を考えると、これは下図のように半径 $1$ の四分円（円の $1/4$）の円周部である。 これは半径 $1$ の円の式 \[ x^{2}+y^{2}=1^{2} \] から \[ y=\pm\sqrt{1-x^{2}} \] となることより分かる。 つまり、$\sqrt{1-x^{2}}$ の区間 $[0,1]$ の積分は、半径 $1$ の四分円の面積 $\pi/4$ である。

これより、$f(x)=4\sqrt{1-x^{2}}$ の区間 $[0,1]$ の積分は半径 $1$ の円の面積 $\pi$ に等しい。

Moodle Server（非公式）で第4回の小テスト①を受験しなさい。

数値積分では、積分しようとする区間を細かく分割し、分割区間ごとに何らかの方法で近似した面積を求めて足し合わせるというのが基本である。

理想的には無限に細かく分割したいところだが、コンピュータでは限界があるため、有限の分割数にせざるを得ない。 有限の分割数では打ち切り誤差が発生するため、できるだけ細かく分割すれば打ち切り誤差を小さくできそうである。 しかし、細かすぎると丸め誤差などの問題が大きくなるため、あまり細かく分割しなくても精度良く積分できる方法が望ましい。

ここではもっとも簡単な区分求積法、台形公式、シンプソンの公式を取り上げる。 さらに変数変換して精度を上げる方法があるので興味がある人は調べてみてほしい。

区間 $[a,b]$ を非常に小さい間隔で等間隔に分割し、それぞれの小区間の領域を長方形とみなせば、各長方形の面積が計算できる。 その長方形の面積を足し合わせて面積を求める方法を区分求積法という。 これは定積分の定義でもある。

区間 $[a,b]$ において等間隔に $n+1$ 個の分点をとり、$n$ 個の区間に分けるとする。 $x$ の小さい方から $i$ 番目（$i=0,1,2,\cdots,n$）の分点は \[ x_{i}=a+ih \] である。 ここで、$h=(b-a)/n$ は小区間の幅で、$x_{0}=a$ と $x_{n}=b$ は端点である。

$i$ 番目の小区間の長方形の面積は \[ f(x_{i})h \] となる。 これらの全ての長方形の面積を足し合わせると \[ S_{n}=\sum_{i=0}^{n-1}f(x_{i})h \] となる。 数値計算では有限の $n$ でこの和を計算する。

そもそも定積分の定義は、区分求積法で $n\rightarrow\infty$（$h\rightarrow 0$）の極限をとって \[ S=\lim_{n\rightarrow\infty}S_{n} \] である。

次の関数を区間 $[0,1]$ で積分を行う。 \[ f(x)=4\sqrt{1-x^{2}} \] この定積分は例 2 で見たように $\pi$ になる。 この積分を区分求積法で行えば、$\pi$ の近似値が得られることになる。

quadrature.R は区分求積法で $f(x)$ の積分を行うプログラムである。 区間数は $n=1000$ にしている。

quadrature.R

# 区分求積法
 
# 被積分関数の定義
f <- function(x) {
  return(4 * sqrt(1 - x ^ 2))
}
 
# 区間数
n <- 1000
 
# 区間
a <- 0
b <- 1
 
# 幅
h <- (b - a) / n
 
s <- 0
for (i in 0 : (n - 1)) {
  x <- a + i * h
  s <- s + f(x) * h
}
print(sprintf("n = %d", n), quote = FALSE)
print(sprintf("h = %f", h), quote = FALSE)
print(sprintf("S = %f", s), quote = FALSE)

Moodle Server（非公式）で第4回の小テスト②を受験しなさい。

区分求積法では小区間の領域の面積を長方形で計算していたが、これを台形で近似して計算する方法を台形公式という。

区分求積法と同様に区間 $[a,b]$ を $n$ 個の区間に分ける。 $i$ 番目の小区間の台形の面積は \[ \frac{f(x_{i})+f(x_{i+1})}{2}h \] となる。 これらの台形の面積の和をとると \[ S_{n}=\sum_{i=0}^{n-1}\frac{f(x_{i})+f(x_{i+1})}{2}h=\left\{\frac{f(a)}{2}+\sum_{i=1}^{n-1}f(x_{i})+\frac{f(b)}{2}\right\}h \] となる。 これが台形公式である。

小区間の台形の面積は、$f(x)$ を一次関数で近似して計算したものに等しい。 それは次のことから分かる。

区間 $[a,b]$ を等間隔に分けるということは、次のように積分区間を分割することである。 \[ \int_{a}^{b}f(x)dx=\sum_{i=0}^{n-1}\int_{x_{i}}^{x_{i+1}}f(x)dx \] ここで $\int_{x_{i}}^{x_{i+1}}f(x)dx$ は $i$ 番目の小区間の面積である。

小区間で $f(x)$ を一次関数 \[ f(x)\simeq Ax+B \] で近似すると、$i$ 番目の小区間の面積は \begin{eqnarray} \int_{x_{i}}^{x_{i+1}}f(x)dx&\simeq&\int_{x_{i}}^{x_{i+1}}(Ax+B)dx\\ &=&\frac{A}{2}(x_{i+1}^{2}-x_{i}^{2})+B(x_{i+1}-x_{i})\\ &=&\left\{\frac{A}{2}(x_{i+1}+x_{i})+B\right\}(x_{i+1}-x_{i})\\ &=&\frac{Ax_{i}+B+Ax_{i+1}+B}{2}(x_{i+1}-x_{i})\\ &=&\frac{f(x_{i})+f(x_{i+1})}{2}h \end{eqnarray} となる。 これは台形の面積と同じである。

台形公式で次の関数を区間 $[0,1]$ で積分を行う。 \[ f(x)=4\sqrt{1-x^2} \]

trapezoid.R は台形公式で $f(x)$ の積分を行うプログラムである。 区間数は $n=1000$ にしている。

trapezoid.R

# 台形公式
 
# 被積分関数の定義
f <- function(x) {
  return(4 * sqrt(1 - x ^ 2))
}
 
# 区間数
n <- 1000
 
# 区間
a <- 0
b <- 1
 
# 幅
h <- (b - a) / n
 
s <- 0.5 * (f(a) + f(b))
for (i in 1 : (n - 1)) {
  x <- a + i * h
  s <- s + f(x)
}
s <- s * h
print(sprintf("n = %d", n), quote = FALSE)
print(sprintf("h = %f", h), quote = FALSE)
print(sprintf("S = %f", s), quote = FALSE)

台形公式では小区間で $f(x)$ を一次関数で近似して計算したものと同じだったが、これを二次関数で近似して計算するのがシンプソンの公式である。

小区間において $f(x)$ を二次関数 \[ f(x)\simeq Ax^{2}+Bx+C \] で近似すると、$i$ 番目の小区間の面積は積分で \begin{eqnarray*} \int_{x_{i}}^{x_{i+1}}f(x)dx&\simeq&\int_{x_{i}}^{x_{i+1}}(Ax^{2}+Bx+C)dx\\ &=&\frac{A}{3}(x_{i+1}^{3}-x_{i}^{3})+\frac{B}{2}(x_{i+1}^{2}-x_{i}^{2})+C(x_{i+1}-x_{i})\\ &=&\left\{\frac{A}{3}(x_{i}^{2}+x_{i}x_{i+1}+x_{i+1}^{2})+\frac{B}{2}(x_{i}+x_{i+1})+C\right\}(x_{i+1}-x_{i})\\ &=&\frac{h}{6}\left\{f(x_{i})+4f(x_{i}+h/2)+f(x_{i+1})\right\} \end{eqnarray*} となる。 これらの面積の和をとると \begin{eqnarray*} S_{n}&=&\sum_{i=0}^{n-1}\frac{h}{6}\left\{f(x_{i})+4f(x_{i}+h/2)+f(x_{i+1})\right\}\\ &=&\frac{h}{6}\left\{f(a)+4\sum_{i=0}^{n-1}f(x_{i}+h/2)+2\sum_{i=1}^{n-1}f(x_{i})+f(b)\right\} \end{eqnarray*} となる。 これがシンプソンの公式である。

シンプソンの公式で次の関数を区間 $[0,1]$ で積分を行う。 \[ f(x)=4\sqrt{1-x^2} \]

simpson.R はシンプソンの公式で $f(x)$ の積分を行うプログラムである。 区間数は $n=1000$ にしている。

simpson.R

# シンプソンの公式
 
# 被積分関数の定義
f <- function(x) {
  return(4 * sqrt(1 - x ^ 2))
}
 
# 区間数
n <- 1000
 
# 区間
a <- 0
b <- 1
 
# 幅
h <- (b - a) / n
 
s1 <- f(a + h * 0.5)
s2 <- 0
for (i in 1 : (n - 1)) {
  x <- a + i * h
  s1 <- s1 + f(x + h * 0.5)
  s2 <- s2 + f(x)
}
s <- (f(a) + 4 * s1 + 2 * s2 + f(b)) * h / 6
print(sprintf("n = %d", n), quote = FALSE)
print(sprintf("h = %f", h), quote = FALSE)
print(sprintf("S = %f", s), quote = FALSE)

Moodle Server（非公式） で第4回の課題を行いなさい。

締め切り：2025年5月7日（水）20時