Python3 Matplotlib で交点を表示したい：交点近似しよう

Matplotlib はかなり使えるグラフ描画 Python パッケージです。一から自分でプログラムを組んでいくことになるので、確かに自由度は高いのですが、例えば交点を表示するのにも一工夫が必要になります。今回はそんなお話です。

環境
手法 : 近似

環境

> python3 --version
Python 3.8.4

> pip install matplotlib
> pip list
Package         Version
--------------- -------
cycler          0.10.0
kiwisolver      1.2.0
matplotlib      3.2.2
numpy           1.19.0
pip             20.1.1
pyparsing       2.4.7
python-dateutil 2.8.1
setuptools      47.1.0
six             1.15.0

Mac OS X Catalina 10.15 で実行しています。Python 3.8.4 は直インストール、その後で venv 環境を立ち上げてその中で実行しています。

手法 : 近似

近似を使います。曲線Aと曲線Bの交点P付近の、曲線AおよびB上のそれぞれ2点から直線を作り、その2直線の交点P'を交点Pとみなすという手法です。

f:id:Axifof:20200715000257j:plain — 交点Pの近似交点P'を求める

Matplotlib の plot メソッドは点を繋いでいるだけなので、この手法とぴったりマッチします。交点はわからなくても大丈夫です。

さて、ここでは高校数学では典型問題の $y = \sin(x) = \cos(x)$ の交点をいくつか描画してみます。 $\sin(x) - \cos(x) = \sqrt{2} \sin(x - \frac{\pi}{4}) = 0$ より $x = \frac{\pi}{4} + n \pi \,\,(n \in \mathbb{N})$ ですから $n=0, 1, 2$ について計算すれば、交点は

$(0.7854, 0.7071), (3.9270, -0.7071), (7.0686, 0.7071)$

と求めることができます(少数第5桁で四捨五入)。

import numpy as np
import matplotlib.pyplot as plt


class Line():
    def __init__(self, point1, point2):
        dx = point2[0] - point1[0]
        dy = point2[1] - point1[1]
        self.a = dy
        self.b = -dx
        self.c = dx * point1[1] - dy * point1[0]

    def get_intersection(self, line):
        denominator = self.a * line.b - line.a * self.b
        x = (self.b * line.c - line.b * self.c) / denominator
        y = (line.a * self.c - self.a * line.c) / denominator
        return (x, y)


def get_intersections(x, y1, y2):
    idx = np.argwhere(np.diff(np.sign(y1 - y2)) != 0)
    line1 = Line((x[idx], y1[idx]), (x[idx+1], y1[idx+1]))
    line2 = Line((x[idx], y2[idx]), (x[idx+1], y2[idx+1]))
    return line1.get_intersection(line2)


def main():
    x = np.linspace(0, 10, 1000)
    y1 = np.cos(x)
    y2 = np.sin(x)
    plt.xlabel('x')
    plt.xlabel('y')
    plt.xlim(0, 10)
    plt.ylim(-1.5, 1.5)
    plt.grid()
    plt.title(f'intersection test')
    plt.plot(x, y1, color = "red", label='y=cos(x)')
    plt.plot(x, y2, color = "blue", label='y=sin(x)')

    ### 交点を描画する処理はここで行います.
    np.set_printoptions(precision=4, floatmode='fixed')
    px, py = get_intersections(x, y1, y2)
    plt.plot(px, py, 'co', ms=5, color="black", label='interception')
    for x, y in zip(px, py):
        plt.text(x, y, f'({x}, {y})', fontsize=10)

    plt.legend(frameon=False, fontsize=10, numpoints=1)
    plt.show()


if __name__ == '__main__':
    main()

f:id:Axifof:20200714234721p:plain — y=sin(x)とy=cos(x)の交点の描画

いきなり

np.argwhere(np.diff(np.sign(y1 - y2)) == 0)

とすれば、確かに交点を求められそうな雰囲気はあるのですが、これは入力する値が十分に小さくなった場合や、複雑な関数になった場合には誤差の影響で成立しません。ですから近似的に求めておく方が、精度を十分に高めてあげれば役に立ちます。今回のコードで言うと

x = np.linspace(0, 10, 1000)

の "1000" を "10000" などにすればいいわけです(区間[0, 10]を10000等分する)。もちろん処理は重たくなりますが、直線がより厳密になるので値は厳密に近づくことでしょう。

### 交点を描画する処理はここで行います.
np.set_printoptions(precision=4, floatmode='fixed')
px, py = get_intersections(x, y1, y2)
plt.plot(px, py, 'co', ms=5, color="black", label='interception')
for x, y in zip(px, py):
    plt.text(x, y, f'({x}, {y})', fontsize=10)

に関してですが

    plt.text(x, y, f'({x}, {y})', fontsize=10)

で交点の座標を描画しています。"[]" が邪魔ですね。ちなみに

np.set_printoptions(precision=4, floatmode='fixed')

の "percision" パラメータは表示される少数の桁数、"floatmode" パラメータはその表示の方法を決定し、今回の場合だと「少数第4桁まで表示(5桁目で四捨五入)し、桁数不足の場合は0埋めする」という意味です。

よし呟こう

思ったことを呟きます。

Python3 Matplotlib で交点を表示したい：交点近似しよう

環境

手法 : 近似