
Curve fitting with polyfit
9 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Hi,
I am completely lost... Not so much with the concept, but still lost. How do I use the coefficients that polyfit is producing? How do I use them to find A and eA?
0 Kommentare
Antworten (2)
Image Analyst
am 26 Okt. 2018
I know that taking the log and then doing an ordinary least squares polynomial fit is the simple way, and the way we are taught, but it may not the the most accurate way. I'll show you how to do it with fitnlm() instead of polyfit. I also compute the value at T=813. It's just to teach you the way many/most people in MATLAB would do this - it's not your homework and not something you can turn in:
% Uses fitnlm() to fit a non-linear model (an exponential growth curve) through noisy data.
% Requires the Statistics and Machine Learning Toolbox, which is where fitnlm() is contained.
% Initialization steps.
clc; % Clear the command window.
close all; % Close all figures (except those of imtool.)
clear; % Erase all existing variables. Or clearvars if you want.
workspace; % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 20;
% Create the X coordinates from 0 to 20 every 0.5 units.
X = [773, 786, 797, 810, 810, 820, 834]
% Define function that the X values obey.
Y = [1.63, 2.95, 4.19, 8.13, 8.19, 14.9, 22.2]
% Plot the noisy initial data.
plot(X, Y, 'b*', 'LineWidth', 2, 'MarkerSize', 15);
grid on;
% Convert X and Y into a table, which is the form fitnlm() likes the input data to be in.
tbl = table(X', Y');
% Define the model as Y = a + exp(-b*x)
% Note how this "x" of modelfun is related to big X and big Y.
% x((:, 1) is actually X and x(:, 2) is actually Y - the first and second columns of the table.
R = 8314;
modelfun = @(b,x) b(1) * exp(b(2)./(R * x(:, 1)));
beta0 = [2e15, -2e8]; % Guess values to start with. Just make your best guess.
% Now the next line is where the actual model computation is done.
mdl = fitnlm(tbl, modelfun, beta0);
% Now the model creation is done and the coefficients have been determined.
% YAY!!!!
% Extract the coefficient values from the the model object.
% The actual coefficients are in the "Estimate" column of the "Coefficients" table that's part of the mode.
coefficients = mdl.Coefficients{:, 'Estimate'}
% Create smoothed/regressed data using the model:
% Make 1000 points between the min x and max x.
xFit = linspace(min(X), max(X), 1000);
yFitted = coefficients(1) * exp(coefficients(2) ./ (R*xFit));
% Now we're done and we can plot the smooth model as a red line going through the noisy blue markers.
hold on;
plot(xFit, yFitted, 'r-', 'LineWidth', 2);
grid on;
title('Exponential Regression with fitnlm()', 'FontSize', fontSize);
xlabel('T(K) [degrees]', 'FontSize', fontSize);
ylabel('k, Reaction Rate', 'FontSize', fontSize);
legendHandle = legend('Noisy Y', 'Fitted Y', 'Location', 'north');
legendHandle.FontSize = 30;
formulaString = sprintf(' k = %g * exp(%g / (R * T))', coefficients(1), coefficients(2))
xl = xlim;
yl = ylim;
text(xl(1), mean(yl), formulaString, 'FontSize', 25, 'FontWeight', 'bold');
% Find the value at X = 813.
y813 = coefficients(1) * exp(coefficients(2) ./ (R*813))
plot(813, y813, 'm*', 'LineWidth', 2, 'MarkerSize', 25);
line([813, 813], [yl(1), y813], 'Color', 'm', 'LineWidth', 2);
% Set up figure properties:
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
% Get rid of tool bar and pulldown menus that are along top of figure.
% set(gcf, 'Toolbar', 'none', 'Menu', 'none');
% Give a name to the title bar.
set(gcf, 'Name', 'Demo by ImageAnalyst', 'NumberTitle', 'Off')

If you're interested, you might compare the sum of the residuals with this method and compare it to the homework method.
1 Kommentar
Star Strider
am 26 Okt. 2018
It depends what ‘A’ and ‘eA’ are, and the polynomial degree you want to fit.
Assuming a linear fit:
p = polyfit(x,y,1)
m = p(1);
b = p(2);
in y=m*x+b.
6 Kommentare
Siehe auch
Kategorien
Mehr zu Linear and Nonlinear Regression finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!