54 Computational Modelling in Hydraulic and Coastal Engineering
mt = 200;
% Initialization of the solution field;
for k=2:mx
f1(k)=0; f1n(k)=0;
f2(k)=0; f2n(k)=0;
f3(k)=0.0001; f3n(k)=0.0001;
end
% Upstream boundary conditions;
f1(1)=1; f1n(1)=1;
f2(1)=1; f2n(1)=1;
f2(2)=1; f2n(2)=1;
f3(1)=1; f3n(1)=1;
% Stability criterion (sc<1);
sc = C*Dt/Dx;
% Integration by the Godunov F.D. scheme;
for n = 1:mt
for n1 = 2:mx-1
f1n(n1)=f1(n1)-sc*(f1(n1)-f1(n1-1));
end
% Boundary conditions;
f1(1)=f1n(1); f1n(mx)=f1n(mx-1)*2-f1n(mx-2);
% Renewal of the f1 values
for n1=1:mx
f1(n1)=f1n(n1);
end
f1p=f1(1:n1);
% Integration by the Fromm F.D. scheme;
for n2=3:mx-1
f2ni(n2)=f2(n2)-(sc/4)*(f2(n2+1)-f2(n2-1)+f2(n2)f2(n2-2));
f2n(n2)=f2ni(n2)+(sc/2)^2*(f2(n2+1)-2*f2(n2)+f2(n21))+((C*Dt)^2-2*C*Dt*Dx)/
(4*Dx^2)*(f2(n2-2)-2*f2(n2-1)+f2(n2));
end
% Renewal of f2 values;
for n2=1:mx
f2(n2)=f2n(n2);
end
% Boundary conditions;
f2n(mx)=f2n(mx-1)*2-f2n(mx-2);
f2(1)=f2n(1); f2(mx)=f2n(mx);
f2p=f2(1:n2);
% Integration by the TVD scheme
% Estimation of the fr parameter;
for n3=2:mx-1
dd=f3(n3+1)-f3(n3);
if dd == 0
r1=-1;
else
r1=(f3(n3)-f3(n3-1))/dd;