---
tags: HW2.4
---
# ***Computational Hydraulics***
## **Group 01:** N88071053 洪豪男 N88091011 黃彥鈞
## **HW2.4:** (Team Work: 2~3 person for one team) Please apply two of the conventional finite difference methods, as listed in the above table, to solving the traveling wave problem (jump and half-sinus wave), where
(a) the used CFL number should be given
(b) for the half-sinus wave problem, please check if the volume is conserved
(PS either Lax-Wendorff or Beam-Warming method should be included

---

(a-1) : HW2.4 : Traveling Wave Problem ( jump ) : Lax - Wendroff

%% HW2.4 : Traveling Wave Problem ( jump ) : Lax - Wendroff
% Sam Y. J. Huang 2020.10.09

if 01

% initial condition :
% Function : u_t + a u_x = 0 ; a = 1.0
%
x = -100:1:100 ; % range of x
dx = x(2) - x(1) ;
nx = length(x) ;
u = zeros(size(x)) ;
a = 1.0 ;
A = 1.2 ; % amplitude

% Function :

for i = 1:length(u)
if x(i)<=0.0
u(i) = 1.0 ;
end
end

% initial :
figure(1);
plot(x,u,'b--');
xlabel('x');
ylabel('u');
axis([-100 100 -1.4 2]);

u0 = u ; % initial velocity
u_new = zeros(size(u));
stepNo = 0 ;
TotalStep = 300 ;
TotalT = 0.0 ;
CFL = 0.99 ; % CFL condition : max CFL < 1.0

figure(2)
for s = 1:TotalStep
stepNo = stepNo +1 ;
% a = max(abs(u)) ;
dt = CFL*dx/a ;
TotalT = TotalT + dt ;
for j = 2:nx-1 ;
% Lax - Wendroff :
upls = u(j+1) ;
umns = u(j-1) ;
AA = (dt/dx)*u(j) ;
u_new(j) = u(j) - 0.5*AA*(upls - umns) + 0.5*(AA^2)*(upls-2*u(j) + umns) ;

end

% Boundary conditions :
u_new(1) = u_new(nx-1);
u_new(nx) = u_new(2) ;
u = u_new ;

% plot figure :
plot(x,u0,'b--',x,u,'r-')
xlabel('x')
ylabel('u')
axis([-100 100 -1.4 2])
text(40,1.5,['Step No.= ',num2str(stepNo)]);
text(40,1.6,['Time = ',num2str(TotalT)]);
title(' Lax - Wendroff ')
pause(0.01)

end

end

---

(a-2) HW2.4 : Traveling Wave Problem ( jump ) : Leapfrog

%% HW2.4 : Traveling Wave Problem ( jump ) : Leapfrog
% Sam Y. J. Huang 2020.10.09

if 01

% initial condition :
% Function : u_t + a u_x = 0 ;
%
x = -100:1:100 ; % range of x
dx = x(2) - x(1) ;
nx = length(x) ;
u = zeros(size(x)) ;
unmens = zeros(size(x)) ;
a = 1.0 ;
A = 1.2 ; % amplitude

% Function :

for i = 1:length(u)
if x(i)<=0.0
u(i) = 1.0 ;
unmens (i) = 1.0 ;
end
end

u = u + 0.2 ;
unmens = unmens + 0.2 ;

% initial :
figure(1);
plot(x,u,'b--');
xlabel('x');
ylabel('u');
axis([-100 100 -1.4 2]);

u0 = u ; % initial velocity
u_new = zeros(size(u));
stepNo = 0 ;
TotalStep = 100 ;
TotalT = 0.0 ;
CFL = 0.35 ; % CFL condition : max CFL < 1.0

figure(2)
for s = 1:TotalStep
stepNo = stepNo +1 ;
% a = max(abs(u)) ;
dt = CFL*dx/a ;
TotalT = TotalT + dt ;
for j = 2:nx-1 ;
% Leapfrog :

AA = a*(dt/dx) ;
upls = u(j+1) ;
umns = u(j-1) ;
u_new(j) = unmens(j) -0.5*AA*(upls - umns)
end

% Boundary conditions :

u_new(1) = u_new(nx-1);
u_new(nx) = u_new(2) ;
unmens = u ;
u = u_new ;

% plot figure :
plot(x,u0,'b--',x,u,'r-')
xlabel('x')
ylabel('u')
axis([-100 100 -1.4 2])
text(40,1.4,['Step No.= ',num2str(stepNo)]);
text(40,1.5,['Time = ',num2str(TotalT)]);
title(' Leapfrog Method ')
pause(0.01)

end

end

---

(a,b) : HW2.4 : Traveling Wave Problem ( half-sinus wave ) : Lax - Wendroff

%% HW2.4 : Traveling Wave Problem ( half-sinus wave ) : Lax - Wendroff
% Sam Y. J. Huang 2020.10.09

if 01

% initial condition :
% Function : u_t + a u_x = 0 ; a = 1.0
%
x = -100:1:100 ; % range of x
dx = x(2) - x(1) ;
nx = length(x) ;
u = zeros(size(x)) ;
a = 1.0 ;
A = 1.2 ; % amplitude

% Function :

for i = 1:length(u)
if x(i) < 30.0 & x(i) > -30.0
agl = (x(i)+30).*pi/60 ;
tmp = sin(agl)
u(i) = A*tmp ;
else
u(i) = 0.0 ;
end
end

u = u + 0.2 ;
init_V = sum( u*dx ) ; % initial Volume

% initial :
figure(1);
plot(x,u,'b--');
xlabel('x');
ylabel('u');
axis([-100 100 -1.4 2]);

u0 = u ; % initial velocity
u_new = zeros(size(u));
stepNo = 0 ;
TotalStep = 1000 ;
TotalT = 0.0 ;
CFL = 0.2 ; % CFL condition : max CFL < 1.0
new_V = zeros(length(TotalStep)) ;

figure(2)
for s = 1:TotalStep
stepNo = stepNo +1 ;
% a = max(abs(u)) ;
dt = CFL*dx/a ;
TotalT = TotalT + dt ;
for j = 2:nx-1 ;
% Lax - Wendroff :
upls = u(j+1) ;
umns = u(j-1) ;
AA = (dt/dx)*a ;
u_new(j) = u(j) - 0.5*AA*(upls - umns) + 0.5*(AA^2)*(upls-2*u(j) + umns) ;

end

% Boundary conditions :

u_new(1) = u_new(nx-1);
u_new(nx) = u_new(2) ;
u = u_new
new_V(stepNo) = sum(u_new*dx) ;

% plot figure :
plot(x,u0,'b--',x,u,'r-')
xlabel('x')
ylabel('u')
axis([-100 100 -1.4 2])
text(40,1.5,['Step No.= ',num2str(stepNo)]);
text(40,1.6,['Time = ',num2str(TotalT)]);
title(' Lax - Wendroff ')
pause(0.01)

end

figure(3)
plot(1:TotalStep,new_V,'ro', [0 TotalStep],[init_V,init_V],'k-');
xlabel(' Step ')
ylabel(' Volume ')
axis([1 TotalStep min(new_V)-0.5 max(new_V)+0.5]);
legend('New volume','Initial volume')
title(' Volume ')

end

---

(a,b) : HW2.4 : Traveling Wave Problem ( half-sinus wave ) : Leapfrog

%% HW2.4 : Traveling Wave Problem ( half-sinus wave ) : Leapfrog

if 01

% initial condition :
% Function : u_t + a u_x = 0 ;
%
x = -100:1:100 ; % range of x
dx = x(2) - x(1) ;
nx = length(x) ;
u = zeros(size(x)) ;
unmens = zeros(size(x)) ;
a = 1.0 ;
A = 1.2 ; % amplitude

% Function :

for i = 1:length(u)
if x(i) < 30.0 & x(i) > -30.0
agl = (x(i)+30).*pi/60 ;
tmp = sin(agl)
u(i) = A*tmp ;
unmens(i) = A*tmp ;
else
u(i) = 0.0 ;
unmens(i) = 0.0 ;
end
end

u = u + 0.2 ;
init_V = sum( u*dx ) ;
unmens = unmens + 0.2 ;
% initial :
figure(1);
plot(x,u,'b--');
xlabel('x');
ylabel('u');
axis([-100 100 -1.4 2]);

u0 = u ; % initial velocity
u_new = zeros(size(u));
stepNo = 0 ;
TotalStep = 100 ;
TotalT = 0.0 ;
CFL = 0.6 ; % CFL condition : max CFL < 1.0
new_V = zeros(length(TotalStep)) ;

figure(2)
for s = 1:TotalStep
stepNo = stepNo +1 ;
% a = max(abs(u)) ;
dt = CFL*dx/a ;
TotalT = TotalT + dt ;
for j = 2:nx-1 ;

% Leapfrog :
AA = a*(dt/dx) ;
upls = u(j+1) ;
umns = u(j-1) ;
u_new(j) = unmens(j) -0.5*AA*(upls - umns)
end

% Boundary conditions :

u_new(1) = u_new(nx-1);
u_new(nx) = u_new(2) ;
unmens = u ;
u = u_new ;
new_V(stepNo) = sum(u_new*dx) ;

% plot figure :
plot(x,u0,'b--',x,u,'ro')
xlabel('x')
ylabel('u')
axis([-100 100 -1.4 2])
text(40,1.4,['Step No.= ',num2str(stepNo)]);
text(40,1.5,['Time = ',num2str(TotalT)]);
title(' Leapfrog Method ')
pause(0.01)

end

figure(3)
plot(1:TotalStep,new_V,'ro', [0 TotalStep],[init_V,init_V],'k-');
xlabel(' Step ')
ylabel(' Volume ')
axis([1 TotalStep min(new_V)-0.5 max(new_V)+0.5]);
legend('New volume','Initial volume')
title(' Volume ')

end