Parameter Estimation - Environmental Modeling - page 198

The corresponding M-file 'paresta.m' is included in the accompanying
software.
The result is c 0 ¼ 0.816 with an improved residual norm of 0.061. The reader

may check the result as an exercise.
It is also possible to combine the last two M-files to estimate both c 0 and

. Then,
two function modules have to be used within the M-file. In one function ( myfun )

l

is
optimized. Within that function, the second one ( myfun2 ) is called, in which c 0 is
optimized.

l

function par_estb
% parameter estimation with derivatives
% for exponential fit with c0 amd lambda as parameters
global tfit cfit c0
c0 = 0.816; lambda = 0.333; % start values

% specify fitting data
tfit = [0.25 1 2 4 8];
cfit = [0.7716 0.5791 0.4002 0.1860 0.1019];

lambda = fzero(@myfun,lambda);
normc = norm(cfit - c0*exp(-lambda*tfit))
display (['Best fit for lambda= ' num2str(lambda)]);
display (['Best fit for c0= ' num2str(c0)]);
display (['Norm of residuals= ' num2str(normc)]);
tmax = tfit(size(tfit,2));
t = [0:0.01*tmax:tmax];
plot (tfit,cfit,'or',t,c0*exp(-lambda*t),'-');
legend ('given','modelled');
text(0.5*tmax,c0*0.8,['\lambda: ' num2str(lambda)]);
text(0.5*tmax,c0*0.7,['c_0: ' num2str(c0)]);
text(0.5*tmax,c0*0.6,['norm of residuals: ' num2str(normc)]);

function f = myfun(lambda);
global tfit cfit c0

options = optimset;
c0 = fzero(@myfun2,c0,options,lambda);
display (['Best fit for c0 = ' num2str(c0)]);

c = c0*exp(-lambda*tfit); % solve linear decay eq. for c with c0
clambda = -c.*tfit; % equation for dc/dlambda
f = (c-cfit)*clambda'; % specify function f to vanish

function f = myfun2(c0,lambda);
global tfit cfit

c = c0*exp(-lambda*tfit); %solve linear decay eq. for c with c(0)=c0
cc0 = exp(-lambda*tfit); % equation for dc/dc0
f = (c-cfit)*cc0'; % specify function f to vanish

The corresponding M-file 'parestb.m' is included in the accompanying
software.
The formal parameter lambda needs to be added in the call of the second
function. In order to do that, options must be transferred, too. The standard options
set is obtained by using the options ¼ optimset command.