Inviscid Burgers equation with trapezoidal boundary data

Question

Im trying to solve $v_t + vv_x = 0$ subject to

$$ v(x,0) = \begin{cases} 0, & 0 \leq x \leq 1 \\ x -1, & 1 \leq x < 2 \\ 1, & 2 \leq x \leq 3 \\ 4 - x, & 3 \leq x \leq 4 \\ 0, &4 \leq x\leq 5 \end{cases} $$

and $v(0,t)=v(5,t)=0$. So, the initial condition is a trapezoid looking function.

We see that we have rarefraction at $x=1$ and $x=4$ and shocks at $x=2,3$. Im trying to find the exact solution just for $0< t \leq 2$, but even in this time interval, it seems a little bit laborious to compute the solutions since the shock waves will intersect with rarefaction waves and so on.

What is the best approach to compute the exact solution? Also, I would like some explanation as to how we could implement godunov scheme in matlab in this situation.

Answer 1

Let us plot the characteristic curves deduced from the method of characteristics. The latter are lines in the $x$-$t$ plane, along which $v$ is constant:

One observes that the curves intersect at the breaking time $t_b = -1/\inf v_x(x,0) = 1$. Before the breaking time, $0 \leq t < 1$, the solution deduced from the method of characteristics reads $$ v(x,t) = \left\lbrace \begin{aligned} &0 & & 0\leq x \leq 1\\ &\tfrac{x-1}{1+t} & & 1\leq x \leq 2+t\\ &1 & & 2+t\leq x \leq 3+t\\ &\tfrac{4-x}{1-t} & & 3+t\leq x \leq 4\\ &0 & & 4\leq x \leq 5\\ \end{aligned} \right. $$ The shock wave generated at $t=1$ has left state $v_l=1$ and right state $v_r=0$. Therefore, the speed of shock deduced from the Rankine-Hugoniot condition is $s = 1/2$. The solution for $t\geq t_b$ reads $$ v(x,t) = \left\lbrace \begin{aligned} &0 & & 0\leq x \leq 1\\ &\tfrac{x-1}{1+t} & & 1\leq x \leq 2+t\\ &1 & & 2+t\leq x \leq (7+t)/2\\ &0 & & (7+t)/2\leq x \leq 5\\ \end{aligned} \right. $$ This solution is maximally valid until $2+t = (7+t)/2$ or $(7+t)/2 = 5$, i.e., $1\leq t<3$.

The Godunov scheme is coded as usual for Burgers' equation, only the initial/boundary conditions must be implemented. Godunov's method is written in conservation form as (see Chap. 12 of (1)) $$ u_i^{n+1} = u_i^n - \frac{\Delta t}{\Delta x}(f_{i+1/2}^n - f_{i-1/2}^n) , $$ with the numerical flux $$ f_{i+1/2}^n = \left\lbrace \begin{aligned} &\tfrac{1}{2}(u_i^n)^2 & &\text{if } u_i^n > 0 \text{ and } \tfrac{1}{2}(u_i^n+u_{i+1}^n) > 0 , \\ &\tfrac{1}{2}(u_{i+1}^n)^2 & & \text{if } u_{i+1}^n < 0 \text{ and } \tfrac{1}{2}(u_i^n+u_{i+1}^n) < 0 , \\ &0 & & \text{if } u_i^n < 0 < u_{i+1}^n . \end{aligned}\right. $$ The initial condition is implemented by a proper initialization of the data vector $(u_i^0)_{0\leq i\leq N_x}$. The boundary conditions are specified in the ghost cells by setting $u_{-k}^n = 0$ and $u_{N_x+k}^n = 0$ for $k \geq 1$ at every time step. A Matlab implementation and its output are provided below.

% numerics
Nx = 80;   % number of points
Co = 0.95;  % Courant number
tmax = 1.5;

% analytical solution for t<3
vth = @(x,t) 0*x + (x-1)./(1+t).*(x>=1).*(x<=2+t)   ...
  + 1.*(x>=2+t).*(x<=min(3+t,(7+t)/2))              ...
  + (4-x)./(1-t).*(x>=3+t).*(x<=4);

% initialization
t = 0;
x = linspace(0,5,Nx);
dx = x(2)-x(1);
x = [x(1)-dx x x(end)+dx];
u = vth(x,t);
f = zeros(1,Nx+1);
for i=1:Nx+1
    s = mean(u(i:i+1));
    f(i) = 0.5*u(i)^2*(u(i)>0)*(s>0) + 0.5*u(i+1)^2*(u(i+1)<0)*(s<0);
end
dt = Co*dx/max(abs(u));

% graphics
figure(1);
xth = linspace(0,5,400);
pth  = plot(xth,vth(xth,t),'k-','LineWidth',2);
hold on
pnum = plot(x,u,'b.-','LineWidth',1);
xlim([0 5]);
ylim([-0.5 1.5]);
ptit = title(strcat('t = ',num2str(t)));
xlabel('x');
ylabel('v')

% iterations
while (t+dt<tmax)
    u(2:end-1) = u(2:end-1) - dt/dx*(f(2:end) - f(1:end-1));
    for i=1:Nx+1
        s = mean(u(i:i+1));
        f(i) = 0.5*u(i)^2*(u(i)>0)*(s>0) + 0.5*u(i+1)^2*(u(i+1)<0)*(s<0);
    end
    dt = Co*dx/max(abs(u));
    t = t + dt;
    set(pth,'YData',vth(xth,t));
    set(pnum,'YData',u);
    set(ptit,'String',strcat('t = ',num2str(t)));
    drawnow;
end

---

(1) R.J. LeVeque, Finite Volume Methods for Hyperbolic Problems, Cambridge University Press, 2002.