【python】【至急】2つのグラフの交点をプロットする

Question

実現したいこと

2つのグラフがあり、その交点をプロットしたい。例えば、y=cos xとy=cos(π cos x)の交点を自動的にもとめて、プロットしたいです。
交点の座標は簡単には求められないものです。

前提

y=cos xとy=cos(π cos x)の交点を自動的に算出して、プロットするようなプログラムを書きたい。Pythonで無理なら、別の言語でも大丈夫です。

np.argwhereを使ってみたものの、xの値はnp.linspaceなどで離散的に取得するので、交点となる無理数の解はリストの中に入らないため、プロットされない。

一方で、absなどで丸め誤差を利用するも、2つのグラフが接近する部分で無駄な部分にプロットされる。

そのため、複雑なグラフでは別のアプローチをとるべきなのかなと思いましたが、当方初心者で何もわかりません。調べてみたものの何もヒットしませんでしたので、ここで質問させてください。

今までは、LaTeXでグラフを書いていたのですが、LaTeXでは自動的にグラフの交点を算出するコマンドがありました。しかし、LaTeXでは複雑なグラフを書くことができず、matplotlibを使おうと思って、最近Pythonを学習し始めました。

発生している問題・エラーメッセージ

xの値はnp.linspaceなどで離散的に取得するので、交点となる無理数の解はリストの中に入らない。

absなどで丸め誤差を利用するも、2つのグラフが接近する部分で無駄な部分にプロットされる。

そのため、別のアプローチで交点を求めることが必要？

該当のソースコード

import matplotlib.pyplot as plt
import numpy as np

N = 99
x_min = 0
x_max = 2 * np.pi

def f(x):
    return np.cos(x)

def g(x):
    return np.cos(np.pi * f(x))

x = np.linspace(x_min, x_max, N)
y1 = f(x)
y2 = g(x)

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

plt.plot(x[idx], y1[idx],'.', markersize=8)

plt.plot(x, y1)
plt.plot(x, y2)

plt.show()

試したこと

・Nの数を増やしてみたものの、交点の数は増えない。
・absを使って丸め誤差を利用するも、不要な点がプロットされてしまう。

補足情報（FW/ツールのバージョンなど）

VS codeを使用

Answer 1

「2つのグラフがあり、その交点をプロットしたい。
例えば、y=cos xとy=cos(π cos x)の交点を自動的にもとめて、プロットしたいです。」
の部分だけをChatGPTさんに丸投げしてプログラムを書いてもらいました。

import numpy as np
import matplotlib.pyplot as plt


# Functions
def f1(x):
    return np.cos(x)

def f2(x):
    return np.cos(np.pi * np.cos(x))

# Finding intersection points
x = np.linspace(-2 * np.pi, 2 * np.pi, 1000)
y1 = f1(x)
y2 = f2(x)
intersection_points = np.argwhere(np.diff(np.sign(y1 - y2))).flatten()
x_points = x[intersection_points]
y_points = y1[intersection_points]

# Plotting
plt.plot(x, y1, label="y = cos(x)")
plt.plot(x, y2, label="y = cos(pi*cos(x))")
plt.plot(x_points, y_points, "o", label="Intersection points")
plt.legend()
plt.show()

Comments

[yowakunaran]

ありがとうございます　ChatGPTすごいですね

重解の場合は、条件式使って書けばいいですね！

Answer 2

すでに解決していますが、グラフ描画のために線分の組のデータは作成済みなので、それらが使えそうです。
すなわちxの開始～終了までfとgの各線分の交点を求めればよいかと思います。
この手法だと、Nが小さい（間隔が粗い）場合でも見た目と同じ場所に交点を打てます。

import numpy as np
import matplotlib.pyplot as plt


N = 99
x_min = 0
x_max = 2 * np.pi

def f(x):
    return np.cos(x)

def g(x):
    return np.cos(np.pi * f(x))

# 交点を求める
# ChatGPTの生成したコードなので注意
def find_intersection(line1, line2):
    xdiff = (line1[0][0] - line1[1][0], line2[0][0] - line2[1][0])
    ydiff = (line1[0][1] - line1[1][1], line2[0][1] - line2[1][1])

    def det(a, b):
        return a[0] * b[1] - a[1] * b[0]

    div = det(xdiff, ydiff)
    if div == 0:
        # 線分が平行であるため交点は存在しない
        return None

    d = (det(*line1), det(*line2))
    x = det(d, xdiff) / div
    y = det(d, ydiff) / div

    if not (min(line1[0][0], line1[1][0]) <= x <= max(line1[0][0], line1[1][0])):
        # 交点が線分の範囲外にあるため交点は存在しない
        return None

    if not (min(line1[0][1], line1[1][1]) <= y <= max(line1[0][1], line1[1][1])):
        # 交点が線分の範囲外にあるため交点は存在しない
        return None

    return x, y


x = np.linspace(x_min, x_max, N)
y1 = f(x)
y2 = g(x)

# 各線分の交点を求める
pts = []
for xs, xe, y1s, y1e, y2s, y2e in zip(x[0:-1], x[1:], y1[0:-1], y1[1:], y2[0:-1], y2[1:]):
    t = find_intersection( ((xs,y1s),(xe,y1e)), ((xs,y2s),(xe,y2e)) )
    if t:
        pts.append(t)
pts = np.array(pts)

plt.plot(x, y1)
plt.plot(x, y2)
plt.plot(pts[:,0], pts[:,1],'.', markersize=8)

plt.show()

Comments

[yowakunaran]

これ精度がものすごくたかいですね！
解読してみたいと思います！

Answer 3

既に解決済みですので御参考です。

先に交点の x 座標を求める方針で考えてみました。

• 交点の x 座標（xr）は scipy.optimize.root() に差分の関数 h(x) とグラフから読み取った近似解 [1.2, 2.5, 3.1, 3.8, 5.1] を与えて求める

• 交点の x 座標（xr）に最も近い x 軸配列（xl）の要素のインデックス（ind_xr）は bisect.bisect_left() を利用して求める

  • 「a は単調増加数列」を条件として，a[n-1] < x <= a[n] のとき bisect_left(a, x) は n を返す

  • 次に x が a[n-1] と a[n] のどちらに近いかを調べてインデックスを決定

下記に記述例を示します。
なお，実際の交点とマークのずれを目立たなくするために N = 629 に増やしています。

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import root
from bisect import bisect_left

N = 629
x_min = 0
x_max = 2 * np.pi

def f(x):
    return np.cos(x)

def g(x):
    return np.cos(np.pi * f(x))

def h(x):
    return f(x) - g(x)

xl = np.linspace(x_min, x_max, N)
y1 = f(xl)
y2 = g(xl)

xr = root(h, [1.2, 2.5, 3.1, 3.8, 5.1]).x
print(xr)
## [1.18427679 2.48133238 3.14159265 3.80185292 5.09890852]

ind_xr = np.empty_like(xr, dtype=int)
for i, xi in enumerate(xr):
    n = bisect_left(xl, xi)
    ind_xr[i] = n - 1 if xl[n] - xi > xi - xl[n - 1] else n

np.set_printoptions(suppress=True)
print(xl[ind_xr] - xr)
## [-0.00367827 -0.00007449  0.00000001  0.00007449  0.00367827]

plt.plot(xl, y1)
plt.plot(xl, y2)
plt.plot(xl[ind_xr], y1[ind_xr], '.', markersize=8)

plt.show()

Answer 4

個々の関数(?)のグラフを描画することができるのであれば，画像的にANDを取ればどうですか？
（質問文に示されたグラフで言えば，オレンジ色のグラフと緑色のグラフが重なっている個所を画像処理的に見つけるという話）

単純には，ANDで残った領域（グラフの太さや交差具合次第になるけど，小さな領域になることを期待）の真ん中（重心とか）の位置に点をプロットすればどうでしょうか．
複数の交点があれば複数の領域が得られるので，それらに同じことをすればよいだけです．

……っていうような絵的な結果だけが欲しいだけじゃなくて「交点の座標の値も示す必要があるんだよ」という場合でも，この方法で

• 交点の個数
• それらの座標の荒い値

が得られるのだとすれば，あとはそこから座標を Refinement すれば良いのではないでしょうか．「荒い」とは言え，その値は何らかの Solver に初期値として与えるには十分なんじゃないかな，と．（例えば数ステップのイテレーションで十分に精度改善できそうな気がします．）

Comments

[fana]

グラフの線の太さを細くするとか，グラフ画像の解像度を上げる，とかすることで「荒い座標」の精度はマシになるかな．

[fana]

> 交点の座標は簡単には求められない

ということらしいので，解析解を求めることが望めないなら，数値計算みたいなことをすることになるのだろうから，その実施に必要な情報を得る手段として手っ取り早いと思う話を書いた．

[fana]

> 真ん中（重心とか）

だとちょっと……みたいな場合には，
領域の各画素位置での関数値をちゃんと調べてやれば「グラフ画像のピクセル精度での」最良位置にできるかと．

[yowakunaran]

そういうアプローチもあるんですね。
画像認識の分野はノータッチなので難しいですが、やってみようと思いました！

[fana]

まぁ，「交点」とは字のごとく「交わる場所　＝　両方のグラフが共通して通る場所」ですから，
グラフの絵の上で両方の線が描画されている場所 を見つけてやろう，というだけの話です．

[fana]

※現時点で（何故か）この回答がベストアンサーになってますが，実際に他の方法があなたの問題解決に寄与した場合にはそっちに変えてもらえれば良いかと思います．