Clem Entwistle
Clem Entwistle

Reputation: 21

Intersection of two data sets, one linear and one not

I need to find the intersection of a straight line and a stress-strain curve to determine the offset yield strength.

Curve

I've been trying to use the following code, which successfully prints the curve and the line, but places the intersection point in the wrong place.

strain1 and stress1 are both rather large arrays extracted from excel and m and c are the corresponding values for the appropriate straight line for strain1 and stress1. Pls help this struggling physics student :)

strain_offset = strain1 + 0.002
stress_offset = m[0]*strain1 + c[0]
yield_strength_index = np.argwhere(np.diff(np.sign(stress1 - stress_offset)))[0][0]

def intersection_point(Ax1, Ay1, Ax2, Ay2, Bx1, By1, Bx2, By2):
    d = (By2-By1)*(Ax2-Ax1)-(Bx2-Bx1)*(Ay2-Ay1)
    if d:
        uA = ((Bx2-Bx1)*(Ay1-By1)-(By2-By1)*(Ax1-Bx1))/d
        uB = ((Ax2-Ax1)*(Ay1-By1)-(Ay2-Ay1)*(Ax1-Bx1))/d
    else:
        return None

    x_intersection = Ax1 + uA * (Ax2 - Ax1)
    y_intersection = Ay1 + uA * (Ay2 - Ay1)
    
    return x_intersection, y_intersection

    
first = yield_strength_index
second = first + 1
# A points from the stress strain curve
Ax1 = strain1[first]
Ay1 = stress1[first]
Ax2 = strain1[second]
Ay2 = stress1[second]
# B points from the offset line
Bx1 = strain_offset[first]
By1 = stress_offset[first]
Bx2 = strain_offset[second]
By2 = stress_offset[second]

# run our function that finds the intersection point
x_intersection, y_intersection = intersection_point(Ax1,Ay1,Ax2,Ay2,Bx1,By1,Bx2,By2)

print(x_intersection,y_intersection)

fig,ax = plt.subplots()


ax.plot(strain1,stress1)
ax.plot(strain_offset,stress_offset)
ax.plot(x_intersection,y_intersection,'go')

ax.set_ylim([0,50000000])
ax.set_xlim([0,0.035])

plt.show()

Upvotes: 1

Views: 55

Answers (1)

9769953
9769953

Reputation: 12261

Using what Stef already described in a comment, you can subtract the two curves and find the point where it crosses the x-axis, y = 0. Or, slightly differently, calculate the (vertical) distance between the two curves and find the minimum of that.

So you would get a function for the vertical distance that can look somewhat like the following:

dist = np.abs(y - func(x, m, c))

with func your linear function, and x and y your data.

and you can then do

i = np.argmin(dist)
xp, yp = x[i], y[i]

with xp and yp the coordinates of the intersection, to a good approximation.

The "good approximation" part is because

  • the actual point may be between two data points. If you want to be more precise, you would have to perform interpolation to your data (or even fitting it, if you have a model)
  • we've minimized for only the vertical distance, y. Not the actual distance, which involves both x and y.

Taken that all together, here's a possible example, using a logistic function to create some idealized data:

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt

def logistic(x, mu, s):
    y = 1 / (1 + np.exp(-(x - mu)/s)))
    return y

def func(x, m, c):  # simple line
    return m * x + c

def minfunc(x, y, m, c):  # function to minimize
    return np.abs(y - func(x, m, c))

# create some sample data
# Pick soimport numpy as np
import scipy as sp
import matplotlib.pyplot as plt

def logistic(x, mu, s):
    y = 1 / (1 + np.exp(-(x - mu)/s))
    return y

def func(x, m, c):  # simple line
    return m * x + c

def minfunc(x, y, m, c):  # function to minimize
    return np.abs(y - func(x, m, c))

# create some sample data
# Pick some values that produces roughly data 
# on the same scale as in the question
x0, x1 = 0.0, 0.035
m, c = 3e9, -5e7
x0, x1 = 0.0, 0.035
mu, scale = x1 / 2, 0.0002 / x1
xdata = np.linspace(x0, x1, 150)
ydata = 5e7 * logistic(xdata, mu, scale) - 2e6

# Find the "closest approach"
diff = minfunc(xdata, ydata, m, c)
i = np.argmin(diff)
xp, yp = xdata[i], ydata[i]

# Plot everything
plt.plot(xdata, ydata)
plt.plot(xdata, func(xdata, m, c))
plt.plot([xp], [yp], 'o', ms=10)
plt.ylim(0, 5e7)

enter image description here

Upvotes: 1

Related Questions