-2

I have data with one independent variable x and two dependent variables y1 and y2 as shown below:

 x      y1      y2
-1.5    16.25   1.02
-1.25   17      1.03
-1      15      1.03
-0.75   9       1.09
-0.5    5.9     1.15
-0.25   5.2     1.17
 0      4.77    1.19
+0.25   3.14    1.35
+0.5    2.5     1.54
+0.75   2.21    1.69
+1      1.91    1.96
+1.25   1.64    2.27
+1.5    1.52    2.56
+1.75   1.37    3.06
+2      1.24    4.12
+2.25   1.2     4.44
+2.5    1.18    4.95
+2.75   1.12    6.49
+3      1.07    10

So, here the value of x where y1 = y2 is somewhere around +1. How do I read the data and calculate this in python?

sfactor
  • 12,592
  • 32
  • 102
  • 152

1 Answers1

2

The naive solution goes like this:

txt = """-1.5    16.25   1.02
-1.25   17      1.03
-1      15      1.03
-0.75   9       1.09
-0.5    5.9     1.15
-0.25   5.2     1.17
 0      4.77    1.19
+0.25   3.14    1.35
+0.5    2.5     1.54
+0.75   2.21    1.69
+1      1.91    1.96
+1.25   1.64    2.27
+1.5    1.52    2.56
+1.75   1.37    3.06
+2      1.24    4.12
+2.25   1.2     4.44
+2.5    1.18    4.95
+2.75   1.12    6.49
+3      1.07    10"""

import numpy as np
# StringIO behaves like a file object, use it to simulate reading from a file
from StringIO import StringIO   

x,y1,y2=np.transpose(np.loadtxt(StringIO(txt)))
p1 = np.poly1d(np.polyfit(x, y1, 1))
p2 = np.poly1d(np.polyfit(x, y2, 1))
print 'equations: ',p1,p2

#y1 and y2 have to be equal for some x, that you solve for :
#  a x+ b = c x + d  --> (a-c) x= d- b

a,b=list(p1)
c,d=list(p2)
x=(d-b)/(a-c)
print 'solution x= ',x

output:

equations:   
-3.222 x + 7.323  
1.409 x + 1.686
solution x=  1.21717324767

But then you plot the 'lines':

import matplotlib.pyplot as p
%matplotlib inline
p.plot(x,y1,'.-')
p.plot(x,y2,'.-')

enter image description here

And you realize you can't use a linear assumption but for a few segments.

x,y1,y2=np.transpose(np.loadtxt(StringIO(txt)))
x,y1,y2=x[8:13],y1[8:13],y2[8:13]
p1 = np.poly1d(np.polyfit(x, y1, 1))
p2 = np.poly1d(np.polyfit(x, y2, 1))
print 'equations: ',p1,p2
a,b=list(p1)
c,d=list(p2)
x0=(d-b)/(a-c)
print 'solution x= ',x0
p.plot(x,y1,'.-')
p.plot(x,y2,'.-')

Output:

equations:   
-1.012 x + 2.968  
1.048 x + 0.956
solution x=  0.976699029126

enter image description here

Even now one could improve by leaving two more points out (looking very linear, but that can be coincidental for a few points).

x,y1,y2=np.transpose(np.loadtxt(StringIO(txt)))
x1,x2=x[8:12],x[9:13]
y1,y2=y1[8:12],y2[9:13]
p1 = np.poly1d(np.polyfit(x1, y1, 1))
p2 = np.poly1d(np.polyfit(x2, y2, 1))
print 'equations: ',p1,p2
a,b=list(p1)
c,d=list(p2)
x0=(d-b)/(a-c)
print 'solution x= ',x0

import matplotlib.pyplot as p
%matplotlib inline
p.plot(x1,y1,'.-')
p.plot(x2,y2,'.-')

Output:

equations:   
-1.152 x + 3.073  
1.168 x + 0.806
solution x=  0.977155172414

enter image description here

Possibly better would be to use more points and apply a 2nd order interpolation np.poly1d(np.polyfit(x,y1,2)) and then solve the equality for two 2nd order polynomials, which I leave as an exercise (quadratic equation) for the reader.

roadrunner66
  • 7,772
  • 4
  • 32
  • 38