Open Educational Resources

Curve Fitting: Nonlinear Regression

Nonlinear Regression

In nonlinear regression, the model function y is a nonlinear function of x and of the parameters a_1, a_2,\cdots,a_m. Given a set of n data points (x_i,y_i) with 1\leq i \leq n, curve fitting starts by assuming a model y(x) with m parameters. The parameters can be obtained by minimizing the least squares:

    \[ S=\sum_{i=1}^n (y(x_i)-y_i)^2\]

In order to find the parameters of the model that would minimize S, m equations of the following form are solved:

    \[ \frac{\partial S}{\partial a_k}=\sum_{i=1}^n\left(2\left(y(x_i)-y_i\right)\frac{\partial y}{\partial a_k}\bigg|_{x=x_i}\right)=0\]

Since y(x) is nonlinear in the coefficients, the m equations formed are also nonlinear and can only be solved using a nonlinear equation solver method such as the Newton Raphson method described before. We are going to rely on the built-in NonlinearModelFit function in Mathematica that does the required calculations.

Example 1

The following data describes the yield strength \sigma_y and the plastic strain \varepsilon_p obtained from uniaxial experiments where the first variable is the plastic strain and the second variable is the corresponding yield strength: (0.0001,204.18),(0.0003,226.4),(0.0015,270.35),(0.0025,276.86),(0.004,296.86),(0.005,299.3),(0.015,334.65),(0.025,346.56),(0.04,371.81),(0.05,377.45),(0.1,398.01),(0.15,422.45),(0.2,434.42).
Use the NonlinearModelFit function in Mathematica to find the best fit for the data assuming the relationship follows the form:

    \[\sigma_y=K\varepsilon_p^n\]

Solution

While this model can be linearized as described in the previous section, the NonlinearModelFit will be used. The array of data is input as a list in the variable Data. The NonlinearModelFit function is then used and the output is stored in the variable: “model”. The form of the best fit function is evaluated using the built-in function Normal[model]. R^2 can be evaluated as well as shown in the code below. The values of K and n that describe the best fit are K=506.64 and n=0.0988. Therefore, the best-fit model is:

    \[\sigma_y=506.64\varepsilon_p^{0.0988}\]


The following is the Mathematica code used along with the output
View Mathematica Code
Data = {{0.0001, 204.18}, {0.0003, 226.4}, {0.0015, 270.35}, {0.0025, 276.86}, {0.004, 296.86}, {0.005, 299.3}, {0.015, 334.65}, {0.025,346.56}, {0.04, 371.81}, {0.05, 377.45}, {0.1, 398.01}, {0.15, 422.45}, {0.2, 434.42}};
model = NonlinearModelFit[Data, Kk*eps^(n), {Kk, n}, eps]
y = Normal[model]
Print["R^2=",R2 = model["RSquared"]]
Plot[y, {eps, 0, 0.21}, Epilog -> {PointSize[Large], Point[Data]}, PlotLegends -> {"Model"}, AxesLabel -> {"Eps_p", "Sigma_y"}, AxesOrigin -> {0, 0}]
View Python Code
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

Data = [[0.0001, 204.18], [0.0003, 226.4], [0.0015, 270.35], [0.0025, 276.86], [0.004, 296.86], [0.005, 299.3], [0.015, 334.65], [0.025,346.56], [0.04, 371.81], [0.05, 377.45], [0.1, 398.01], [0.15, 422.45], [0.2, 434.42]]
def f(eps, k, n): return k*eps**n
coeff, covariance = curve_fit(f, [point[0] for point in Data],
                                 [point[1] for point in Data])
print("coeff: ",coeff)
x_val = np.arange(0,0.21,0.001)
plt.title('%.5feps**(%.5f)' % tuple(coeff))
plt.plot(x_val, f(x_val, coeff[0], coeff[1]))
plt.scatter([point[0] for point in Data], [point[1] for point in Data], c='k')
plt.xlabel("Eps_p"); plt.ylabel("Sigma_y")
plt.grid(); plt.show()

# R squared
x = np.array([point[0] for point in Data])
y = np.array([point[1] for point in Data])
y_fit = f(x, coeff[0], coeff[1])
ss_res = np.sum((y - y_fit)**2)
ss_tot = np.sum((y - np.mean(y))**2)
r2 = 1 - (ss_res / ss_tot)
print("R Squared: ",r2)

nonlinearregres

The following links provide the MATLAB codes for implementing a nonlinear model fit for yield strength vs plastic strain. The function for the model is provided in File 2.

Example 2

Exponential and Logarithmic Models offer a wide range of model functions that can be used for various engineering applications. An example is the exponential model that describes an initial rapid increase in the dependent variable y which then levels off to become asymptotic to an upper limit a_1:

    \[y(x)=a_1(1-e^{-a_2x})\]

where a_1 and a_2 are the model parameters. For plastic materials, Armstrong and Frederick proposed this model in 1966 (republished here in 2007) to describe how a quantity called the “Centre of the Yield Surface” changes with the increase in the plastic strain. If \alpha is the centre of the yield surface and \varepsilon_p is the plastic strain, then, the model for \alpha has the form:

    \[\alpha = C\left(1-e^{-\gamma \varepsilon_p}\right)\]

Find the values of C and \gamma for the best fit to the experimental data ({\varepsilon_p}_i,\alpha_i) given by: (0.0001,2.067),(0.005,3.96),(0.0075,6.22),(0.008,8.07),(0.009,9.32),(0.01,12.02),(0.02,13.2),(0.03,15.22),(0.04,15.51),(0.05,15.44),(0.06,15.1),(0.075,15.76),(0.1,15.11),(0.17,15.53),(0.2,15.26)

Solution

Unlike the model in the previous question, this model cannot be linearized. The NonlinearModelFit built-in function in Mathematica will be used. The array of data is input as a list in the variable Data. The NonlinearModelFit function is then used and the output is stored in the variable: “model”. The form of the best fit function is evaluated using the built-in function Normal[model]. R^2 can be evaluated as well as shown in the code below. The values of C and \gamma that describe the best fit are C=15.5624 and \gamma=94.3836. Therefore, the best-fit model is:

    \[\alpha = 15.5624\left(1-e^{-94.3836 \varepsilon_p}\right)\]


The following is the Mathematica code used along with the output

View Mathematica Code
Data = {{0.0001, 2.067}, {0.005, 3.96}, {0.0075, 6.22}, {0.008,8.07}, {0.009, 9.32}, {0.01, 12.02}, {0.02, 13.2}, {0.03,15.22}, {0.04, 15.51}, {0.05, 15.44}, {0.06, 15.1}, {0.075,15.76}, {0.1, 15.11}, {0.17, 15.53}, {0.2, 15.26}};
model = NonlinearModelFit[Data, Cc (1 - E^(-gamma*eps)), {Cc, gamma}, eps]
y = Normal[model]
Print["R^2=", R2 = model["RSquared"]]
Plot[y, {eps, 0, 0.21}, Epilog -> {PointSize[Large], Point[Data]}, PlotLegends -> {"Model"}, AxesLabel -> {"Eps_p", "Alpha"}, AxesOrigin -> {0, 0}]
View Python Code
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

Data = [[0.0001, 2.067], [0.005, 3.96], [0.0075, 6.22], [0.008,8.07], [0.009, 9.32], [0.01, 12.02], [0.02, 13.2], [0.03,15.22], [0.04, 15.51], [0.05, 15.44], [0.06, 15.1], [0.075,15.76], [0.1, 15.11], [0.17, 15.53], [0.2, 15.26]]
def f(eps, c, gamma): return c*(1 - np.exp(-gamma*eps))
coeff, covariance = curve_fit(f, [point[0] for point in Data],
                                 [point[1] for point in Data])
print("coeff: ",coeff)
x_val = np.arange(0,0.21,0.001)
plt.title('%.5f(1 - e**(-%.5feps))' % tuple(coeff))
plt.plot(x_val, f(x_val, coeff[0], coeff[1]))
plt.scatter([point[0] for point in Data], [point[1] for point in Data], c='k')
plt.xlabel("Eps_p"); plt.ylabel("Alpha")
plt.grid(); plt.show()

# R squared
x = np.array([point[0] for point in Data])
y = np.array([point[1] for point in Data])
y_fit = f(x, coeff[0], coeff[1])
ss_res = np.sum((y - y_fit)**2)
ss_tot = np.sum((y - np.mean(y))**2)
r2 = 1 - (ss_res / ss_tot)
print("R Squared: ",r2)

graph

The following links provide the MATLAB codes for implementing a Nonlinear model fit for centre of the yield surface. The function for the model is provided in File 2.

Lecture Video

Leave a Reply

Your email address will not be published.