Environmental Engineering Reference
In-Depth Information
The second condition follows from the chain rule
@c=@D ¼ @c=@x
ð
Þ @x=@D
ð
Þ
and the second boundary condition given in (
10.31
). It follows:
D
0
¼
0 and
D
1
¼Q=D
2
.
The following M-file demonstrates the procedure for an example data set (see
Rom
2005
):
function
par_est2
% parameter estimation with derivatives
% Idea from FEMLAB - there R instead of Q
% see COMSOL News, Nr. 1, 2005, page 15
global
xfit cfit Q
% specify fitting data
xfit = [0.05:0.1:0.95];
cfit = [0.9256859756097451 0.7884908536585051 0.6665396341462926
0.559832317073104 0.4683689024389414 0.39214939024380824
0.33117378048770196 0.28544207317062964 0.25495426829258294
0.23971036585356142];
Q = -2;
D = fzero(@myfun,2);
display (['Best fit for D = ' num2str(D)]);
x = [0:0.01:1];
plot (xfit,cfit,'o',x,-(Q/D/2)*x.*x + (Q/D)*x + 1,'-');
legend ('given','modelled');
xlabel ('x'); ylabel ('c');
function
f = myfun(D);
global
xfit cfit Q
% solve diffusion equation for c with c(0)=1 and dc/dx(1)=0
c = -(Q/D/2)*xfit.*xfit + (Q/D)*xfit + 1;
% solve Poisson equation for dc/dD (cD) with boundary conditions
cD = (Q/D/D/2)*xfit.*xfit - (Q/D/D)*xfit;
% specify function de to vanish
f = 2*(c-cfit)*cD';
The corresponding M-file
'parest2.m'
is included in the accompanying
software.
The main program follows the same line given by the
'par_est'
M-files presented
in the previous subchapters. The
fzero
function is called with an initial guess
D ¼
2. The function
myfun
consists of three commands. The first evaluates
c
as
formulated in (
10.30
);
@c=@D
following
(
10.32
); the third specifies function f according to formula (
10.27
).
For
Q ¼
the second evaluates the derivative
2 the procedure delivers the correct result of
D ¼
1.312 ! The best fit
is visualized in Fig.
10.8
:
Hitherto, the analytical solution of the problem has been utilized again. The first
two commands in the function
myfun
can be replaced by calls of differential
equation solvers. In the first command, the (
10.25
) needs to be solved with regard
to the boundary conditions (
10.31
); in the second command, it is the differential
(
10.29
) with regard to conditions (
10.33
). The function f then becomes: