How can I fit an exponential curve?
66 views (last 30 days)
Show older comments
fengsen huang
on 15 Nov 2018
Answered: Arturo Gonzalez
on 8 Sep 2020
Below is an example of finding a fit with only one term of exponential term but I dont know how to find the fit of the curve when it has 2 degree of exponential term, i.e.[y = a*e^(bx) + c*e^(dx)]
example for y = a*e^(bx)
phi =[ones(size(xx)),xx];
aa=phi\log(yy);
yfit = exp(phi*aa);
plot(xx, yy, ’ro’, xx, yfit, ’k-’) ;
s=sprintf(’y=%8fexp(%8fx)’,exp(aa(1)),aa(2));
0 Comments
Accepted Answer
Star Strider
on 15 Nov 2018
Edited: Star Strider
on 15 Nov 2018
Try this:
filename1 = 'x2.mat';
D1 = load(filename1);
x = D1.x2;
filename2 = 'y2.mat';
D2 = load(filename2);
y = D2.y2;
fcn = @(b,x) b(1).*exp(b(2).*x) + b(3).*exp(b(4).*x);
[B,fval] = fminsearch(@(b) norm(y - fcn(b,x)), ones(4,1));
figure
plot(x, y, 'p')
hold on
plot(x, fcn(B,x), '-')
hold off
grid
The fitted parameters are:
B =
1.13024777245481
-2.75090020576997
-2.09865110252493
-5.48051415288241
EDIT — Added plot figure.
9 Comments
Star Strider
on 15 Nov 2018
Edited: Star Strider
on 15 Nov 2018
This line:
fcn = @(b,x) b(1).*exp(b(2).*x) + b(3).*exp(b(4).*x);
is the objective function, the expression that describes the function to fit to the data.
This line:
[B,fval] = fminsearch(@(b) norm(y - fcn(b,x)), ones(4,1))
calls the fminsearch function to fit the function to the data. The norm function compares the function output to the data and returns a single scalar value (the square root of the sum of squares of the difference between the function evaluation and the data here), that fminsearch uses. I refer you to the documentation on fminsearch (link) for details on how it works.
The value would be calculated as:
Rsq = 1 - sum((y - fcn(B,x)).^2) / sum((y - mean(y)).^2)
returning:
Rsq =
0.980214434988184
We are not comparing models, so this is the only statistic available. There are several ways to compare models, a subject much more involved than I will go into here. See any good text on nonlinear parameter estimation for details.
@fengsen huang —
If my Answer helped you solve your problem, please Accept it!
More Answers (2)
Arturo Gonzalez
on 8 Sep 2020
Per this answer, you can do it with the following matlab code
clear all;
clc;
% get data
dx = 0.001;
x = (dx:dx:1.5)';
y = -1 + 5*exp(0.5*x) + 4*exp(-3*x) + 2*exp(-2*x);
% calculate n integrals of y and n-1 powers of x
n = 3;
iy = zeros(length(x), n);
xp = zeros(length(x), n+1);
iy(:,1) = cumtrapz(x, y);
xp(:,1) = x;
for ii=2:1:n
iy(:, ii) = cumtrapz(x, iy(:, ii-1));
xp(:, ii) = xp(:, ii-1) .* x;
end
xp(:, n+1) = ones(size(x));
% get exponentials lambdas
Y = [iy, xp];
A = pinv(Y)*y;
Ahat = [A(1:n)'; [eye(n-1), zeros(n-1, 1)]];
lambdas = eig(Ahat);
lambdas
% get exponentials multipliers
X = [ones(size(x)), exp(lambdas'.*x)];
P = pinv(X)*y;
P
% show estimate
y_est = X*P;
figure();
plot(x, y); hold on;
plot(x, y_est, 'r--');
0 Comments
See Also
Categories
Find more on Interpolation in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!