I want to find the 2 closest points on 2 equations
I will be using:
python and sympy
example of my problem:
Say I have 2 equations:
y=x^2y=-(x-3)^2
I would like to find the 2 closest points between these two equations (one point from each equation). How would I do it in sympy?
Note: I want to be able to solve this problem for any nth degree polynomial, such as y=x^4 x^2 x and y=x^3 3x^2etc.
It is ok if I get multiple solutions.
I've searched around and looked at https://math.stackexchange.com/questions/731207/closest-distance-between-two-quadratic-curves#:~:text=First simply take the derivatives,line between the two curves. and a few others, but I don't understand how to implement those solutions in such a way so that it would work for any pair of nth degree polynomial.
My understanding is that I have to do some calculations related to the lines perpendicular to the tangents of the functions, but I have no idea how I can do / implement this in sympy.
If this is impossible in sympy, please point me to some other solution that allows me to find closest points between 2 functions.
TLDR: given 2 equations in sympy, find 2 points on the 2 equations where the 2 points are closest.
Thank you
CodePudding user response：
You can set up the equations to be solved like this:
In [94]: x1 = symbols('x1')
In [95]: x2 = symbols('x2')
In [96]: y1 = x1**2
In [97]: y2 = -(x2 - 3)**2
In [98]: distance_squared = (y1 - y2)**2 (x1 - x2)**2
In [99]: distance_squared
Out[99]:
2
2 ⎛ 2 2⎞
(x₁ - x₂) ⎝x₁ (x₂ - 3) ⎠
In [100]: eq1, eq2 = [distance_squared.diff(x) for x in [x1, x2]]
In [101]: eq1
Out[101]:
⎛ 2 2⎞
4⋅x₁⋅⎝x₁ (x₂ - 3) ⎠ 2⋅x₁ - 2⋅x₂
In [102]: eq2
Out[102]:
⎛ 2 2⎞
-2⋅x₁ 2⋅x₂ ⎝x₁ (x₂ - 3) ⎠⋅(4⋅x₂ - 12)
Now you have a pair of polynomial equations in x1 and x2. You just need to find the solutions of those equations. Note that there are limits on what you can do here depending on the degree of the polynomials due to the Abel-Ruffini theorem. A general approach that can give the solutions numerically is by using a Groebner basis:
In [106]: g = groebner([eq1, eq2], [x1, x2])
In [107]: g
Out[107]:
⎛⎡ 4 3 2 5 4 3
GroebnerBasis⎝⎣181⋅x₁ - 24⋅x₂ 212⋅x₂ - 624⋅x₂ 737⋅x₂ - 432, 8⋅x₂ - 96⋅x₂ 472⋅x₂ - 1206⋅
2 ⎤ ⎞
x₂ 1656⋅x₂ - 999⎦, x₁, x₂, domain=ℤ, order=lex⎠
In [108]: r = real_roots(g[-1])
In [109]: r
Out[109]:
⎡ ⎛ 3 2 ⎞⎤
⎣CRootOf⎝4⋅x - 36⋅x 110⋅x - 111, 0⎠⎦
In [112]: solve([g[0], x2 - r[0].n()], [x1, x2]) # numerical result
Out[112]: [(0.728082123067953, 2.27191787693205)]
In [111]: solve([g[0], x2 - r[0]], [x1, x2]) # exact result in terms of rootof
Out[111]:
⎡⎛ 3
⎢⎜ ⎛ 3 2 ⎞ ⎛ 3 2 ⎞
⎢⎜- 212⋅CRootOf⎝4⋅x - 36⋅x 110⋅x - 111, 0⎠ - 737⋅CRootOf⎝4⋅x - 36⋅x 110⋅x - 111, 0⎠ 432
⎢⎜─────────────────────────────────────────────────────────────────────────────────────────────────
⎣⎝ 181
4 2
⎛ 3 2 ⎞ ⎛ 3 2 ⎞
24⋅CRootOf⎝4⋅x - 36⋅x 110⋅x - 111, 0⎠ 624⋅CRootOf⎝4⋅x - 36⋅x 110⋅x - 111, 0⎠
──────────────────────────────────────────────────────────────────────────────────────────, CRootOf
⎞⎤
⎟⎥
⎛ 3 2 ⎞⎟⎥
⎝4⋅x - 36⋅x 110⋅x - 111, 0⎠⎟⎥
⎠⎦
In simple cases like this it is possible to solve the equations symbolically using solve but the answers can be very complicated:
In [114]: solve([eq1, eq2], [x1, x2])
Out[114]:
⎡⎛ __
⎢⎜ ╱ 8
⎢⎜ 3 _____________ 5/6 3 ___________ 3 ╱ ─
⎢⎜- 15066259⋅╲╱ 81 3⋅√753 - 549045⋅√251⋅3 ⋅╲╱ 27 √753 410109939 14945237⋅√753 ╲╱ 8
⎢⎜───────────────────────────────────────────────────────────────────────────────────────, - ──────
⎢⎜ 2/3 ⎛ 6 ___ 2/3⎞
⎢⎜ (27 √753) ⋅⎝1647135⋅√251⋅╲╱ 3 15066259⋅3 ⎠
⎢⎜
⎣⎝
___________ ⎞ ⎛
1 3⋅√753 ⎟ ⎜
─ ────── ⎟ ⎜ 3 _____________ 5/6 3 ___________
8 1 ⎟ ⎜1647135⋅√251⋅╲╱ 81 3⋅√753 15066259⋅3 ⋅╲╱ 27 √753
─────────── ─────────────────── 3⎟, ⎜──────────────────────────────────────────────────────────
3 _____________ ⎟ ⎜ 2/3 ⎛
╱ 81 3⋅√753 ⎟ ⎜ (27 √753) ⋅⎝1647135⋅√25
2⋅3 ╱ ── ────── ⎟ ⎜
╲╱ 8 8 ⎠ ⎝
5/6 3 ___________ 3 _____________
549045⋅√251⋅3 ⋅ⅈ⋅╲╱ 27 √753 15066259⋅ⅈ⋅╲╱ 81 3⋅√753 820219878⋅ⅈ 29890474⋅√753⋅ⅈ
────────────────────────────────────────────────────────────────────────────────────────────────, 3
2/3 6 ___ 2/3 6 ___ ⎞
1⋅3 45198777⋅╲╱ 3 - 15066259⋅3 ⋅ⅈ - 1647135⋅√251⋅╲╱ 3 ⋅ⅈ⎠
_____________ ⎞ ⎛
⎛ 1 √3⋅ⅈ⎞ ╱ 81 3⋅√753 ⎟ ⎜
⎜- ─ ────⎟⋅3 ╱ ── ────── ⎟ ⎜ 3 ____________
⎝ 2 2 ⎠ ╲╱ 8 8 1 ⎟ ⎜1647135⋅√251⋅╲╱ 81 3⋅√753
- ────────────────────────────── ────────────────────────────────⎟, ⎜───────────────────────────
3 _____________⎟ ⎜
⎛ 1 √3⋅ⅈ⎞ ╱ 81 3⋅√753 ⎟ ⎜
2⋅⎜- ─ ────⎟⋅3 ╱ ── ────── ⎟ ⎜
⎝ 2 2 ⎠ ╲╱ 8 8 ⎠ ⎝
_ 5/6 3 ___________ 3 _____________
15066259⋅3 ⋅╲╱ 27 √753 - 29890474⋅√753⋅ⅈ - 820219878⋅ⅈ - 15066259⋅ⅈ⋅╲╱ 81 3⋅√753 - 5490
───────────────────────────────────────────────────────────────────────────────────────────────────
2/3 ⎛ 2/3 6 ___ 6 ___ 2/3 ⎞
(27 √753) ⋅⎝1647135⋅√251⋅3 45198777⋅╲╱ 3 1647135⋅√251⋅╲╱ 3 ⋅ⅈ 15066259⋅3 ⋅ⅈ⎠
_____________
⎛ 1 √3⋅ⅈ⎞ ╱ 81 3⋅√753
5/6 3 ___________ ⎜- ─ - ────⎟⋅3 ╱ ── ──────
45⋅√251⋅3 ⋅ⅈ⋅╲╱ 27 √753 1 ⎝ 2 2 ⎠ ╲╱ 8 8
────────────────────────────, 3 ──────────────────────────────── - ──────────────────────────────
_____________ 3
⎛ 1 √3⋅ⅈ⎞ ╱ 81 3⋅√753
2⋅⎜- ─ - ────⎟⋅3 ╱ ── ──────
⎝ 2 2 ⎠ ╲╱ 8 8
⎞ ⎤
⎟ ⎥
⎟ ⎥
⎟ ⎛3 3⋅ⅈ 3 3⋅ⅈ⎞ ⎛3 3⋅ⅈ 3 3⋅ⅈ⎞⎥
⎟, ⎜─ - ───, ─ - ───⎟, ⎜─ ───, ─ ───⎟⎥
⎟ ⎝2 2 2 2 ⎠ ⎝2 2 2 2 ⎠⎥
⎟ ⎥
⎟ ⎥
⎠ ⎦