How to do a nonlinear fit using least squares

I have a set of data points giving me the values for the second virial coefficient, for various values of T, of the virial expansion which is an equation that corrects the ideal gas law for empirical deviations:

I'm trying to do a least squares fit to determine how well the van der Waals equation predicts using MATLAB.

The equation you derive for B using van der Waals' equation ends up being:

where a and b are unknown constants I need to determine and R is just the ideal gas constant.

Now my data points for true values are here:

B = [-160 -35 -4.2 9 16.9 21.3]*10^-6 ;

T = [100 200 300 400 500 600]; (kelvin)

How can I write some code which will use the least squares method to generate estimates for a and b using the given data points?

My problem is that this is essentially a equation which I'm not sure how to represent in MATLAB.

I can only find options for quadratic, cubic, polynomial fit etc.

Answer 1

[Star Strider]

This is actually a linear problem, so a linear approximation will estimate the parameters correctly:

R = 8.314462;                                                     % J K^−1 mol
B = [-160 -35 -4.2 9 16.9 21.3];
T = [100 200 300 400 500 600];

ab = [ones(size(B(:))) - 1./(R*T(:))] \ B(:);                     % Linear Approximation Parameter Estimation
Bfit = [ones(size(B(:))) - 1./(R*T(:))] * ab;                     % Linear Fit

figure
plot(T, B, 'pg')
hold on
plot(T, Bfit,'-r')
hold off
grid
xlabel('T')
ylabel('B')
legend('Data','Linear Least-Squares Fit', 'Location','E')

The parameters are:

a = 64.2001320479734
b = 182307.574287957

EDIT —

Added plot figure:

[Andrew]

Yup, that's exactly right. Thanks!

[Star Strider]

My pleasure!

[Andrew]

Yeah that's what my graph looks like. I forgot to add that all B data points should be multiplied by a factor of 10^-6 but it was easy to figure out that the order of magnitude difference on your values was due to my mistake.

How were you able to determine the parameters a and b?

[Star Strider]

All you need to do is scale ‘B’, and the parameters will scale accordingly.

The parameters are the result of configuring the design matrix corectly, then using the mldivide,\ operator (or function) to do the linear fit. (If you want confidence intervals on the parameters as well, you can use essentially the same code with the regress function or its friends.)

Essentially the parameters are the x solution of where A is the design matrix:

[ones(size(B(:))) - 1./(R*T(:))]

and b here, is appropriately, ‘B’. Both ‘T’ and ‘B’ must be column vectors for this to work correctly. (Capital letters in general use denote matrices, and lower-case letters denote vectors.)

The solution x corresponds to the constant (or intercept) term (the first column of the design matrix) and the slope term, corresponding to the second column of the design matrix. Here, x is the column vector ‘ab’, so the first row is ‘a’, corresponding to the first column of the design matrix, and the second row is ’b’, corresponding to the second column of the desigh matrix. The calculation of ‘ab’ is essentially:

without the instabilities associates with calculating . (The mldivide function does not explicitly calculate the inverse, so the operation is more efficient and more mathematically stable.) Other than that, it is essentially straightforward linear algebra. (I also could have used polyfit, however this is more efficient.)

So

ab(1) = a
ab(2) = b

Explore the documentation on mldivide to understand in detail how it works. (The documentation explains it much better than I ever could.)

[Andrew]

gosh, MATLAB has endless capabilities. Thank you for the detailed explanation. Extremely informative. Thank you!

[Star Strider]

My pleasure!

Note that your equation is linear because it is linear in the parameters. That is because the partial derivatives of ‘B’ with respect to either parameter are not functions of that parameter or other parameters. So a linear approximation is appropriate, and provides the correct result.

[Andrew]

got it. Thanks for the clarification!

[Star Strider]

My pleasure!

Answer 2

[Jackson Burns]

Attached is code generated by cftool using a custom equation. Hope this helps!

[Andrew]

That's a very nifty tool! Thanks! It definitely generates the right curve but I'm not sure how to determine what equation defines the curve that best fits the data points. I need the equation to determine the constants a and b. I'm now wondering if the best way is to make a substitution to linearize the graph and then find the slope and y-intercept

[Jackson Burns]

If you use the function like this:

[result,~] = createFit(T,B)

result contains the coefficients of the function a and b.

Glad it helped!

[Andrew]

Pardon my ignorance but I'm having trouble implimenting the command:

[result,~] = createFit(T,B)

its telling me that "result" is unused. Is there something I need to link to that spot? Thanks again

[Jackson Burns]

If you're making this call in a function, it will tell you that result is unused if you don't later on reference result or return it from the function.

[Andrew]

got it. Thanks!