---
# System prepended metadata

title: '***Computational Hydraulics***'
tags: [HW2.4]

---

---
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

