scipy.optimize.minimize to find intersection of curves

Question

to find intersection between 2 polynomials I tried to minimize sum of least-squares like in this code - first plotted correct solution, then printed result of minimization.

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(0,1,100)    # ???????? [0,1]

def polyD(x):
    return  1.115355004199118 - 1.597163991790283* x**1 + 0.6539311181514963* x**2
def polyS(x):
    return -0.03070291735792586 + 0.1601011622660309* x**1 + 0.8530319920733438* x**2

root= 0.61      # x - suppose found somewhere earlier...
y_root= polyD(root) # y
print(root, y_root) # x=0.61 y=0.38441273827121725

plt.plot(root, y_root, 'yo',   x, polyD(x), 'r-', x, polyS(x), 'b-', ms=20, )
plt.show()

####################

from scipy.optimize import minimize

x0=  0
res= minimize(lambda t: sum((polyD(x) - polyS(x))**2), x0)

print(res)
print(res.fun)  # 36.59096676853359 ????????????

I can assume a little difference due to the algorithm, but here I see the difference of 10^2 order.

What is my mistake?? How can I correct it? If I should define some constraints to x=[0,1]? or how to get correct result about the intersection_point?

P.S. or, perhaps, it will be easier to solve with derivatives (minimizing to find extremum)? EDIT: partial_derivatives of constrained function(s) (Jacobian matrix) can be used to speed-up the code execution - here

Answer 1

Don't use minimize and don't use fsolve; use Polynomial and its built-in .roots:

poly_d = Polynomial((1.115355004199118, -1.597163991790283, 0.6539311181514963))
poly_s = Polynomial((-0.03070291735792586, 0.1601011622660309, 0.8530319920733438))

for x in (poly_d - poly_s).roots():
    print(f'd({x:.2f}) = {poly_d(x):.2f}')
    print(f's({x:.2f}) = {poly_s(x):.2f}')

d(-9.44) = 74.41
s(-9.44) = 74.41
d(0.61) = 0.38
s(0.61) = 0.38

Or since this is only order 2, you may also solve it analytically.

c, b, a = poly_d - poly_s
u = 0.5/a
d = np.sqrt(b*b - 4*a*c)
x0 = u*(d - b)
x1 = u*(-d - b)

Answer 2

The reason your code doesn't work is because of this line:

lambda t: sum((polyD(x) - polyS(x))**2)

Notice anything wrong? You used x instead of t. Swap those xs for ts and you get your desired results.

Alternatively, you could use scipy.optimize.root. Below I show both methods, though I only apply bounds to the minimize method.

import numpy as np
from scipy.optimize import minimize, root


def polyD(x):
    return 1.115355004199118 - 1.597163991790283*x**1 + 0.6539311181514963*x**2


def polyS(x):
    return -0.03070291735792586 + 0.1601011622660309*x**1 + 0.8530319920733438*x**2


def f(x):
    return polyD(x) - polyS(x)


def g(x):
    return np.sum(f(x)**2)


x0 = 0
min_res = minimize(g, x0, bounds=[(0, 1)])
root_res = root(f, x0)

print(min_res)
print(root_res)

Output:

  message: CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL
  success: True
   status: 0
      fun: 7.999974469252813e-13
        x: [ 6.100e-01]
      nit: 5
      jac: [-3.538e-06]
     nfev: 12
     njev: 6
 hess_inv: <1x1 LbfgsInvHessProduct with dtype=float64>

 message: The solution converged.
 success: True
  status: 1
     fun: [ 6.772e-15]
       x: [ 6.100e-01]
    nfev: 7
    fjac: [[-1.000e+00]]
       r: [ 2.000e+00]
     qtf: [ 6.131e-09]

Answer 3

another solution, probably: to use difference of functions as eq constraint:

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize

x = np.linspace(0,1,100)

def fd(x):
    return  1.115355004199118 - 1.597163991790283* x**1 + 0.6539311181514963* x**2
def fs(x):
    return -0.03070291735792586 + 0.1601011622660309* x**1 + 0.8530319920733438* x**2

plt.plot(  x, fs(x), 'r-', x, fd(x), 'b-')
plt.show()

######################

cons = ( { 'type' : 'eq', 'fun': lambda x: fs(x)-fd(x)} )

res = minimize(fs, [0., ],   constraints= cons,  tol=1.4e-100)
if res.success== True: print(res.x , ' ~ ', fd(res.x), '\n')
else: res.message

plt.plot(res.x, fd(res.x), 'yo', x, fd(x), 'r-', x, fs(x), 'b-', ms=20, )
plt.show()

really, code_design can vary, giving needed result, I suppose.

As I understand: we can use constraints for ys & bounds for xs. (?)