pythonmathequation-solvingquadratic

Solve a 2 equations systems of quadratic equations with 2 unknows

I would like to solve this system of equation for X and Y : enter image description here

Do you know if their exist a mathematical solution to this system ? And how to solve it ? Or a numerical solution ?

If yes, do you think their is a way to solve it with Python ?

I would like to apply this solution to a gas sensing problem. My problem is that I measure 2 photoacoustic signal with 2 two different sensors. Each signal is a function of the concentration of CO2 and relative humidity. The equation for each signal is : enter image description here

C_gaz is a CO2 gas concentration in ppm. For the moment I'm measuring CO2 in a lab room. So it should be a positive real number between 0 and 5000 (usually around 400). X_RH is the realtive humidity in the room. It's a percentage. It's a real number between 0 and 100 (Not 0 and 1). It's very unlikely to be above 80 and under 20.

If I substract the signal to the first coefficient a_0, I end up with the system of equation I meantionned in the beginning.

The coefficients for the first sensor are :

and for the second sensor :

The PA signal I'm measuring is a real positive number. For the first sensor, It usually is between 1e-03 and 20e-3 V. You can consider it to be 2.4e-03 V. For the second sensor it usually is between 1e-03 and 20e-03 V. You can consider it to be around 13e-03 V .

I'm monitoring this signal for a couple hours. So I have a numpy array of a couple thousand elements containing the measured signal every 2 seconds approximately.

The data is here : https://www.dropbox.com/scl/fo/vzyk3r38shp3z1sccu6z5/AFeLzb-Rndw6deqmBrjtcBQ?rlkey=9xg0q357r12zk6io9h8gox6xs&st=ubah66bk&dl=0

Solution

  • Well, this is a numerical solution and it's not vectorised, but at least you can get some solution.

    First the maths. Use each equation to write Y as a function of X from the quadratic formula. (Note that I'm using b rather than a' for the second equation - it minimises typos.)

    enter image description here

    Equate the two expressions, cross-multiply to remove the denominators and rearrange with square roots on the left:

    enter image description here

    where we have polynomials for later:

    enter image description here

    Square this once:

    enter image description here

    Rearrange again:

    enter image description here

    Square again to remove the square root:

    enter image description here

    Looking at the degree of each term, this IS a polynomial equation of degree 4 and can be dealt with by numpy’s polynomial package. All square roots have been squared away.

    Once you have X then you can get Y from either of the original quadratics.

    Note that the process of squaring removes sign-dependence, so may introduce some spurious roots. Just check your (X,Y) pairs against the original equations.

    Python Code

    There are two tests - intersection of two circles where I know precisely where the intersection points are and then your latest data. Note that some of the roots are complex - you will have to extract the pure real ones yourself (e.g. eliminating those with non-zero imaginary part).

    import cmath
import numpy as np
from numpy.polynomial import Polynomial

def evaluate( a, a1, a2, a12, a11, a22, x, y ):
    return a + a1*x + a2*y + a12*x*y  + a11*x**2 + a22*y**2

def biquadratic( a, a1, a2, a12, a11, a22, b, b1, b2, b12, b11, b22 ):
    Da = Polynomial( [ a2**2 - 4*a22*a, 2*a12*a2 - 4*a1*a22, a12**2-4*a11*a22 ] )
    Db = Polynomial( [ b2**2 - 4*b22*b, 2*b12*b2 - 4*b1*b22, b12**2-4*b11*b22 ] )
    R = Polynomial( [ b22*a2 - a22*b2, b22*a12-a22*b12 ] )
    P = 4 * a22 ** 2 * b22 ** 2 * Da * Db - ( R**2 - b22**2 * Da - a22**2 * Db ) ** 2

    xy = []
    SMALL = 1.0e-10
    for x in P.roots():
        y1 = ( -( a12 * x + a2 ) + cmath.sqrt( ( a12 * x + a2 ) ** 2 - 4 * a22 * ( a11 * x ** 2 + a1 * x + a ) ) ) / ( 2 * a22 )
        y2 = -( a12 * x + a2 ) / a22 - y1
        if abs( evaluate(a,a1,a2,a12,a11,a22,x,y1 ) ) < SMALL and abs( evaluate(b,b1,b2,b12,b11,b22,x,y1) ) < SMALL:
            xy.append( [ x, y1 ] )
        if abs( evaluate(a,a1,a2,a12,a11,a22,x,y2 ) ) < SMALL and abs( evaluate(b,b1,b2,b12,b11,b22,x,y2) ) < SMALL:
            xy.append( [ x, y2 ] )
    return xy


print( "TEST 1: intersection of two circles" )
a, a1, a2, a12, a11, a22 = ( -5,  4, 0, 0, 1, 1 )     # circle, radius 3, centre (-2,0)
b, b1, b2, b12, b11, b22 = ( -5, -4, 0, 0, 1, 1 )     # circle, radius 3, centre ( 2,0)
solution = biquadratic( a, a1, a2, a12, a11, a22, b, b1, b2, b12, b11, b22 )
print( "Intersection should be at (0,+-sqrt(5))" )
for x, y in solution: print( x, y )

print()

print( "TEST 2: given data" )
a, a1, a2, a12, a11, a22 = ( 2.210e-04, 6.071e-07, -1.942e-05, 9.131e-08,  1.258e-10, 2.269e-07 )
b, b1, b2, b12, b11, b22 = (-2.658e-04, 4.486e-06,  1.824e-05, 4.013e-07,  1.767e-11,-2.623e-07 )
solution = biquadratic( a, a1, a2, a12, a11, a22, b, b1, b2, b12, b11, b22 )
print( "Solution for x, y:" )
for x, y in solution: print( x, y )

    Output:

    TEST 1: intersection of two circles
Intersection should be at (0,+-sqrt(5))
0.0 (2.23606797749979+0j)
0.0 (-2.23606797749979+0j)
0.0 (2.23606797749979+0j)
0.0 (-2.23606797749979+0j)

TEST 2: given data
Solution for x, y:
(5.675136187595641+0j) (14.34454156969241+0j)
(7.478653532183469+0j) (67.94975301818636-0j)
(1660.6834584401354-1338.7179286020366j) (-10.795970969115467+0.4116308115063788j)
(1660.6834584401354+1338.7179286020366j) (-10.795970969115467-0.4116308115063788j)

    You can also plot your data to check it:

    enter image description here