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:
The method of characteristics gives the implicit equation $u = \phi(x-ut)$. Therefore, if $x-ut\leq 0$, we have $u=1$. If $x-ut\geq 1$, we have $u=0$. Otherwise, we have $u = \frac{1-x}{1-t}$. To summarize the solution obtained via the method of characteristics,
$$
u(x,t) =
\left\lbrace
\begin{aligned}
&1 && x\leq t\\
&\tfrac{1-x}{1-t} && t<x<1\\
&0 && 1\leq x
\end{aligned}
\right.
$$
which is valid for $t<1$. At $t=1$ (breaking time), characteristics intersect and a shock wave occurs. This is illustrated on the following plot of the characteristic curves in the $x$-$t$ plane
See this post for a book reference.
Best Answer
The method of characteristics gives $u=f(x-ut)$ where $f(x) = \sin(x)$ represents the initial data $u(x,0)$. This implicit equation can be used to approach the solution, until a shock wave is met by the characteristics. Here is a sketch of the characteristics in the $x$-$t$ plane:
As shown in this post, the shock is formed when characteristics intersect for the first time, i.e. at the breaking position and time \begin{aligned} x_B &= \arg \min f'(x) = (2n+1)\pi , \qquad (n\in \Bbb{Z})\\ t_B &= \frac{-1}{\min f'(x)} =1 \end{aligned} as shown in the figure. Setting $v=-u$, one notes the relationship $$ v(x_B - x, t) = \sin \left[x - v(x_B - x, t)\, t\right] $$ which proves the symmetry property $u(x_B - x, t) = -u(x, t)$.
By solving the implicit equation $u = f(x-ut)$, the method of characteristics yields the values $u_L$, $u_R$ of $u$ on the left and the right of the shock. The position $(x_s(t), t)$ of the shock wave is given by the Rankine-Hugoniot condition $$x'_s(t) = \tfrac{1}{2} (u_L+u_R)$$ with initial position $x_s(t_B) = x_B$. Therefore, using symmetry $u_L = -u_R$, the shock wave remains at the same place (static shock). Its magnitude $u_L$ increases and then vanishes in time, as shown in the figure below, where $u_L$ --obtained numerically-- is plotted with respect to $t$: