Determining Constants with Optimization Toolkit

I am attempting to determine two constants (E0 and sigma) such that a given equation:

D=F*(tau^2)*(t1/tau + exp(-t1/tau) -1) + (tau*(t2-t1))/lambda + tau*(sqrt(sigma*tau))*((-1/lambda)*(sqrt(tau/sigma)) - atanh((1/lambda)*sqrt(tau/sigma)) - sqrt(1 + ((tau/(lambda^2)*sigma)-1)*exp((-2/tau)*(T - t2))) + atanh(sqrt(1 + ((tau/(lambda^2)*sigma)-1)*exp((-2/tau)*(T-t2)))));

approximates a set of points as closely as possible. I was told this might be accomplished through the optimization toolkit, but I cannot see how. The difficulty is that E0 and sigma are required in this code to find the values of D, the value attempting to be matched to the data set:

x = [43.03,100.9,132,206,284.8,440.7,757.4,1577.5,3386,4345.4,5207.4,7377];
y = [400,800,1000,1500,2000,3000,5000,10000,20000,25000,30000,42195];

The y-values are the values of D. My original thought was to somehow iterate through both all values of T and E0 and sigma simultaneously, then stripping out the best fit at the end. But I am not sure how to accomplish this. This is what I have written, though I have only included initial values for E0 and sigma and have not attempted to iterate-through and sort for best values of these two constants.

sigma = 0 % for initial iteration only
E0 = 0 % for initial iteration only
F = 7.5168 
tau = 1.5011
%F and tau are known values calculated using nlinfit for 0-200 m distances

[Tc,Ts]=findTcTs(E0,sigma,F,tau);  %from Keller (3.3)

T = Ts
T0 = Tc
for k = 1:83001 %to determine t1 and t2 for all possible intervals of seconds from 20 seconds to 2 hours
    [t1,t2,lambda,D] = findt1t2lambdaD(E0,sigma,F,tau,T,T0)
    T0 = t2 %to update accurate initial guess for t2 for next iteration
    T = T + 0.1
end

%Equations for determining D (y-values) from t1, t2, lambda (with T being
%given on each iteration)
lfn=@(t)(sqrt(tau/sigma*(1-4*exp(-2*(T-t)/tau))));%Eq. (3.15)
t1fn=@(t)(-tau*log(1-1/(lfn(t)*F)));%Eq. (3.12) (solved in terms of t1)
Et2=@(t)(real(E0+sigma*t-tau^2/(2*lfn(t)^2)...
-F^2*tau*(-3*tau/2+t1fn(t)+2*tau*exp(-t1fn(t)/tau)-tau/2*exp(-2*t1fn(t)/tau))...
-tau/lfn(t)^2*(t-t1fn(t))));%Eq. (3.14). Note real part is taken to avoid complex values.

t2=fzero(Et2,t20);%solve for t2.
t1=t1fn(t2);
lambda=lfn(t2);
D=F*(tau^2)*(t1/tau + exp(-t1/tau) -1) + (tau*(t2-t1))/lambda + tau*(sqrt(sigma*tau))*((-1/lambda)*(sqrt(tau/sigma)) - atanh((1/lambda)*sqrt(tau/sigma)) - sqrt(1 + ((tau/(lambda^2)*sigma)-1)*exp((-2/tau)*(T - t2))) + atanh(sqrt(1 + ((tau/(lambda^2)*sigma)-1)*exp((-2/tau)*(T-t2)))));
%D from Keller (3.17)

modelfn = @(b,x)(%D eqn above...);
x = [43.03,100.9,132,206,284.8,440.7,757.4,1577.5,3386,4345.4,5207.4,7377];
y = [400,800,1000,1500,2000,3000,5000,10000,20000,25000,30000,42195];
beta0 = [] %initial guess dependent on number of variables
beta = nlinfit(x,y,modelfn,beta0)

Answer 1

[Alan Weiss]

Perhaps this example will help. Or, if you have an Optimization Toolbox™ license, see this Nonlinear Data-Fitting example.

Alan Weiss

MATLAB mathematical toolbox documentation