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()

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()

• python
• interpolation
• physics
• intersection
• curve

Share Improve this question Follow edited Mar 10 at 9:58 Stef 15.6k22 gold badges2222 silver badges3838 bronze badges asked Mar 9 at 13:55 Clem EntwistleClem Entwistle 2111 silver badge11 bronze badge 4

• maybe take some small data so you could calculate it on paper - and compare your calculation with every step of calculations in your code. – furas Commented Mar 9 at 14:30
• Hi. There's a lot of variables in your code and I'm a bit too lazy to understand each of their meanings. Could you please tell me whether I understood your problem correctly: you have a curve, represented by a 2-column numpy array of (x,y) values, and a straight line, represented by the two coefficients m and c of equation y=mx+c. You want to find the coordinates of intersection of the line and the curve (assuming there is only one intersection point) – Stef Commented Mar 10 at 9:46
• If your curve data is presented in two 1-dimensional arrays X and Y, and you have the two coefficients m and c of the line, you can write Y2 = Y - (m*X + c) and search the zeroes of Y2 using this method – Stef Commented Mar 10 at 9:54
• @Stef Thank you for replying, yep there are a lot and I haven't understood them either if I'm honest., I got it from professorkazarinoff.github.io/Engineering-Materials-Programming/…, which explains the exact process I am attempting to complete. My curve data is two 1-dimensional arrays, I'll give that method you linked a go thank you! – Clem Entwistle Commented Mar 10 at 12:34

Add a comment  | 

1 Answer 1

Sorted by: Reset to default 1

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)