I am trying to solve the equation below using fminsearch
, but I think the objective function is wrong.
How should I write the objective function or modify any other part of the code? This is basically a fitting problem, where the optimization procedure should fit the equation to the given data.
% Consider the following data:
Data = ...
[0.0000 5.8955
0.1000 3.5639
0.2000 2.5173
0.3000 1.9790
0.4000 1.8990
0.5000 1.3938
0.6000 1.1359
0.7000 1.0096
0.8000 1.0343
0.9000 0.8435
1.0000 0.6856
1.1000 0.6100
1.2000 0.5392
1.3000 0.3946
1.4000 0.3903
1.5000 0.5474
1.6000 0.3459
1.7000 0.1370
1.8000 0.2211
1.9000 0.1704
2.0000 0.2636];
% Let's plot these data points.
t = Data(:,1);
y = Data(:,2);
plot(t,y,'ro')
title('Data points')
hold on
% fit the function: y = c(1)*exp(-lam(1)*t) + c(2)*exp(-lam(2)*t)
%
% define the parameters in terms of one variable x:
% x(1) = c(1)
% x(2) = lam(1)
% x(3) = c(2)
% x(4) = lam(2)
%
% Then define the curve as a function of the parameters x and the data t:
F = @(x,t)(x(1)*exp(-x(2)*t) + x(3)*exp(-x(4)*t));
% We arbitrarily set our initial point x0 as follows: c(1) = 1,
% lam(1) = 1, c(2) = 1, lam(2) = 0:
x0 = [1 1 1 0];
% We run the solver and plot the resulting fit
options = optimset('TolFun',1e-5,'TolX',1e-5,'MaxFunEvals',10,'MaxIter',4000,'Display','iter');
[x,fval,exitflag,output] = fminsearch(F,x0,options)
plot(t,F(x,t))
hold off
You're right, your objective function doesn't make sense. You could do a least squares fitting. Then the objective function should be defined as:
F = @(x,t) (x(1)*exp(-x(2)*t) + x(3)*exp(-x(4)*t));
Obj = @(x) (norm(F(x, Data(:,1))-Data(:,2)));
Then
x0 = [1 1 1 0];
options = optimset('TolFun',1e-5, 'TolX', 1e-5, 'MaxFunEvals',1000, 'MaxIter', 4000,'Display','iter');
[x,fval,exitflag,output] = fminsearch(Obj,x0,options);
tp = 0:0.01:2;
plot(Data(:,1), Data(:,2), 'ro');
title('Data points')
hold on
plot(tp,F(x,tp));
hold off
gives me:
EDIT:
Assuming you already know a parameter and you want to use function handles, you could do it like this
p = ... % Your calculation to get the parameter. In your case x(3) from the previous F
F = @(x, p, t) (x(1)*exp(-x(2)*t) + p*exp(-x(3)*t));
helper = @(x, p) (norm(F(x, p, Data(:,1))-Data(:,2)));
Obj = @(x) (helper(x, p));
x0 = [1 1 0]; % Note that there's now one variable/parameter less
options = optimset('TolFun',1e-5, 'TolX', 1e-5, 'MaxFunEvals',1000, 'MaxIter', 4000,'Display','iter');
[x,fval,exitflag,output] = fminsearch(Obj,x0,options);
tp = 0:0.01:2;
plot(Data(:,1), Data(:,2), 'ro');
title('Data points')
hold on
plot(tp,F(x,p,tp)); % Note that you need to pass p to F
hold off