Numerical method for steady-state solution to viscous Burgers’ equation

finite difference methodsfinite differencesnumerical methodsordinary differential equationspartial differential equations

I am reading a paper in which a specific partial differential equation (PDE)
on the space-time domain $[-1,1]\times[0,\infty)$ is studied. The authors are
interested in the steady-state solution. They design a finite difference method (FDM) for the PDE. As usual, there are certain discretizations in time-space,
$U_j^n$, that approximate the solution $u$ at the mesh points, $u(x_j,t_n)$. The
authors conduct the FDM method on $[-1,1]\times [0,T]$, for $T$ sufficiently large
such that
$$ \left|\frac{U_j^N-U_j^{N-1}}{\Delta t}\right|<10^{-12},\quad \forall j, $$
where $t_N=T$ is the last point in the time mesh and $\Delta t$ is the
distance between the points in the time mesh. The approximations for the steady-state
solution are given by $\{U_j^N\}_j$.

I wonder why the authors rely on the PDE to study the steady-state solution. As
far as I know, the steady-state solution comes from equating the derivatives
with respect to time to $0$ in the PDE. The remaining equation is thus an ordinary
differential equation (ODE) in space. To approximate the steady-state solution,
one just needs to design a FDM for this ODE, which is easier than dealing with
the PDE for sure. Is there anything I am not understanding properly?

For completeness, I am referring to the paper Supersensitivity due to
uncertain boundary conditions. The authors deal with the PDE $u_t+uu_x=\nu u_{xx}$, $x\in (-1,1)$, $u(
-1,t)=1+\delta$, $u(1,t)=-1$, where $\nu,\delta>0$. They employ a FDM for this
PDE for large times until the steady-state is reached. Why not considering the
ODE $uu'=\nu u''$, $u(-1)=1+\delta$, $u(1)=-1$, instead?

Best Answer

You are right. Indeed, solving the ODE for steady-state (time-independent) solutions $(\tfrac12 u^2)' = \nu u''$ is rather straightforward, and we have $$ \tfrac12 u^2 = \nu u' + k_1 \, , $$ where $k_1$ is an integration constant. Thus, integrating the previous ODE gives $$ u(x) = \sqrt{2 k_1} \tanh \left( \sqrt{\frac{k_1}{2}} \frac{k_2 - x}{\nu} \right) , $$ where $k_2$ is an integration constant. The two integration constants are obtained by imposing the boundary conditions, and the sensitivity of the solution with respect to $\delta$ can be investigated.

With such an approach, the thing is that the time information found in the original PDE is lost. In particular, there is no more mapping with initial conditions. In numerical analysis for PDE, steady-state solutions are a good test case for the validation of time-domain numerical methods. Firstly because most of them can be computed analytically, secondly because it is quite challenging for numerical methods to converge towards the equilibrium solution as time becomes infinite. So I believe that validation was the first purpose of the authors. And then, by making the boundary data "uncertain", they remarked a sensitivity of the numerical results with respect to the boundary data.

Related Solutions

[Math] Help deriving Lax-Wendroff scheme for advection equation $u_t+c(x)u_x = 0$

We follow the steps in this post. The first and second-order time derivatives write $u_t = -c(x) u_x$ and \begin{aligned} u_{tt} &= -c(x) u_{tx} \\ & = c(x)\, (c(x) u_x)_x \\ &= c(x)\,( c'(x) u_x + c(x) u_{xx}) \, , \end{aligned} respectively. Using central finite differences in space, we therefore have $$ u_t(x_j,t_n) \approx -c_j \frac{D_x u_j^n}{2 \Delta x} $$ and $$ u_{tt}(x_j,t_n) \approx c_j \left(\frac{D_x c_j}{2 \Delta x}\frac{D_x u_j^n}{2\Delta x} + c_j \frac{\delta^2_x u_j^n}{\Delta x^2}\right) . $$ Injecting this Ansatz in the Taylor series $$ u(x_j,t_{n+1}) = u(x_j,t_{n}) + \Delta t\, u_t(x_j,t_{n}) + \frac{\Delta t^2}{2} u_{tt}(x_j,t_{n}) +\dots $$ gives the proposed numerical method.

Partial Differential Equations – Method of Characteristics for Burgers’ Equation

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:

Best Answer

Related Solutions

[Math] Help deriving Lax-Wendroff scheme for advection equation $u_t+c(x)u_x = 0$

Partial Differential Equations – Method of Characteristics for Burgers’ Equation

Related Question