Reputation: 21
I need to find the intersection of a straight line and a stress-strain curve to determine the offset yield strength.
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
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
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)
Upvotes: 1