積分の答え合わせをしたい
この夏は訳あって数学の院試を一人で解き続けていた。院試の解答はもってないので、基本的に答え合わせはできないのだが、積分の計算問題は数値計算で結果を求めれば検算ができると気がついた。例えば次のような問題。(名古屋大学多元数理2022年度2次募集1日目大問4)
D = \{(x,y)\in \mathbb{R}^2 \mid \sqrt{x} + \sqrt{y} \leq 1, x \geq0, y \geq0 \}とするとき，
\iint_{D}|x-y|dxdyの値を求めよ．
これを解く方法はいろいろあると思うが、
\begin{pmatrix}
x\\y
\end{pmatrix}
=
\begin{pmatrix}
(s+t)^2\\(s-t)^2
\end{pmatrix}
と変形して、積分区間をDと1体1対応するように
D' = \{(s,t)\in \mathbb{R} \mid 2s \leq 1, s+t \geq0, s-t \geq0 \}
被積分関数を
|x-y| = 4s|t|
と変形して積分するのが一つのやり方だろう。計算すると、答えは
\frac{1}{24}
と出た。これの答え合わせをしたい。
素朴なやり方
もともと積分区間は$[0,1]\times[0,1]$の中に入るので、それに対応するN×Nの行列を用意して、その(i, j)成分に座標(i/N, j/N)での被積分関数の値を代入し、積分区間内でその値の和を取れば良い。
import numpy as np
import matplotlib.pyplot as plt
N = 1000
F = np.zeros([N, N])
def domain(x, y):
return x**0.5 + y**0.5 <= 1
def function(x, y):
return abs(x - y)
for i in range(len(F)):
for j in range(len(F[i])):
x, y = i/N, j/N
if domain(x, y):
F[i, j] = function(x, y)
plt.imshow(F.T, cmap="Blues", origin="lower") # 積分区間と被積分関数を可視化してみる。
F.flatten().sum() / (N * N), 1 / 24 # dxdy = 1/(N*N)、結果は1/24に一致するはず
結果
(0.042201526, 0.041666666666666664)
数値計算と手計算の結果は大体同じになった。おそらくあっているだろう。N = 5000にすると、
(0.04176988024000004, 0.041666666666666664)
とさらに精度が上がる。ただし自分のパソコンで5秒ほど計算に時間がかかってしまった。
改良
そういえば、Numpyの行列計算でfor文を使うのはアンチパターンであった。for文は使わずにnp.ndarrayのキャスト機能等を駆使して、行列をあたかも数値であるかのように扱うのがコツである。そうすることで変数や関数を、その定義域上でとりうる値を要素としてもった行列として実現できるのだった。
import numpy as np
import numpy.ma as ma
import matplotlib.pyplot as plt
r = 1 # 領域[0,r]*[0,r]
N = 5000 # 0-1のサンプル点数
index = np.linspace(0, r, int(r * N)+1)
X, Y = np.meshgrid(index, index) # x座標とy座標を持った行列
# 積分区間以外はマスクする
domain = X**0.5 + Y**0.5 <= 1
Y = ma.masked_where(~domain, Y)
X = ma.masked_where(~domain, X)
# 被積分関数に対応する行列
F = abs(X - Y)
plt.imshow(F, cmap="Blues", origin="lower")
F.flatten().sum() / (N * N), 1/24 # 一致するはず
計算時間は5秒→0.5秒に短縮された。
次の問題
次の問題にいこう。積分区間がさっきより複雑である。(名古屋大学多元数理2022年度1日目大問3)
\{(x,y)\in \mathbb{R}^2 \mid x \gt0, y \gt0 \}において，4つの曲線$xy = 1,\ xy = 2,\ y = x,\ y = 3x$
によって囲まれる領域を$D$とする．このとき\iint_{D}\frac{y}{x}dxdyを求めよ．
これは、積分区間が簡単になるように
\begin{pmatrix}
u\\v
\end{pmatrix}
=
\begin{pmatrix}
xy\\y/x
\end{pmatrix}
と変形するのが良いだろう。答えはちょうど1とでた。
この場合も、さっきとほとんど同様にできる。あらかじめ簡単にグラフを書いて、積分区間が$[0,3]\times[0,3]$に入ることを確認する。
import numpy as np
import numpy.ma as ma
import matplotlib.pyplot as plt
r = 3 # 領域[0,r]*[0,r]
N = 1000 # 0-1のサンプル点数
index = np.linspace(0, r, int(r*N)+1)
X, Y = np.meshgrid(index, index)
# 積分区間
d1 = Y >= X
d2 = Y <= 3 * X
d3 = Y*X >= 1
d4 = Y*X <= 2
domain = d1*d2*d3*d4
# 積分区間以外はマスクする
Y = ma.masked_where(~domain, Y)
X = ma.masked_where(~domain, X)
plt.imshow(Y/X, cmap="Blues", origin="lower")
(Y/X).flatten().sum()/(N*N) # 1に一致するはず
結果(0.2秒)
1.0004635842841911
終わりに
やり方がわかったので、他の積分も検算ができそうである。
ちゃんと調べられていないが、scipyを使ったら、行列のことなど考えずに抽象的に高い精度の数値計算が実行できるようだ。自分で作った方が仕組みがわかって面白いような気もする。