Solving a linear-quadratic system of equations, both graphically and numerically, in numpy/matplotlib?

Question

I have a system of a linear equation and a quadratic equation that I can set up with numpy and scipy so I can get a graphical solution. Consider the example code:

#!/usr/bin/env python
# Python 2.7.1+

import numpy as np #
import matplotlib.pyplot as plt #

# d is a constant;
d=3
# h is variable; depends on x, which is also variable

# linear function:
# condition for h: d-2x=8h; returns h
def hcond(x):
  return (d-2*x)/8.0

# quadratic function:  
# condition for h: h^2+x^2=d*x ; returns h
def hquad(x):
  return np.sqrt(d*x-x**2)

# x indices data
xi = np.arange(0,3,0.01)

# function values in respect to x indices data
hc = hcond(xi)
hq = hquad(xi)

fig = plt.figure() 
sp = fig.add_subplot(111)

myplot = sp.plot(xi,hc)
myplot2 = sp.plot(xi,hq)

plt.show()

That code results with this plot:

It's clear that the two functions intersect, thus there is a solution.

How could I automatically solve what is the solution (the intersection point), while keeping most of the function definitions intact?

this can be done analytically too ? or is your real problem a bit harder than just polynomials of first n second degree? — usethedeathstar, Jun 16 '14 at 09:21
Thanks for the comment, @usethedeathstar - this was indeed my real problem; if by analytically you mean I could have solved the equations manually myself, you are correct; I just wanted an "auto-solution" as a "sanity check" device (to double-check my manual solving). Cheers! — sdaau, Jun 16 '14 at 10:31

score 4 · Answer 1 · edited May 23 '17 at 10:29

It turns out one can use scipy.optimize.fsolve to solve this, just need to be careful that the functions in the OP are defined in the y=f(x) format; while fsolve will need them in the f(x)-y=0 format. Here is the fixed code:

#!/usr/bin/env python
# Python 2.7.1+

import numpy as np #
import matplotlib.pyplot as plt #
import scipy
import scipy.optimize

# d is a constant;
d=3
# h is variable; depends on x, which is also variable

# linear function:
# condition for h: d-2x=8h; returns h
def hcond(x):
  return (d-2*x)/8.0

# quadratic function:
# condition for h: h^2+x^2=d*x ; returns h
def hquad(x):
  return np.sqrt(d*x-x**2)

# for optimize.fsolve;
# note, here the functions must be equal to 0;
# we defined h=(d-2x)/8 and h=sqrt(d*x-x^2);
# now we just rewrite in form (d-2x)/16-h=0 and sqrt(d*x-x^2)-h=0;
# thus, below x[0] is (guess for) x, and x[1] is (guess for) h!
def twofuncs(x):
  y = [ hcond(x[0])-x[1], hquad(x[0])-x[1] ]
  return y

# x indices data
xi = np.arange(0,3,0.01)

# function values in respect to x indices data
hc = hcond(xi)
hq = hquad(xi)

fig = plt.figure()
sp = fig.add_subplot(111)

myplot = sp.plot(xi,hc)
myplot2 = sp.plot(xi,hq)


# start from x=0 as guess for both functions
xsolv = scipy.optimize.fsolve(twofuncs, [0, 0])
print(xsolv)
print("xsolv: {0}\n".format(xsolv))

# plot solution with red marker 'o'
myplot3 = sp.plot(xsolv[0],xsolv[1],'ro')


plt.show()

exit

... which results with:

xsolv: [ 0.04478625  0.36380344]

... or, on the plot image:

Refs:

Solving a linear-quadratic system of equations, both graphically and numerically, in numpy/matplotlib?

1 Answers1