The method of Characteristics for Burgers' equation

Question

I'm trying to solve numerically the inviscid Burgers' equation $u_t + u u_x = 0$ with the method of characteristics. Most of all, I want to see how the numerical solution gets "multiple-values" for times $t>1$ as shown here in figure 3.5.

The initial condition is $u(x,0)=1-\cos(x)$ for $x \in [0,2\pi]$

I wanted to use central finite differences and Euler's method.

The discretized equation becomes $$u'(t)=-u_i(t) \cdot \frac{u_{i+1}(t)-u_{i-1}(t)}{2h} .$$ Now the PDE has become an ODE and I'll use Euler's method, as done in the following runable MatLab code:

mx=100; %number of nodes in x
x=linspace(0,2*pi,mx)';
h=(2*pi)/mx; %step size 
mt=200; %number of time steps
tend=1.5; %final time
k=tend/mt; %step 

%Build the matrix of finite centered differences
B = toeplitz(sparse(1,2,-1/(2*h),1,mx),sparse(1,2,1/(2*h),1,mx));

%initial condition
u0=@(x) 1-cos(x);

u=NaN(mx,mt+1);
u(:,1)=u0(x);
t=0;


for n=1:mt
    u(:,n+1)=u(:,n) - k*(u(:,n).*(B*u(:,n))); %Euler's method
    t=t+k; %update the current time
    plot(x,u(:,n))
    axis([0,2*pi,0,7])
    title(sprintf('t = %0.2f',t));
    xlabel('x')
    ylabel('u(t,x)')
    pause(0.01)
end

This code produces a plot for each time $t$, where I can't see a multi-valued solution! I just see that things go bad for $t>1$, as expected. In fact, for $t \approx 1.1$, I got the following graphic:

Maybe the problem is that I'm using forward Euler, which is not a good choice?

Answer 1

The solution deduced by the method of characteristics becomes multi-valued, as shown in the plot of characteristics below,

and in the linked figure 3.5 of the OP. The code for it can be found in the comments by @LutzL, and a similar case is examined in this post. In facts, discontinuous weak solutions (shock waves) arise at the breaking time $$t_b = -\frac{1}{\min u'_0(x)} = 1 \, .$$ The shock speed is then given by the Rankine-Hugoniot condition, allowing an analytical expression of the shock profile in many cases.

As outlined in the comments by @LutzL, the method proposed here is called method of lines (MOL). Using explicit Euler time-integration, one gets the scheme $$ u_i^{n+1} = u_i^{n} - \frac{\Delta t}{2\Delta x}\, u_i^n (u_{i+1}^n - u_{i}^n) \, , $$ where $u_i^n \approx u(i\Delta x,n\Delta t)$. This method may be unstable due to incorrect upwinding, which is the cause of the oscillations observed here. A similar upwind-biased version of the method is adequate for smooth solutions but will not, in general, converge to a discontinuous weak solution of Burgers' equation as the grid is refined. A classical choice is the Lax-Friedrichs scheme $$ u_i^{n+1} = \frac{1}{2}(u_{i+1}^{n}+u_{i-1}^{n}) - \frac{\Delta t}{2\Delta x} \left(\tfrac{1}{2}(u_{i+1}^n)^2 - \tfrac{1}{2}(u_{i-1}^n)^2\right) , $$ which is stable under the CFL condition $\max u_j^n \frac{\Delta t}{\Delta x} \leq 1$ and converges to the correct shock profile as the grid is refined.

This is implemented in the Matlab code below

mx = 500; % number of nodes in x
CFL = 0.95; % Courant number
g = @(x) (1-cos(x)).*(0<=x).*(x<=2*pi); % initial condition
tend = 1.3; % final time

% initialization
t = 0;
x = linspace(0,2*pi,mx)';
dx = (x(end)-x(1))/(mx-1);
u = g(x);
utemp = u;
dt = CFL*dx/max(u);

figure;
hch = plot(x+t*g(x), g(x), 'k--');
hold on
hlf = plot(x, g(x), 'b.');
xlabel('x');
ylabel('u');
xlim([x(1) x(end)]);
ylim([0 2.1]);
ht = title(sprintf('t = %0.2f',t));

while (t+dt<tend)
    % Lax-Friedrichs
    for i=2:mx-1
        dflux = 0.5*u(i+1)^2 - 0.5*u(i-1)^2;
        utemp(i) = 0.5*(u(i+1) + u(i-1)) - 0.5*dt/dx* dflux;
    end
    utemp(1) = utemp(mx-1);
    utemp(mx) = utemp(2);

    u = utemp;
    t = t + dt;
    dt = CFL*dx/max(u);
    set(hch,'XData',x+t*g(x));
    set(hlf,'YData',u);
    set(ht,'String',sprintf('t = %0.2f',t));
    drawnow;
end
legend('Char.','LF','Location','northwest');

with periodic boundary conditions and the following output:

An interesting feature is the decrease of total energy $E(t) = \int u^2(x,t)\, \text d x$ after the shock formation, as observed numerically in the following picture: