Appendix A
Nonlinear Curve Fitting
A sample consists of a layer of aluminum and a layer of a composite coating is tested in a vacuum chamber by measuring its temperature as a function of time. The behavior of the sample temperature has a nonlinear dependence on the emissivity of the sample and the heat transfer coefficient h between the coating and the vacuum air.
Figure A-1. A sample enclosed within a testing chamber.
The unknown parameters h and may be obtained by fitting the model equation to experimental data as shown in Figure A.2 where the curve represents the model equation and the circles represent the data.
Figure A-2. Transient temperature of a typical sample.
To illustrate how this is done, first consider a portion of the graph in Figure A-2 that is re-plotted in Figure A-3. The relationship between the temperature Ti obtained from the model equation and the experimental value Ti,exp can be expressed generally as
Ti,exp = Ti (t; , h) + ei(A-2)
where ei is a random error that can be negative or positive. Ti is a function of the independent variable ti and the parameters h and . The random error is also called the residual, which is the difference between the calculated and measured values.
Figure A-3. Relationship between the model equation and the data
Nonlinear regression is based on determining the values of the parameters that minimize the sum of the squares of the residuals called an objective function Fobj.
Fobj = = =(A-3)
Where N is the number of data points or measured temperatures in this case. The temperature from equation A-1 can be expanded in a Taylor series around h and and curtailed after the first derivative.
Ti,j+1 = Ti,j + + h(A-4)
Where j is the guess and j+1 is the prediction, = j+1j, and h = hj+1hj. We have linearized the original model with respect to the parameters h and . Equation (A-4) can be substituted into Eq. (A-2) to yield
Ti,exp Ti,j = + h + ei(A-5a)
or in matrix form
{D} = [Zj]{A} + {E}(A-5b)
where [Zj] is the matrix of partial derivatives of the function (called the Jacobian matrix) evaluated at the guess j, the vector {D}contains the differences between the measure temperature and the calculated temperature at the guess j, the vector {A}contains the changes in the parameter values, and the vector {E} contains the residuals. It should be noted that as the final values of the parameters are obtained after the iterations vector {D} is the same as vector {E}.
[Zj] = , {D} = , {A} = , {E} =
We minimize the objective function
Fobj = = (A-3)
by taking its derivative with respect to each of the parameters and setting the resulting equation to zero.
= 2 = 0(A-6a)
= 2 = 0(A-6b)
This algorithm is Gauss-Newton method for minimizing the sum of the squares of the residuals between data and nonlinear functions. Equations (A-6a) and (A-6b) can be combined in a matrix form
[Zj]T{E} = 0(A-7)
where [Zj]T is the transpose of [Zj]. Let consider N = 3 so we can see the combination from (A-6a) and (A-6b) to (A-7).
= 0
Substitute {E} = {D} [Zj]{A} from Eq. (A-5b) into (A-7)
[Zj]T{{D} [Zj]{A}} = 0
or
[Zj]T[Zj]{A} = {[Zj]T{D}}(A-8)
The Jacobian matrix [Zj] may be evaluated numerically for the model equation (A-1).
(A-9a)
(A-9b)
Typically, can be chosen to be 0.01 and h can be chosen to be 0.01 W/m2K. Thus, the Gauss-Newton method consists of solving Eq. (A-8) for {A}, which can be employed to compute improved values for the parameters h and .
j+1 = j + (from {A})
hj+1 = hj + h (from {A})
This procedure is repeated until the solution converges that is until and h fall below an acceptable criterion. The Gauss-Newton method is a common algorithm that can be found in many numerical method texts. However, this description follows the notations and development by Chapra and Canale.
Example A-1
Fit the function T(t; , h) = (1 eht) to the data.
t / 0.25 / 0.75 / 1.25 / 1.75 / 2.25T / 0.28 / 0.75 / 0.68 / 0.74 / 0.79
Use initial guesses of h = 1 and = 1 for the parameters.
Solution
The partial derivatives of the function with respect to the parameters h and are
= 1 eht and = teht
[Zj] = =
The matrix multiplied by its transpose results in
[Zj]T[Zj] =
[Zj]T[Zj] =
The vector {D} consists of the differences between the measurements T and the model predictions T(t; , h) = (1 eht)
{D} = =
The vector {D} is pre-multiplied by [Zj]T to give
[Zj]T{D} = =
The vector {A} can be calculated by using MATLAB statement dA=ZjTZj\ZjTD({A} = {[Zj]T[Zj] \[Zj]T{D}})
{A} =
The next guesses for the parameters and h are
= 1 0.2715 = 0.7285
h = 1 + 0.5019 = 1.5019
Table A-1 lists the MATLAB program with the results of two iterations.
Table A-1 ______
% Gauss-Newton method
%
t=[0.25 0.75 1.25 1.75 2.25]';
T=[0.28 0.57 0.68 0.74 0.79]';
e=1;h=1;
Tmodel='e*(1-exp(-h*t))';
dTde='1-exp(-h*t)';dTdh='e*t.*exp(-h*t)';
for i=1:2
Zj=[eval(dTde) eval(dTdh)];
ZjTZj=Zj'*Zj
D=T-eval(Tmodel)
ZjTD=Zj'*D
dA=ZjTZj\ZjTD
e=e+dA(1);
h=h+dA(2);
fprintf('Iteration #%g: e = %8.4f, h = %8.4f\n',i,e,h)
end
> Gauss
ZjTZj =
2.3194 0.9489
0.9489 0.4404
D =
0.0588
0.0424
-0.0335
-0.0862
-0.1046
ZjTD =
-0.1534
-0.0366
dA =
-0.2715
0.5019
Iteration #1: e = 0.7285, h = 1.5019
ZjTZj =
3.0660 0.4162
0.4162 0.0780
D =
0.0519
0.0777
0.0629
0.0641
0.0863
ZjTD =
0.2648
0.0397
dA =
0.0625
0.1758
Iteration #2: e = 0.7910, h = 1.6777
The Matlab function fminsearch can also be used to fit the data to an expression with more than one parameter, T(t; , h) = (1 eht). Table A-2 lists the function required by fminsearch. Table A-3 lists the program that calls fminsearch to find the two parameters = 1 and h = 1. The program also plots the fitted results with the experimental data shown in Figure A-4.
______Table A-2 Matlab program to define the objective function ______
function y=nlin(p)
t=[.25 .75 1.25 1.75 2.25];
T=[.28 .75 0.68 0.74 0.79];
e=p(1);h=p(2);
Tc=e*(1-exp(-h*t));
y=sum((T-Tc).^2);
______Table A-3 Matlab program to find andh ______
clf
t=[.25 .75 1.25 1.75 2.25];
T=[.28 .75 0.68 0.74 0.79];
p=fminsearch('nlin',[1 1])
tp=.25:.1:2.25;
e=p(1);h=p(2);
Tc=e*(1-exp(-h*tp));
plot(tp,Tc,t,T,'o')
grid on
xlabel('t');ylabel('T')
legend('Fitted','Data')
> nlinear
p =
0.7754 2.3752
Figure A-4 Nonlinear Regression of T(t; , h) = (1 eht).
1