TVD numerical solution of Burgers' equation with sinusoidal initial data

Question

I have solved Burgers' equation $u_t + u u_x = 0$ using a Total Variation Diminishing method for the initial conditions $u(x,0) = \sin(2\pi x)$. I have used forward Euler time-integration with CFL of 0.2. The TVD method is MUSCL with Kurganov-Tadmor central scheme and Van Leer limiter. The results are as shown here:

My questions about it:

• Is the CFL of 0.2 very small?
• I was getting oscillations for higher CFL, though its a TVD scheme. That's weird right? Do the results look okay?
• How would I find the exact solution? I am aware of this post, but how can I do error estimate?

Answer 1

Let us derive the Kurganov-Tadmor central scheme for Burgers' equation $u_t + f(u)_x = 0$ which physical flux $f(u) = \frac{1}{2}u^2$ is convex. The numerical method is written in semi-discrete form (see Eq. (4.2) of (1)) $$ \frac{\text d u_i}{\text d t} = -\frac{H_{i+1/2} - H_{i-1/2}}{\Delta x} , $$ where the numerical flux reads $$ \begin{aligned} H_{i+1/2} &= \frac{1}{2} \left(f(u^L_{i+1/2}) + f(u^R_{i+1/2}) - a_{i+1/2}\, (u^R_{i+1/2}-u^L_{i+1/2})\right) ,\\ a_{i+1/2} &= \max \left\lbrace |f'(u^L_{i+1/2})|, |f'(u^R_{i+1/2})|\right\rbrace . \end{aligned} $$ The slope-limited extrapolated interface values of $u$ are given by $$ \begin{aligned} u_{i+1/2}^{L} &= u_i + \frac{\Delta x}{2} (u_x)_{i} ,\\ u_{i+1/2}^{R} &= u_{i+1} - \frac{\Delta x}{2} (u_x)_{i+1} ,\\ (u_x)_{i} &= \text{minmod}\left(\frac{u_i-u_{i-1}}{\Delta x}, \frac{u_{i+1}-u_{i}}{\Delta x}\right) , \end{aligned} $$ where the minmod limiter function is $(a,b)\mapsto \frac{1}{2}(\text{sign}\, a + \text{sign}\, b)\min(|a|,|b|)$. Numerical results obtained with the following piece of Matlab code are shown below. Here, second-order Runge-Kutta integration is used. The Courant number is set to $\text{Co} = 0.2$, and periodic boundary conditions are implemented. Numerically, it seems that TV-stability is even preserved for larger Courant numbers, e.g. $\text{Co} = 0.9$, but this is no longer true if forward Euler time-integration is used instead of the present RK2 method. Note that the modification $(u_x)_i = 0$ of the method yields the local Lax-Friedrichs (a.k.a Rusanov) method if forward Euler time-integration is used.

function y = minmod(a,b)
    y = 0.5*(sign(a)+sign(b)).*min(abs(a),abs(b));
end

function [y,yp] = f(u)
    y = 0.5*u.^2;
    yp = u;
end

function y = RHS(u,dx)
    ux = minmod((u-circshift(u,[0 1]))/dx,(circshift(u,[0 -1])-u)/dx);
    uL = circshift(u + 0.5*dx*ux,[0 1]);
    uR = u - 0.5*dx*ux;
    [fL,fpL] = f(uL);
    [fR,fpR] = f(uR);
    a = max(abs(fpL),abs(fpR));
    H = 0.5*(fL + fR - a.*(uR-uL));
    y = -(circshift(H,[0 -1]) - H)/dx;
end

%%%%%%%%%%%%%% main program %%%%%%%%%%%%%%
%parameters
Nx = 200;
Co = 0.2;
tmax = 2;

%initialization
x = linspace(0,1,Nx);
t = 0;
u = sin(2*pi*x);
dx = x(2) - x(1);
dt = Co*dx/max(abs(u));

figure(1);
clf;
plot(x,u);
hold on
h = plot(x,u,'k.');

%iterations
while (t<tmax)
    u1 = u + dt*RHS(u,dx);
    u = 0.5*u + 0.5*(u1 + dt*RHS(u1,dx));

    set(h,'YData',u);
    drawnow;

    dt = Co*dx/max(abs(u));
    t  = t + dt;
end

The theoretical solution can be obtained quasi-analytically by following the steps in this post. Indeed, a static shock is located at $x=0.5$, and the solution on each side can be deduced from the method of characteristics. Hence, we need to solve numerically the equation $u = \sin (2\pi (x - ut))$ to find the value of $u(x,t)$ at $x\neq 0.5$, which can be done by using root finding methods:

fun = @(u) u - sin(2*pi*(x-u*t));
u0 = x/t.*(x<0.5) + (x-1)/t.*(0.5<x);
uth = fsolve(fun,u0);
plot(x,uth,'r');
disp(norm(uth-u));

Repeating these steps for various mesh sizes $\Delta x$ leads to error measurements. However, it should be noted that we are limited by the precision of the root finding method. Lastly, the order of convergence equals $\approx 0.5$ when the solution is discontinuous (i.e., when $t>\frac{1}{2\pi}$ is larger than the breaking time).

---

(1) A. Kurganov, E. Tadmor (2000): "New High-Resolution Central Schemes for Nonlinear Conservation Laws and Convection-Diffusion Equations", J. Comput. Phys. 160(1), 241–282. doi:10.1006/jcph.2000.6459