Environmental Engineering Reference
In-Depth Information
function par_est
% parameter estimation with derivatives
% for exponential fit for lambda
global tfit cfit c0
% specify fitting data
tfit = [0.25 1 2 4 8];
cfit = [0.7716 0.5791 0.4002 0.1860 0.1019];
c0 = 0.8; lambda0 = 0.5;
lambda = fzero(@myfun,lambda0);
normc = norm(cfit-c0*exp(-lambda*tfit));
display (['Best fit for lambda = ' num2str(lambda)]);
display (['Norm of residuals =' num2str(normc)]);
tmax = tfit(size(tfit,2));
t = [0:0.01*tmax:tmax];
figure; plot (tfit,cfit,'or',t,c0*exp(-lambda*t),'-');
legend ('given','modelled');
text(0.5*tmax,c0*0.7,['\lambda: ' num2str(lambda)]);
text(0.5*tmax,c0*0.8,['norm of residuals: ' num2str(normc)]);
function f = myfun(lambda);
global tfit cfit c0
c = c0*exp(-lambda*tfit); %solve linear decay eq. for c with c(0)=c0
f = (cfit-c)*(c.*tfit)'; % specify function f to vanish
The corresponding M-file
'parest.m'
is included in the accompanying
software.
The vectors of fitting data
tfit
and
cfit
as well as the initial concentration
c0
are
now specified in the input section of the M-file. Also the initial value for
(in the
M-file called
lambda0
) is specified at the beginning.
normc
, computed directly after
the calculation of
l
, represents the norm of the residuals. The next two lines initiate
output of
lambda
and
normc
in the MATLAB
l
command window. The last four
commands in the main module produce a plot showing given values and best fit
curve, the best fit value and the norm of residuals. The following figure is the plot
obtained for the microcystins example data introduced at the beginning of this
chapter (Fig.
10.5
).
The optimum l is 0.3329 and the norm of residuals is 0.0646. Obviously, the
function with the new value of
®
is a much better fit than the exponential fit of the
previous sub-chapter that was obtained by polynomial curve fitting.
Using the calculated
l
one may have the idea to improve the fit further by
changing
c
0
in formula (
10.2
). For this purpose, the
'par_est.m'
file has to be
changed slightly in order to optimize for
c
0
and not for
l
l
. Instead of condition
(
10.8
) one obtains:
@c
0
¼
2
X
cðt
fit
Þc
fit
@c
@c
0
ðt
fit
Þ¼
0
@e
(10.12)
The factor 2 can be omitted. Evaluation of the last factor delivers the condition:
X
cðt
fit
Þc
fit
exp
ð
l
t
fit
Þ¼
0
(10.13)
