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?
Best Answer
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
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: