The code worked for one of the datasets but not the other. I receive errors (shown at the bottom of my provided code).
You can see that I have two commented out data sections by type: One for work density, one for strain. I am trying to plot the 6 datasets, present in each type, as independent curves against the x data.
I don't need all 12 plotted- I am just exchanging the strain data for the work density data as I need.
I am unable to succeed. I am sure the code is correct and should work but something is too calculation-intensive...?
The expected results: A plot of the fitting functions (logistic sigmoid) for 6 datasets, then a subplot included of the derivative of their theoretical curves. Returned values of R^2 squared values (I cannot use mean-squared error).
# This code produces a 99-100% fitness of fit to Normal Fibril strain v. temperature data, or, work
density v. temperature data, using a technical
# 4-parameter logistic sigmoid function, but with the centroid parameter pre-defined at 77.5 degrees
Celsius on the x-axis. This 'b' parameter is
# left undefined for the derivative function.
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from sklearn.metrics import r2_score
# 1. Normal Fibrils user dataset:
x_data = np.array([25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100, 105, 110, 115,
120, 125, 130])
y_data50 = np.array([0, 0, 0.250155, 1.00062, 1.751085, 2.50155, 3.252015, 4.752945, 6.253875,
7.00434, 8.50527, 10.756665, 13.00806, 16.260075, 18.51147, 19.261935, 19.51209, 19.762245, 19.762245,
20.0124, 20.0124, 20.0124])
y_data100 = np.array([0, 0, 0, 0, 1.98162, 3.467835, 5.94486, 7.431075, 9.9081, 13.87134, 16.348365,
18.82539, 22.78863, 28.73349, 34.67835, 40.62321, 42.60483, 43.59564, 44.58645, 45.081855, 45.081855,
45.081855])
y_data150 = np.array([0, 0, 0, 0, 2.221965, 4.44393, 7.40655, 10.36917, 13.33179, 17.035065, 22.21965,
25.18227, 31.10751, 37.03275, 44.4393, 51.84585, 62.21502, 68.14026, 71.10288, 74.0655, 75.54681,
75.54681])
y_data200 = np.array([0, 0, 0, -1.97181, -1.97181, 3.94362, 7.88724, 13.80267, 18.732195, 25.63353,
31.54896, 39.4362, 45.35163, 51.26706, 62.112015, 73.942875, 86.75964, 96.61869, 106.47774, 108.44955,
110.42136, 112.39317])
y_data250 = np.array([0, 0, 0, 2.46231, 4.92462, 9.84924, 14.77386, 20.929635, 27.08541, 34.47234,
41.85927, 54.17082, 64.02006, 73.8693, 86.18085, 103.41702, 118.19088, 130.50243, 137.88936,
140.35167, 142.81398, 145.27629])
y_data300 = np.array([0, 0, 1.476405, 2.95281, 4.429215, 7.382025, 13.287645, 20.66967, 28.051695,
38.38653, 51.674175, 67.91463, 82.67868, 100.39554, 121.06521, 137.305665, 150.59331, 157.975335,
165.35736, 168.31017, 171.26298, 171.26298])
# Normal Fibrils' Work Density Data:
# m(w) = 50m(f): y_data50 = np.array([0, 0, 0.250155, 1.00062, 1.751085, 2.50155, 3.252015, 4.752945,
6.253875, 7.00434, 8.50527, 10.756665, 13.00806, 16.260075, 18.51147, 19.261935, 19.51209, 19.762245,
19.762245, 20.0124, 20.0124, 20.0124])
# m(w) = 100m(f): y_data100 = np.array([0, 0, 0, 0, 1.98162, 3.467835, 5.94486, 7.431075, 9.9081,
13.87134, 16.348365, 18.82539, 22.78863, 28.73349, 34.67835, 40.62321, 42.60483, 43.59564, 44.58645,
45.081855, 45.081855, 45.081855])
# m(w) = 150m(f): y_data150 = np.array([0, 0, 0, 0, 2.221965, 4.44393, 7.40655, 10.36917, 13.33179,
17.035065, 22.21965, 25.18227, 31.10751, 37.03275, 44.4393, 51.84585, 62.21502, 68.14026, 71.10288,
74.0655, 75.54681, 75.54681])
# m(w) = 200m(f): y_data200 = np.array([0, 0, 0, -1.97181, -1.97181, 3.94362, 7.88724, 13.80267,
18.732195, 25.63353, 31.54896, 39.4362, 45.35163, 51.26706, 62.112015, 73.942875, 86.75964, 96.61869,
106.47774, 108.44955, 110.42136, 112.39317])
# m(w) = 250m(f): y_data250 = np.array([0, 0, 0, 2.46231, 4.92462, 9.84924, 14.77386, 20.929635,
27.08541, 34.47234, 41.85927, 54.17082, 64.02006, 73.8693, 86.18085, 103.41702, 118.19088, 130.50243,
137.88936, 140.35167, 142.81398, 145.27629])
# m(w) = 300m(f): y_data300 = np.array([0, 0, 1.476405, 2.95281, 4.429215, 7.382025, 13.287645,
20.66967, 28.051695, 38.38653, 51.674175, 67.91463, 82.67868, 100.39554, 121.06521, 137.305665,
150.59331, 157.975335, 165.35736, 168.31017, 171.26298, 171.26298])
# Normal Fibrils' Strain Data:
# x_data = np.array([25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100, 105, 110, 115,
120, 125, 130])
# m(w) = 50m(f): y_data50 = np.array([0, 0, 0.00617284, 0.024691358, 0.043209877, 0.061728395,
0.080246914, 0.117283951, 0.154320988, 0.172839506, 0.209876543, 0.265432099, 0.320987654,
0.401234568, 0.456790123, 0.475308642, 0.481481481, 0.487654321, 0.487654321, 0.49382716, 0.49382716,
0.49382716])
# m(w) = 100m(f): y_data100 = np.array([0, 0, 0, 0, 0, 0, 0.021978022, 0.038461538, 0.065934066,
0.082417582, 0.10989011, 0.153846154, 0.181318681, 0.208791209, 0.252747253, 0.318681319, 0.384615385,
0.450549451, 0.472527473, 0.483516484, 0.494505495, 0.5])
# m(w) = 150m(f): y_data150 = np.array([0, 0, 0, 0, 0.015151515, 0.03030303, 0.050505051, 0.070707071,
0.090909091, 0.116161616, 0.151515152, 0.171717172, 0.212121212, 0.252525253, 0.303030303,
0.353535354, 0.424242424, 0.464646465, 0.484848485, 0.505050505, 0.515151515, 0.515151515])
# m(w) = 200m(f): y_data200 = np.array([0, 0, 0, -0.009345794, -0.009345794, 0.018691589, 0.037383178,
0.065420561, 0.088785047, 0.121495327, 0.14953271, 0.186915888, 0.214953271, 0.242990654, 0.294392523,
0.35046729, 0.411214953, 0.457943925, 0.504672897, 0.514018692, 0.523364486, 0.53271028])
# m(w) = 250m(f): y_data250 = np.array([0, 0, 0, 0.009090909, 0.018181818, 0.036363636, 0.054545455,
0.077272727, 0.1, 0.127272727, 0.154545455, 0.2, 0.236363636, 0.272727273, 0.318181818, 0.381818182,
0.436363636, 0.481818182, 0.509090909, 0.518181818, 0.527272727, 0.536363636])
# m(w) = 300m(f): y_data300 = np.array([0, 0, 0.004504505, 0.009009009, 0.013513514, 0.022522523,
0.040540541, 0.063063063, 0.085585586, 0.117117117, 0.157657658, 0.207207207, 0.252252252,
0.306306306, 0.369369369, 0.418918919, 0.459459459, 0.481981982, 0.504504505, 0.513513514,
0.522522523, 0.522522523])
# 2. Four-parameter sigmoid fitting function:
def sigmoid(x, a, b, c, U):
return (U / (1 + np.exp(-c * (x - b))))**a
# 3. Sigmoid derivative function:
def sigmoid_derivative(x, a, b, c, U):
return (a * U * np.exp(-c * (x - b)) * (1 + np.exp(-c * (x - b)))**(-a - 1) * c) / (1 + np.exp(-c
* (x - b)))**(2*a)
# 4. Initial guess for the fitting parameters:
initial_guess = [1, 77.5, 0.1, 0.5]
# 5. Fit the sigmoid function to each dataset and calculate R-squared value:
popt50, _ = curve_fit(sigmoid, x_data, y_data50, p0=initial_guess)
y_pred50 = sigmoid(x_data, *popt50)
r_squared_50 = r2_score(y_data50, y_pred50)
popt100, _ = curve_fit(sigmoid, x_data, y_data100, p0=initial_guess)
y_pred100 = sigmoid(x_data, *popt100)
r_squared_100 = r2_score(y_data100, y_pred100)
popt150, _ = curve_fit(sigmoid, x_data, y_data150, p0=initial_guess)
y_pred150 = sigmoid(x_data, *popt150)
r_squared_150 = r2_score(y_data150, y_pred150)
popt200, _ = curve_fit(sigmoid, x_data, y_data200, p0=initial_guess)
y_pred200 = sigmoid(x_data, *popt200)
r_squared_200 = r2_score(y_data200, y_pred200)
popt250, _ = curve_fit(sigmoid, x_data, y_data250, p0=initial_guess)
y_pred250 = sigmoid(x_data, *popt250)
r_squared_250 = r2_score(y_data250, y_pred250)
popt300, _ = curve_fit(sigmoid, x_data, y_data300, p0=initial_guess)
y_pred300 = sigmoid(x_data, *popt300)
r_squared_300 = r2_score(y_data300, y_pred300)
# 6. Generate x values for the fitted curve:
x_fit = np.linspace(min(x_data), max(x_data), 1000)
# 7. Calculate the y values using the fitted parameters for each dataset:
y_fit50 = sigmoid(x_fit, *popt50)
y_fit100 = sigmoid(x_fit, *popt100)
y_fit150 = sigmoid(x_fit, *popt150)
y_fit200 = sigmoid(x_fit, *popt200)
y_fit250 = sigmoid(x_fit, *popt250)
y_fit300 = sigmoid(x_fit, *popt300)
# 8. Calculate the derivatives of each fitted curve:
dy_dx_fit50 = sigmoid_derivative(x_fit, *popt50)
dy_dx_fit100 = sigmoid_derivative(x_fit, *popt100)
dy_dx_fit150 = sigmoid_derivative(x_fit, *popt150)
dy_dx_fit200 = sigmoid_derivative(x_fit, *popt200)
dy_dx_fit250 = sigmoid_derivative(x_fit, *popt250)
dy_dx_fit300 = sigmoid_derivative(x_fit, *popt300)
# 9. Plotting:
plt.figure(figsize=(15, 5))
# 10. Plot sigmoid fits:
plt.subplot(1, 2, 1)
plt.plot(x_fit, y_fit300, '-', color = 'black', label = f'm(w) = 300m(f)|(R^2={r_squared_300:.2f})')
plt.plot(x_fit, y_fit250, '-', color = 'purple', label = f'm(w) = 250m(f)|(R^2={r_squared_250:.2f})')
plt.plot(x_fit, y_fit200, '-', color = 'red', label = f'm(w) = 200m(f)|(R^2={r_squared_200:.2f})')
plt.plot(x_fit, y_fit150, '-', color = 'orange', label = f'm(w) = 150m(f)|(R^2={r_squared_150:.2f})')
plt.plot(x_fit, y_fit100, '-', color = 'gold', label = f'm(w) = 100m(f)|(R^2={r_squared_100:.2f})')
plt.plot(x_fit, y_fit50, '-', color = 'green', label = f'm(w) = 50m(f)|(R^2={r_squared_50:.2f})')
plt.xlabel('Temperature (°C)')
plt.ylabel('Work Density WD')
plt.title('Sigmoidal Fit to Normal Fibril Data [WD v. T]')
#plt.ylabel('Strain ε, (Δl/$\mathregular{l_{0}}$)')
#plt.title('Sigmoidal Fit to Normal Fibril Data [ε v. T]')
plt.legend()
plt.grid(True)
# 11. Plot derivatives:
plt.subplot(1, 2, 2)
plt.plot(x_fit, dy_dx_fit300, '--', color = 'black', label = 'Rate for m(w) = 300m(f)')
plt.plot(x_fit, dy_dx_fit250, '--', color = 'purple', label = 'Rate for m(w) = 250m(f)')
plt.plot(x_fit, dy_dx_fit200, '--', color = 'red', label = 'Rate for m(w) = 200m(f)')
plt.plot(x_fit, dy_dx_fit150, '--', color = 'orange', label = 'Rate for m(w) = 150m(f)')
plt.plot(x_fit, dy_dx_fit100, '--', color = 'gold', label = 'Rate for m(w) = 100m(f)')
plt.plot(x_fit, dy_dx_fit50, '--', color = 'green', label='Rate for m(w) = 50m(f)')
plt.xlabel('Temperature (°C)')
plt.ylabel('ΔWD/ΔT, (Nm/kg°C)')
#plt.ylabel('Δε/ΔT, ((Δl/$\mathregular{l_{0}}$)/°C)')
plt.title('Contraction Rate')
plt.legend()
plt.grid(True)
plt.show()
The errors:
'Runtime Warning: divide by zero encountered in power'
'Runtime Warning: invalid value encountered in power'
'Optimize Warning: Covariance of the parameters could not be estimated.'
'line 52, in <module> popt100, _ ..... line 982, in curve_fit raise Runtime Error("Optimal parameters not found:)'
'Optimal parameters not found: Number of calls to function has reached maxfev = 1000.'
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from sklearn.metrics import r2_score
# 1. User datasets:
x_data = np.array([25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100, 105, 110, 115, 120, 125, 130])
y_dataWD50 = np.array([0, 0, 0.250155, 1.00062, 1.751085, 2.50155, 3.252015, 4.752945, 6.253875, 7.00434, 8.50527, 10.756665, 13.00806, 16.260075, 18.51147, 19.261935, 19.51209, 19.762245, 19.762245, 20.0124, 20.0124, 20.0124])
y_dataWD100 = np.array([0, 0, 0, 0, 1.98162, 3.467835, 5.94486, 7.431075, 9.9081, 13.87134, 16.348365, 18.82539, 22.78863, 28.73349, 34.67835, 40.62321, 42.60483, 43.59564, 44.58645, 45.081855, 45.081855, 45.081855])
y_dataWD150 = np.array([0, 0, 0, 0, 2.221965, 4.44393, 7.40655, 10.36917, 13.33179, 17.035065, 22.21965, 25.18227, 31.10751, 37.03275, 44.4393, 51.84585, 62.21502, 68.14026, 71.10288, 74.0655, 75.54681, 75.54681])
y_dataWD200 = np.array([0, 0, 0, -1.97181, -1.97181, 3.94362, 7.88724, 13.80267, 18.732195, 25.63353, 31.54896, 39.4362, 45.35163, 51.26706, 62.112015, 73.942875, 86.75964, 96.61869, 106.47774, 108.44955, 110.42136, 112.39317])
y_dataWD250 = np.array([0, 0, 0, 2.46231, 4.92462, 9.84924, 14.77386, 20.929635, 27.08541, 34.47234, 41.85927, 54.17082, 64.02006, 73.8693, 86.18085, 103.41702, 118.19088, 130.50243, 137.88936, 140.35167, 142.81398, 145.27629])
y_dataWD300 = np.array([0, 0, 1.476405, 2.95281, 4.429215, 7.382025, 13.287645, 20.66967, 28.051695, 38.38653, 51.674175, 67.91463, 82.67868, 100.39554, 121.06521, 137.305665, 150.59331, 157.975335, 165.35736, 168.31017, 171.26298, 171.26298])
# (A.) Normal Fibrils' Work Density Data:
# m(w) = 50m(f): y_data50 = np.array([0, 0, 0.250155, 1.00062, 1.751085, 2.50155, 3.252015, 4.752945, 6.253875, 7.00434, 8.50527, 10.756665, 13.00806, 16.260075, 18.51147, 19.261935, 19.51209, 19.762245, 19.762245, 20.0124, 20.0124, 20.0124])
# m(w) = 100m(f): y_data100 = np.array([0, 0, 0, 0, 1.98162, 3.467835, 5.94486, 7.431075, 9.9081, 13.87134, 16.348365, 18.82539, 22.78863, 28.73349, 34.67835, 40.62321, 42.60483, 43.59564, 44.58645, 45.081855, 45.081855, 45.081855])
# m(w) = 150m(f): y_data150 = np.array([0, 0, 0, 0, 2.221965, 4.44393, 7.40655, 10.36917, 13.33179, 17.035065, 22.21965, 25.18227, 31.10751, 37.03275, 44.4393, 51.84585, 62.21502, 68.14026, 71.10288, 74.0655, 75.54681, 75.54681])
# m(w) = 200m(f): y_data200 = np.array([0, 0, 0, -1.97181, -1.97181, 3.94362, 7.88724, 13.80267, 18.732195, 25.63353, 31.54896, 39.4362, 45.35163, 51.26706, 62.112015, 73.942875, 86.75964, 96.61869, 106.47774, 108.44955, 110.42136, 112.39317])
# m(w) = 250m(f): y_data250 = np.array([0, 0, 0, 2.46231, 4.92462, 9.84924, 14.77386, 20.929635, 27.08541, 34.47234, 41.85927, 54.17082, 64.02006, 73.8693, 86.18085, 103.41702, 118.19088, 130.50243, 137.88936, 140.35167, 142.81398, 145.27629])
# m(w) = 300m(f): y_data300 = np.array([0, 0, 1.476405, 2.95281, 4.429215, 7.382025, 13.287645, 20.66967, 28.051695, 38.38653, 51.674175, 67.91463, 82.67868, 100.39554, 121.06521, 137.305665, 150.59331, 157.975335, 165.35736, 168.31017, 171.26298, 171.26298])
# (B.) Normal Fibrils' Strain Data:
# x_data = np.array([25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100, 105, 110, 115, 120, 125, 130])
# m(w) = 50m(f): y_data50 = np.array([0, 0, 0.00617284, 0.024691358, 0.043209877, 0.061728395, 0.080246914, 0.117283951, 0.154320988, 0.172839506, 0.209876543, 0.265432099, 0.320987654, 0.401234568, 0.456790123, 0.475308642, 0.481481481, 0.487654321, 0.487654321, 0.49382716, 0.49382716, 0.49382716])
# m(w) = 100m(f): y_data100 = np.array([0, 0, 0, 0, 0, 0, 0.021978022, 0.038461538, 0.065934066, 0.082417582, 0.10989011, 0.153846154, 0.181318681, 0.208791209, 0.252747253, 0.318681319, 0.384615385, 0.450549451, 0.472527473, 0.483516484, 0.494505495, 0.5])
# m(w) = 150m(f): y_data150 = np.array([0, 0, 0, 0, 0.015151515, 0.03030303, 0.050505051, 0.070707071, 0.090909091, 0.116161616, 0.151515152, 0.171717172, 0.212121212, 0.252525253, 0.303030303, 0.353535354, 0.424242424, 0.464646465, 0.484848485, 0.505050505, 0.515151515, 0.515151515])
# m(w) = 200m(f): y_data200 = np.array([0, 0, 0, -0.009345794, -0.009345794, 0.018691589, 0.037383178, 0.065420561, 0.088785047, 0.121495327, 0.14953271, 0.186915888, 0.214953271, 0.242990654, 0.294392523, 0.35046729, 0.411214953, 0.457943925, 0.504672897, 0.514018692, 0.523364486, 0.53271028])
# m(w) = 250m(f): y_data250 = np.array([0, 0, 0, 0.009090909, 0.018181818, 0.036363636, 0.054545455, 0.077272727, 0.1, 0.127272727, 0.154545455, 0.2, 0.236363636, 0.272727273, 0.318181818, 0.381818182, 0.436363636, 0.481818182, 0.509090909, 0.518181818, 0.527272727, 0.536363636])
# m(w) = 300m(f): y_data300 = np.array([0, 0, 0.004504505, 0.009009009, 0.013513514, 0.022522523, 0.040540541, 0.063063063, 0.085585586, 0.117117117, 0.157657658, 0.207207207, 0.252252252, 0.306306306, 0.369369369, 0.418918919, 0.459459459, 0.481981982, 0.504504505, 0.513513514, 0.522522523, 0.522522523])
# 2. Sigmoid function defined:
def sigmoid(x, U, k, x0, a):
return U / (1 + np.exp(-k * (x - x0)))**a
# 3. Initial guesses for calculating fitting parameters:
p0_WD50 = [max(y_dataWD50), -0.1, 77.5, 1.0]
p0_WD100 = [max(y_dataWD100), -0.1, 77.5, 1.0]
p0_WD150 = [max(y_dataWD150), -0.1, 77.5, 1.0]
p0_WD200 = [max(y_dataWD200), -0.1, 77.5, 1.0]
p0_WD250 = [max(y_dataWD250), -0.1, 77.5, 1.0]
p0_WD300 = [max(y_dataWD300), -0.1, 77.5, 1.0]
# 3b. Derivative of 'sigmoid' defined:
def sigmoid_derivative(x, U, k, x0, a):
exp_term = np.exp(-k * (x - x0))
sigmoid_term = (1 + exp_term) ** (-k)
return a * (U / (1 + exp_term)) ** (a - 1) * U * k * exp_term * sigmoid_term
# 4. Sigmoid fitting operation:
x_fit = np.linspace(25, 130, 100)
# 4a. Estimated parameters for fitting & R^2 fittness of fit values:
# (i.) m(w) = 50m(f)
poptWD50, pcov = curve_fit(sigmoid, x_data, y_dataWD50, p0_WD50, maxfev = 10000)
y_estWD50 = sigmoid(x_data, *poptWD50)
r_squaredWD50 = r2_score(y_dataWD50, y_estWD50)
print(f"R^2 Value: {r_squaredWD50:.3f}")
# (ii.) m(w) = 100m(f)
poptWD100, pcov = curve_fit(sigmoid, x_data, y_dataWD100, p0_WD100, maxfev = 10000)
y_estWD100 = sigmoid(x_data, *poptWD100)
r_squaredWD100 = r2_score(y_dataWD100, y_estWD100)
print(f"R^2 Value: {r_squaredWD100:.3f}")
# (iii.) m(w) = 150m(f)
poptWD150, pcov = curve_fit(sigmoid, x_data, y_dataWD150, p0_WD150, maxfev = 10000)
y_estWD150 = sigmoid(x_data, *poptWD150)
r_squaredWD150 = r2_score(y_dataWD150, y_estWD150)
print(f"R^2 Value: {r_squaredWD100:.3f}")
# (iv.) m(w) = 200m(f)
poptWD200, pcov = curve_fit(sigmoid, x_data, y_dataWD200, p0_WD200, maxfev = 10000)
y_estWD200 = sigmoid(x_data, *poptWD200)
r_squaredWD200 = r2_score(y_dataWD200, y_estWD200)
print(f"R^2 Value: {r_squaredWD200:.3f}")
# (v.) m(w) = 250m(f)
poptWD250, pcov = curve_fit(sigmoid, x_data, y_dataWD250, p0_WD250, maxfev = 10000)
y_estWD250 = sigmoid(x_data, *poptWD250)
r_squaredWD250 = r2_score(y_dataWD250, y_estWD250)
print(f"R^2 Value: {r_squaredWD250:.3f}")
# (vi.) m(w) = 300m(f)
poptWD300, pcov = curve_fit(sigmoid, x_data, y_dataWD300, p0_WD300, maxfev = 10000)
y_estWD300 = sigmoid(x_data, *poptWD300)
r_squaredWD300 = r2_score(y_dataWD300, y_estWD300)
print(f"R^2 Value: {r_squaredWD300:.3f}")
# 4b. Calculation of Sigmoid Fit y-values & derivative y-values using the estimated parameters per dataset:
# (i.) [m(w) = 50m(f)]:
y_fitWD50 = sigmoid(x_fit, *poptWD50)
dy_dx_fitWD50 = sigmoid_derivative(x_fit, *poptWD50)
# (ii.) [m(w) = 100m(f)]:
y_fitWD100 = sigmoid(x_fit, *poptWD100)
dy_dx_fitWD100 = sigmoid_derivative(x_fit, *poptWD100)
# (iii.) [m(w) = 150m(f)]:
y_fitWD150 = sigmoid(x_fit, *poptWD150)
dy_dx_fitWD150 = sigmoid_derivative(x_fit, *poptWD150)
# (iv.) [m(w) = 200m(f)]:
y_fitWD200 = sigmoid(x_fit, *poptWD200)
dy_dx_fitWD200 = sigmoid_derivative(x_fit, *poptWD200)
# (v.) [m(w) = 250m(f)]:
y_fitWD250 = sigmoid(x_fit, *poptWD250)
dy_dx_fitWD250 = sigmoid_derivative(x_fit, *poptWD250)
# (vi.) [m(w) = 300m(f)]:
y_fitWD300 = sigmoid(x_fit, *poptWD300)
dy_dx_fitWD300 = sigmoid_derivative(x_fit, *poptWD300)
# 6. Print the logistic sigmoid fitting parameters:
print("Estimated Parameters:", poptWD50)
print("Estimated Parameters:", poptWD100)
print("Estimated Parameters:", poptWD150)
print("Estimated Parameters:", poptWD200)
print("Estimated Parameters:", poptWD250)
print("Estimated Parameters:", poptWD300)
# 7. Plot data, logistic sigmoid function fit & the derivative:
plt.figure(figsize = (9, 10))
# 7a. Datas' & Sigmoid Fits' subplot:
plt.subplot(2, 1, 1)
plt.plot(x_fit, y_fitWD50, '-', color = 'green', label=f'm(w) = 50m(f) (R^2 Value = {r_squaredWD50:.3f})')
plt.plot(x_fit, y_fitWD100, '-', color = 'yellow', label=f'm(w) = 100m(f) (R^2 Value = {r_squaredWD100:.3f})')
plt.plot(x_fit, y_fitWD150, '-', color = 'gold', label=f'm(w) = 150m(f) (R^2 Value = {r_squaredWD150:.3f})')
plt.plot(x_fit, y_fitWD200, '-', color = 'darkorange', label=f'm(w) = 200m(f) (R^2 Value = {r_squaredWD200:.3f})')
plt.plot(x_fit, y_fitWD250, '-', color = 'red', label=f'm(w) = 250m(f) (R^2 Value = {r_squaredWD250:.3f})')
plt.plot(x_fit, y_fitWD300, '-', color = 'purple', label=f'm(w) = 300m(f) (R^2 Value = {r_squaredWD300:.3f})')
plt.xlabel('Temperature (°C)')
plt.ylabel('Work Density ' + '$\it{WD}$ ' + '(Nm/kg)')
plt.title('$\mathregular{[l_{i}}$ = $\mathregular{L_{0}]}$ '+ 'Logistic Sigmoidal Fit to Work Density Data')
plt.legend()
plt.grid(True)
# 7b. Derivatives of the sigmoids' subplot:
plt.subplot(2, 1, 2)
plt.plot(x_fit, dy_dx_fitWD50, '--', color = 'green', label='Rate for m(w) = 50m(f)')
plt.plot(x_fit, dy_dx_fitWD100, '--', color = 'yellow', label='Rate for m(w) = 100m(f)')
plt.plot(x_fit, dy_dx_fitWD150, '--', color = 'gold', label='Rate for m(w) = 150m(f)')
plt.plot(x_fit, dy_dx_fitWD200, '--', color = 'darkorange', label='Rate for m(w) = 200m(f)')
plt.plot(x_fit, dy_dx_fitWD250, '--', color = 'red', label='Rate for m(w) = 250m(f)')
plt.plot(x_fit, dy_dx_fitWD300, '--', color = 'purple', label='Rate for m(w) = 300m(f)')
plt.xlabel('Temperature (°C)')
plt.ylabel('Rate of Change, ' + '$\it{ΔWD/ΔT}$ ' + '(Nm/kg°C)')
plt.legend()
plt.title('Work Density Rates (Derivative of Characteristic Sigmoids)')
plt.show()