User:Sb~enwikibooks/Playground1

Matlab implementation
This Matlab implementation section is intended to be read while having a Matlab window opened. The code is supposed to be copied step by step into the Matlab environment by copy and paste.

Measurement data
As measurement data, the height (in meter) of the upper interface versus time (in seconds) is taken. The batch settling process has the homogenious initial concentration u0=0.123. The container has the height L=0.405, the measurements are observed up to the final time T=34800. global T H T=[0 120 240 360 480 600 720 840 960 1080 1200 1320 1440 1560 1680 1800 1920 2040 2160 2280 2400 2640 2880 3120 3360 3600 3840 4080 4320 4560 4800 5400 6000 6600 7200 7800 8400 9000 9600 10200 10800 11400 12000 13200 14400 15600 16800 20400 27600 34800]; H=[0.4050 0.3855 0.3700 0.3540 0.3390 0.3235 0.3075 0.2920 0.2760 0.2600 0.2440 0.2280 0.2140 0.2050 0.2000 0.1970 0.1945 0.1920 0.1900 0.1885 0.1870 0.1840 0.1815 0.1790 0.1770 0.1750 0.1735 0.1720 0.1705 0.1690 0.1680 0.1650 0.1620 0.1600 0.1580 0.1560 0.1545 0.1530 0.1510 0.1495 0.1485 0.1470 0.1460 0.1440 0.1425 0.1410 0.1400 0.1375 0.1350 0.1348]; plot(T,H,'x'); axis([0 34800 0 0.4050])

Mfiles
Subroutines are encapsulated by Mfiles (Matlab files with ending .m).

Direct problem
The closure functions are prescribed by function f = f(u) global uinf umax C f=uinf.*u.*(1-u).^C.*(u>0).*(uuc).*(u=0).*(u<=1)*length(Au))+(u<1)); b=dt/dx*([F(1:J-1);0]-[0;F(1:J-1)])-dt/dx^2*([Ax(1:J-1);0]-[0;Ax(2:J)]-[Ax(1:J-1);0]+[0;Ax(1:J-1)]); ax=a(u); epsilon=1e-5; F1eps=feo(g,u(2:J)+epsilon,u(1:J-1)); F2eps=feo(g,u(2:J),u(1:J-1)+epsilon); Jminus=(F1eps-F)/epsilon; Jplus=(F2eps-F)/epsilon; M1=-dt/dx*[Jminus(1:J-1)] - dt/dx^2*[ax(1:J-1)]; M2=1+dt/dx*([Jminus(1:J-1);0] - [0;Jplus(1:J-1)]) + 2*dt/dx^2*ax; M3=dt/dx*[Jplus(1:J-1)]-dt/dx^2*[ax(2:J)]; M=diag(M1,-1)+diag(M2,0)+diag(M3,1); such that the direct problem can be solved by function U = S(g,u,T) global dt U=u; for dt=diff(T) [M b]=Mb(g,u); du=M\b; u=u+du; U=[U,u]; end The time steps of the numerical method are adapted to the measurment times. They can be refined by the factor 2 by the following mfile: function T = refinetime(T) dThalf=diff(T)/2; taxis=2:2:2*length(T)-2; dT(taxis)=dThalf; dT(taxis-1)=dThalf; Tnew=cumsum(dT); T=Tnew;

Optimization
As result of the numerical simulation, one obtains a grid of computed volume fractions. From this grid one can extract a computed upper interface.

First, an upper and a lower bound on the grid, consisting of two neighboring cells has to be found for each time. The computational interface is yielt by linear interpolation. function h=interface(U,u) [uupper iupper]=max(U-u - (U-u>0)); [ulower ilower]=min(U-u + (U-u<0)); h=X(iupper)+(X(ilower)-X(iupper)).*(0-uupper)./(ulower-uupper); Qdq computes Q and dq is used in the main optimization step function [Q dq]=Qdq(e) global epsilon q=cost(e,0); for K=[1 2], q1(K)=cost(e,K); end for K=[1 1 2 2; 1 2 1 2]; q2(K(1),K(2))=cost(e,K); end dq=(q1-q)/epsilon; Q=(q2'-2*q1'*[1,1]+q)/epsilon^2; % the last computation is a short form for % % Q11=(q2(1,1) - 2*q1(1)+q)/epsilon^2; % Q22=(q2(2,2) - 2*q1(2)+q)/epsilon^2; % Q12=(q2(2,1) - 2*q1(1)+q)/epsilon^2; % Q21=(q2(1,2) - 2*q1(2)+q)/epsilon^2; % % Q=[Q11 Q12; Q21 Q22]; The cost function is evaluated as q=function cost(e,K) global T H u0 parametershift(e,K) U=S('f',T); q=norm(h-H); The small variation over the parameters is done by a parametershift function parametershift(e,K) global uinf C a0 uc k epsilon emask eder=[]; for j=K, eder=[eder, emask(j)]; end I=diag(ones(size(e))); eplus=e+sum(I(eder,:))*epsilon; uinf=eplus(1); C=eplus(2); a0=eplus(3); uc=eplus(4); k=eplus(5);

Testing of direct scheme
The constitutive assumptions can be plotted by global umax uinf C a0 uc k umax=0.66; uinf=1; C=5; a0=1; uc=0.1; k=3; u=(0:0.01:1)'; au=a(u); fu=f(u); subplot(2,1,1); plot(u,au); subplot(2,1,2); plot(u,fu) The maximum of the flux function which is used in the Engquist-Osher numerical flux is computed as global um; [gm im]=max(fu); um=u(im); subplot(2,1,2); hold on; plot(um,gm,'gx') The integral A is calculated by the trapezoidal rule global Au Au=[0;cumsum(([au(1:length(u)-1);0]+[au(2:length(u));0])/2)]; subplot(2,1,2); cla; plot(u,Au) Testing Mb global dt dx; dx=max(H)/5; X=(0:dx:max(H)); u=0.123*ones(size(X')); global J; J=length(u); [M b]=Mb('f',u) Testing S T=[0:0.01:0.04]; U=S('f',u,T); figure(1); hold on; for u=U, plot(u,1:length(u)),pause; end

Testing of optimization scheme
Testing interface . Testing Qdq . Testing cost function .

Testing parametershift global uinf C a0 uc k global epsilon; epsilon=0.1; global emask; emask=[2,5]; e=[1 2 3 4 5]; K=[1 1 2]; parametershift(e,K); [uinf, C, a0, uc, k]

One optimization step
A single optimization step is performed by global epsilon; epsilon=1e-5; global u0; u0=0.123*ones(J,1);

e=[1 2 3 4 5]; [Q dq] = Qdq(e); de=Q\dq; e=e+de; where the above defined subroutines are employed.

Time estimate for reasonable accuracy
Time estimate for reasonable accuracy, where J=100 and T=100, the main loop Mb has to be called 700 times for one optimization step. tic, for k=1:1000; [M b]=Mb('f',u); end; toc gives: "Elapsed time is 14.469154 seconds." This is rather reasonable.