# Solving Burger's equation with NDSolve at large time

I want to solve the following Burger's equation $$\partial_tu+\partial_x\left(\frac{u^2}{2}\right)=0,~~x\in[0,2\pi],~t>0\\u(x,0)=\frac{1}{3}+\frac{2}{3}\sin(x)$$ with mathematica. Here's my code:

Clear["Global*"];
{SOL} = NDSolve[{D[u[x, t], t] + u[x, t] D[u[x, t], x] == 0,
u[x, 0] == 1/3 + 2 Sin[x]/3}, u, {x, 0, 2 Pi}, {t, 0, 2 Pi},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MaxStepSize" -> 2 Pi/4000}}];
Plot[u[x, 1.5] /. SOL, {x, 0, 2 Pi}, AxesLabel -> Automatic]


The result looks good at time $$t=1.5$$:

I wanna solve this equation numerically for larger $$t$$, for example $$t=2\pi$$, however when I plot $$t=2\pi$$ everything goes wrong, Obviously the exact solution cannot be like that, it seems the above code is no longer working for large $$t$$, but I don't know how to edit the code to make the numerical results reasonable. May I ask what should I add/change to the code to solve the results for $$t=2\pi$$? Thanks in advance!

Edit: If we don't use FEM, we can try some other ways, such as the Godunov's scheme, and the code goes like below (The ExactSol in the code is the solution by using other methods, here I used mesh size $$\Delta x = \frac{2\pi}{10} ,\frac{2\pi}{20} , \cdots\frac{2\pi}{640}$$ to show that this is actually working~)

color[1] = Red; color[2] = Orange; color[3] = Yellow; color[4] = Green; color[5] = Cyan; color[6] = Blue; color[7] = Purple;
For[o = 1, o <= 7, o++,
f[u_] = u^2/2; T = 2 Pi; h = 0.4 Pi/2^o; k = 2 Pi/1280;
Godunov[ul_, ur_] = If[ul > ur, Max[f[ul], f[ur]], If[ul <= 0 <= ur, 0, Min[f[ul], f[ur]]]];
For[m = -2000, m <= 2000, m++, x[m] = h m; u[m, 0] = 1/3 + 2 Sin[x[m]]/3];
For[n = 0, n < 1280, n++,
For[m = -1999 + n, m < 1999 - n, m++,
u[m, n + 1] = N[u[m, n] - k*(Godunov[u[m, n], u[m + 1, n]] - Godunov[u[m - 1, n], u[m, n]])/h];
];
];
s[o] = ListLinePlot[Table[u[m, 1280], {m, 0, 5*2^o}], DataRange -> {0, 2 Pi}, PlotStyle -> {color[o], Thickness[0.001]}, Frame -> True, PlotLegends -> {5*2^o "Grids Numerical Solution"}, FrameTicks -> {{{-0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7}, True}, {{0, 0.5 Pi, Pi, 1.5 Pi, 2 Pi}, True}}];
];
Show[s[1], s[2], s[3], s[4], s[5], s[6], s[7], Plot[ExactSol[x, 2 Pi], {x, 0, 2 Pi}, PlotRange -> {-0.5, 1}, PlotStyle -> {Black, Thickness[0.0015]}, Frame -> True, PlotLegends -> {"Exact Solution"}, FrameTicks -> {{{-0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7}, True}, {{0, 0.5 Pi, Pi, 1.5 Pi, 2 Pi}, True}}], PlotRange -> All]


and the result looks nice too :)

Similar problems have been discussed in this site for several times:

Solution of Burgers equation with some initial data

1D Euler equations (fluid dynamics) with NDSolve

Circumvent NDSolve::bdord: error for 1-D Euler Equations

In short, PDE related to hyperbolic conservation law whose solution involves shock wave is not well handled by NDSolveat the moment. We can of course add artificial viscosity (as shown in N0va's answer and linked questions above):

mol[n:_Integer|{_Integer..}, o_:"Pseudospectral"] := {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n,
"MinPoints" -> n, "DifferenceOrder" -> o}}

{ufunc, vfunc} =
NDSolveValue[{With[{μ = 0.001, u = u[x, t], v = v[x, t]},
{D[u, t] + 1/2 D[v, x] == μ D[u, x, x],
v == u^2,
{v == (1/3 + 2 Sin[x]/3)^2,
u == 1/3 + 2 Sin[x]/3} /. t -> 0}],
u[0, t] == u[2 Pi, t]}, {u, v}, {x, 0, 2 Pi}, {t, 0, 20},
Method -> mol[600, 8]]
(* {1.94787, Null} *)


but as we can see, it's rather tricky.

Remark

I've introduced the intermediate variable v to make the solution respect the conservation law, otherwise the obtained solution will probably be undesired. See here for more info.

So I'd like to show one more solution based on Julia package Trixi. Thanks to ExternalEvaluate, nowadays we can easily call Julia inside Mathematica (Remember to configure your Julia and install the needed packages first):

session = StartExternalSession["Julia"];
eval = ExternalEvaluate[session, #] &;

eval@"using OrdinaryDiffEq, Trixi"; // AbsoluteTiming
(* {33.4249, Null} *)

sol = eval@"
equations = InviscidBurgersEquation1D()
ic(x, t, equation::InviscidBurgersEquation1D)=1/3 + 2/3*sin(x[1])

volume_flux = flux_lax_friedrichs
solver = DGSEM(polydeg=3,surface_flux=flux_lax_friedrichs,
volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

xmin = 0
xmax = 2*pi
mesh = TreeMesh(xmin, xmax,
initial_refinement_level=8,
n_cells_max=10^4)
semi = SemidiscretizationHyperbolic(mesh, equations, ic, solver)

tspan = (0, 20)
ode = semidiscretize(semi, tspan)
sol = solve(ode,  SSPRK43())"; // AbsoluteTiming
(* {18.8699, Null} *)

xcoord =
eval@"collect(Iterators.flatten(semi.cache.elements.node_coordinates))";
tcoord = eval@"sol.t";
solfunc = ListInterpolation[sol[[;; ;; 4]], Flatten /@ {xcoord[[;; ;; 4]], tcoord}
(*, InterpolationOrder -> 1*)];
Manipulate[Plot[solfunc[x, t], {x, 0, 2 Pi}, PlotRange -> {-0.5, 1.5}],
{t, 0, 2 Pi}]


Notice b.c. isn't explicitly set in Julia code, because the default setting of SemidiscretizationHyperbolic is periodic b.c., which is exactly what we need in this case.

To understand the Julia code, you may want to read

https://github.com/trixi-framework/Trixi.jl/issues/939

https://github.com/trixi-framework/Trixi.jl/issues/958

Finally a comparison of the two solutions above:

Plot[{ufunc[x, 2 Pi], solfunc[x, 2 Pi]}, {x, 0, 2 Pi},
PlotStyle -> {Automatic, Dashed}]


## Update

In 1st solution above TensorProductGrid sub-method is used for spatial discretization. One can use FiniteElement method instead. The advantage of FiniteElement is we don't need to introduce intermediate variable like above (maybe because the ODE (DAE?) solver has done a better job in this case), but the disadvantage is, the behavior of periodic b.c. of FiniteElement is undesired. This issue is discussed in this post. To circumvent this, we need to introduce a ghost layer and set 2 PeriodicBoundaryCondition:

molfem[measure_: Automatic] := {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> MaxCellMeasure -> measure}};

With[{μ = 0.001, eps = 1/1000},
solfem = NDSolveValue[{D[u[x, t], t] +
1/2 D[u[x, t]^2, x] == μ D[u[x, t], x, x], u[x, 0] == 1/3 + 2 Sin[x]/3,
PeriodicBoundaryCondition[u[x, t], x == 0, TranslationTransform[{2 Pi}]],
PeriodicBoundaryCondition[u[x, y], x == 2 Pi + eps,
TranslationTransform[{-2 Pi}]]}, u, {x, 0, 2 Pi + eps}, {t, 0, 20},
Method -> molfem[0.01]]]; // AbsoluteTiming
(* {34.2291, Null} *)


Let's compare the wave speed at $$t=20$$.

With[{μ = 0.001},
solundesired =
NDSolveValue[{D[u[x, t], t] + 1/2 D[u[x, t]^2, x] == μ D[u[x, t], x, x],
u[x, 0] == 1/3 + 2 Sin[x]/3, u[0, t] == u[2 Pi, t]},
u, {x, 0, 2 Pi}, {t, 0, 20},
Method -> mol[1000, 8]]]; // AbsoluteTiming
(* {1.78104, Null} *)

Plot[{ufunc[x, 20], solfem[x, 20], solfunc[x, 20], solundesired[x, 20]},
{x, 0, 2 Pi}, PlotStyle -> {Automatic, Dashed, DotDashed}]


• Nice answer, very interesting. I tried to find a Finite Element solution, which shows a similar wave moving to the right. But this solution keeps the boundary x==0 constant, probably due to implicit flux condition (NeumannValue). I tried to introduce periodic boundary conditions(TensorProductGrid sets this bc type automatically ? ) but didn't succeed. Any idea how to find a FEM solution? Thanks! May 8 at 7:06
• @UlrichNeumann Yeah, there's a issue about PeriodicBoundaryCondition. Actually I learned the solution from a question asked by you :D . See my update. May 8 at 8:55
• Thank you for your fast reply and the hint concerning my own question ;-). Could you please add molfem to your answer? May 8 at 9:10
• @UlrichNeumann Oops, sorry, I forgot to add its definition. Revised. Thx for pointing out. May 8 at 9:12
• @xzczd, as an alternative to the ghost layer, it might work to use a triangle based mesh, try the option "MeshElementType"->"TriangleElement", this way you may get rid of the ghost layer. May 9 at 12:18

The answer of Nasser is only in parts correct. The analytic/implicit solution $$u(x,t) = f(\xi) = f(x-ut), \quad \xi = x - f(\xi) t$$ is only the physical solution in the interval $$0\le t \le t_b$$ with $$t_b = \inf_x\left(\frac{-1}{ f'(x)}\right),$$ for more details consult the already mentioned wikipedia article. For the initial condition at hand $$t_b=3/2$$ and it is no coincidence that NDSolve fails at $$t\approx 1.5$$ since at this point a shockwave is formed at $$x\approx3.6$$. The physical solution to the Burger's equation for $$t>t_b$$ is no longer smooth/analytic and only exists in a weak sense. Naive discretization schemes (like finite differences) are notorious for being ill suited for problems like the invicit Burgers equation. I do not know if there are options for NDSolve which can cope with shock formation but one way to weaken the numerical difficulty without changing the qualitative features of the physical solution too much is to introduce a small viscosity term and to study the viscous Burgers equation $$\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} = \nu\frac{\partial^2 u}{\partial x^2},$$ with the viscosity $$\nu>0$$. The introduced diffusion smears out the shockwave and makes the solution differentiable this combined with a high resolution allows the "MethodOfLines" option for NDSolve to work. I included proper (periodic) boundary condition to get a well defined problem:

eqs={D[u[x,t],t]+u[x,t] D[u[x,t],x]==\[Nu]*D[u[x,t],x,x],u[x,0]==1/3+2 Sin[x]/3,u[0,t]==u[2\[Pi],t]}


With small diffusion ($$\nu=0.01$$) it is possible with NDSolve to compute a solution after shock formation:

{sol}=NDSolve[eqs/.\[Nu]->0.01,u,{x,0,2 Pi},{t,0,2},Method->{"MethodOfLines","SpatialDiscretization"->{"TensorProductGrid","MinPoints"->1000, "DifferenceOrder"->"Pseudospectral"}},AccuracyGoal->10,PrecisionGoal->10];

ti={0,1,1.5,2.};
Table[u[x,t],{t,ti}]/.sol;
Plot[%,{x,0,2 Pi},FrameLabel->{"x","u[x,t]"},Frame->True,PlotRange->{{0,2\[Pi]},{-0.5,1.5}},PlotLegends->Placed[Row[{"t=",#}]&/@ti,{Right,Top}]]


without diffusion ($$\nu=0$$) it is not beyond $$t_b=1.5$$:

{sol}=NDSolve[eqs/.\[Nu]->0.0,u,{x,0,2 Pi},{t,0,2},Method->{"MethodOfLines","SpatialDiscretization"->{"TensorProductGrid","MinPoints"->1000, "DifferenceOrder"->"Pseudospectral"}},AccuracyGoal->10,PrecisionGoal->10];

ti={0,1,1.5,1.6};
Table[u[x,t],{t,ti}]/.sol;
Plot[%,{x,0,2 Pi},FrameLabel->{"x","u[x,t]"},Frame->True,PlotRange->{{0,2\[Pi]},{-0.5,1.5}},PlotLegends->Placed[Row[{"t=",#}]&/@ti,{Right,Top}]]


The solution to the PDE beyond $$t_b$$ is no longer naively differentiable after the formation of a shock (for $$\nu=0$$) which gives rise to the spurious oscillations we can see in the plot above at $$t=1.6$$. As time progresses the magnitude of those oscillations increase and ultimately destabilize the numerical scheme.

A proper way to solve the invicit Burger's equation numerically would be to deploy shock capturing numerical schemes for conservation equations like specific finite volume or finite element methods which work with a weak formulation and which compute weak, physical solutions to this kind of PDEs.

• Still, this isn't the end. The wave speed in this case is probably undesired. See: nbviewer.org/github/numerical-mathematics/… Turn to FiniteElement method will help, because its formal form is based on the conservation law, see e.g. mathematica.stackexchange.com/a/197305/1871 May 6 at 10:47
• I have to retract the line "Turn to FiniteElement method will help, because its formal form is based on the conservation law", see discussion under mathematica.stackexchange.com/q/267966/1871. But it's still true that when the artificial viscosity is small enough and TensorProductGrid is chosen for spatial discretization and the conservation law isn't respected, the wave speed can easily go wrong. My current guess is, this is related to the ODE (DAE?) solver. May 10 at 5:46

Following the example set by xzczd's answer: an alternative to an external library would be to implement a simple shock capturing finite volume scheme like the one put forward in New High-Resolution Central Schemes for Nonlinear Conservation Laws and Convection–Diffusion Equations, 2000, A. Kurganov and E. Tadmor -- KT FV scheme -- see also this wikipedia article. I have implemented this scheme in Mathematica a few years ago. The scheme is purpose build for hyperbolic conservation laws and includes extensions for the treatment of parabolic conservation laws (e.g. PDEs including diffusion terms). The scheme can be applied directly to non-linear PDEs (initial value problems) of the type \begin{align} \partial_t u(t,x)+\partial_x F(t,x,u)&=\partial_x Q(t,x,u,\partial_x u)\\ u(0,x)&=u_0(x) \end{align} with a broad range of boundary conditions (ones which can be implemented easily with ghost-cells). Below you find an implementation of the scheme using CompiledFunction for performance.

(** Periodic boundary conditions for ghost cells **)
BCperiodic=Function[{vin},Module[{v},
v=vin;
v[[1]]=v[[-5]];(* Subscript[u, -2 ]= Subscript[u, n-3] *)
v[[2]]=v[[-4]];(* Subscript[u, -1 ]= Subscript[u, n-2] *)
v[[-3]]=v[[3]];(* Subscript[u, n-1 ]= Subscript[u, 0] *)
v[[-2]]=v[[4]];(* Subscript[u, n] = Subscript[u, 1] *)
v[[-1]]=v[[5]];(* Subscript[u, n+1] = Subscript[u, 2] *)
v
]];

(** Selected MUSCL scheme Flux limiter **)
ϕMinMod=Function[{Δuj,Δujp1},If[Abs@Δujp1==0.,1.,Max[0.,Min[1.,Δuj/Δujp1]]]]; (* MinMod Limiter [https://en.wikipedia.org/wiki/Flux_limiter] *)

(** Main step method **)
ClearAll[KTstepper];
KTstepper[nx_][Δx_,x_,x12_][F_,dFdu_,Q_,BC_,pars___]:=
Compile[{{t,_Real},{v,_Real,1}},
Module[{
u,(* Conserved quanties, dim=nx+4 *)
Δu,(* Differences of conserved quanties, dim=nx+3 *)

ux,(* Reconstructed slopes, dim =nx+2 *)
up12T,(* Right intermediate values at the cell interface, dim =nx+1 *)
um12T,(* Left intermediate values at the cell interfaces, dim=nx+1 *)
Fp12T,(* Advection flux through the cell interfaces, dim=nx+1 *)
Fm12T,(* Advection flux through the cell interfaces, dim=nx+1 *)
λp12T, (* Jacobian eigenvalues (single value) the cell interfaces, dim=nx+1 *)
λm12T, (* Jacobian eigenvalues (single value)the cell interfaces, dim=nx+1 *)
a12T,(* Approxmiate velocities at the cell interfaces, dim=nx+1 *)
H12 ,(* Numerical advection fluxes at the cell interfaces, dim=nx+1 *)

P12(* Diffusion fluxe through the cell interfaces, dim =nx+1 *)
},
(* **** Boundary condition **** *)
u=BC[Join[{0.,0.},Take[v,{1,nx}],{0.,0.}]];

(* **** PDE Convection Flux **** *)
Δu=Differences[u];
ux=MapThread[0.5*#2(ϕMinMod[#1,#2])&,{Take[Δu,{1,-2}],Take[Δu,{2,-1}]}];(* [KTO2-0, eq. (2.4)*0.5*Δx]: modified to be compatible with generic flux limiters *)
up12T=Take[u,{3,-2}]-Take[ux,{2,-1}];(* [KTO2-0, eq. (4.5)]: modified to be compatible with generic flux limiters *)
um12T=Take[u,{2,-3}]+Take[ux,{1,-2}];(* [KTO2-0, eq. (4.5)]: modified to be compatible with generic flux limiters *)

a12T=MapThread[Max[Max@Abs@#1,Max@Abs@#2]&,{λp12T,λm12T}]; (* [KTO2-0, eq. (3.2)] [KTO2-0, footnote 2] *)

H12=0.5*(Fp12T+Fm12T-a12T*(up12T-um12T)); (* [KTO2-0, eq. (4.4)] *)

(* **** PDE Diffusion Flux **** *)
Join[{x12[[1]]-Δx*0.5},x],
Join[x,{x12[[-1]]+Δx*0.5}],
Take[u,{2,-3}],
Take[u,{3,-2}],
Take[Δu,{2,-2}]/Δx
}]; (* [KTO2-0, eq. (4.14)] *)

(* **** Result **** *)
Return[Flatten[(-Differences[H12]+Differences[P12])/Δx,1]] (* [KTO2-0,eq.(4.13)] *)
], RuntimeOptions->"Speed", CompilationOptions->{"InlineExternalDefinitions"->True,"ExpressionOptimization"->True}
]

(** Solver/time-stepper using NDSolve to solve the ODE (MOL) system of KTstepper *)
ClearAll[KTsolver]
KTsolver/:Options[KTsolver]={AccuracyGoal->Automatic,PrecisionGoal->Automatic,WorkingPrecision->MachinePrecision,Method->Automatic,MaxSteps->Automatic};
KTsolver[OptionsPattern[]][x0_,x1_,n_Integer][F_,dFdu_,Q_,BC_,params___][u0_][t0_,t1_]:=Module[{
no,
Δx,xi,xi12,
v0,
PDEcompileFkt,PDEfkt,PDEsystem,PDEv,PDEt,PDEsolver,PDEsolution,
nt,tw,tm},

(* x0 and x1 are the first and last cell centers *)
Δx=(x1-x0)/(n-1);
xi=Table[x0+Δx*i,{i,0,n-1}];
xi12=Table[x0+Δx*(i-0.5),{i,0,n}];

v0=u0[#]&/@xi;(* approimate cell averages with mid point value *)

PDEcompileFkt=KTstepper[n][Δx,xi,xi12][F,dFdu,Q,BC,params];
PDEfkt[t_?NumericQ,v_]:=PDEcompileFkt[t,v];
PDEsystem={Equal[PDEv[t0],v0],Equal[Derivative[1][PDEv][PDEt],PDEfkt[PDEt,PDEv[PDEt]]]};

nt=0;
tw=-AbsoluteTime[];
PDEsolver=Inactive[NDSolveValue][PDEsystem,PDEv,{PDEt,t0,t1},
StepMonitor:>{nt++,tm=PDEt},
Method->OptionValue[Method],
AccuracyGoal->OptionValue[AccuracyGoal],
PrecisionGoal->OptionValue[PrecisionGoal],
WorkingPrecision->OptionValue[WorkingPrecision],
MaxSteps->OptionValue[MaxSteps]
];
PDEsolution=Activate@PDEsolver;

tw+=AbsoluteTime[];
tm=PDEsolution["Domain"]//Last//Last;

KTsolution[<|"solution"->PDEsolution,"n"->n,"no"->no,"Δx"->Δx,"xi"->xi,"nt"->nt,"tw"->tw,"t0"->t0,"t1"->t1,"tm"->tm|>]
]

(** Output wrapper **)
ClearAll[KTsolution];
KTsolution[asoc_][key_]:=asoc[key] /; MemberQ[Keys[asoc], key]
KTsolution[asoc_][key_Symbol]:=asoc[ToString@key] /; MemberQ[Keys[asoc], ToString@key]
KTsolution[asoc_][t_]:=With[{sol=asoc["solution"],xi=asoc["xi"]},{xi,sol[t]}\[Transpose]/;!MissingQ[sol]]

MakeBoxes[KTsolution[asoc_Association],StandardForm]/;BoxFormUseIcons:=Module[{xi,n,no,Δx,nt},
(* [https://mathematica.stackexchange.com/a/99914] *)
(* [https://mathematica.stackexchange.com/a/79891] *)
xi=With[{a=asoc["xi"]},If[MissingQ[a],{},{a[[1]],a[[2]],Skeleton[Length[a]-3],Last@a}]];
n=With[{a=asoc["n"]},If[MissingQ[a],0,a]];
nt=With[{a=asoc["nt"]},If[MissingQ[a],0,a]];
Δx=With[{a=asoc["Δx"]},If[MissingQ[a],Missing,a]];
BoxFormArrangeSummaryBox[KTsolution,
KTsolution[asoc],
None,
{{"n = "<>ToString[n],"Δx = "<>ToString[Δx]},{Row[{"Subscript[x, i] = ",xi}],SpanFromLeft},{"nt="<>ToString[nt],SpanFromLeft}},
{},
StandardForm,
"Interpretable"->True
]
]


The only interesting part of the code above is the one for KTstepper as it implements the semi-discrete method-of-lines (MOL) finite volume discretion of Kurganov and Tadmore. KTsolver generates and solves the MOL ODE system using NDSolve and KTsolution is just an output wrapper.

In the above notation the Burger's equation is given by $$F(t,x,u)=F(u)=u^2$$ and $$Q(t,x,u,\partial_x u)=Q(\partial_x u)=\nu \partial_x u$$ in case of non-vanishing diffusion $$\nu>0$$. Apart from initial and boundary (in this example periodic ones) conditions the scheme requires the Eigenvalues of the Jacobian $$(\partial F/\partial u)$$, which in case of the Burger's equation is simply $$\partial F/\partial u = u$$. The following code computes three solutions at variable $$\nu$$:

burgersEq=Sequence[{t,x,u}|->0.5u^2,{t,x,u}|->u,{t,x,u,dudx,nu}|->nu*dudx,BCperiodic];

νi={0,0.01,0.05};
νiSols=Table[KTsolver[AccuracyGoal->6,PrecisionGoal->6][0.,2π,512][burgersEq,ν][x|->1/3+2 Sin[x]/3][0,2π],{ν,νi}]



The method works well even on small spatial grids:

ni={64,128,256,512};
niSols=Table[KTsolver[AccuracyGoal->6,PrecisionGoal->6][0.,2π,n][burgersEq,0][x|->1/3+2 Sin[x]/3][0,2π],{n,ni}]



Note that the solutions are completely free of spurious oscillations around the shock wave. The reason for this is the fact that the KT FV scheme is shock capturing and second order accurate in $$\Delta x$$. It employs a MUSCL reconstruction. In my experience with the numerical solution of this kind of conservation laws the KT FV is a very good compromise between simplicity and performance. It works shockingly well and can be implemented in only a few lines of code (most of the code above is just setup, NDSolve options and the output wrapper) using mainly list manipulations. The performance in Mathematica (especially when using the CompiledFunction) is also not bad. If you want to see the method flex its "MUSCLs" have a look into the examples studied in the paper or on wikipedia.

• "It works shockingly well" - I see what you did there May 10 at 6:44
• It's rare to see someone implementing PDE solver with certain generality, +11 :) . May 15 at 12:11
• Do you manage to reproduce the EXAMPLE 11 (Two-Dimensional Burgers-Type Equations) and EXAMPLE 12 (Two-Dimensional Buckley–Leverett Equation) in the paper? I implemented the scheme myself, so far it seems to work well for 1D PDE and 1D PDE system. But for 2D problem (EXAMPLE 11 and 12), I have to choose a coarse grid (about 25 for each dimension), otherwise NDSolve fails half way or the solution becomes unstable. Not sure if it's because I've made mistake(s) in implementation. May 19 at 2:10
• @xzczd I have not spent much time with the 2D generalization but I found some old code and I am able to run EXAMPLE 11 and 12. Have you tested the examples without diffusion? The diffusion term in EXAMPLE 11 is numerically very challenging (in terms of time-stepping) while the canonical diffusion term of EXAMPLE 12 is unproblematic. EXAMPLE 11 without diffusion takes 2s to run on 60x60 with Precision and AccuracyGoal of 3. EXAMPLE 11 with diffusion takes 1s to run on 25x25 with PG=AP=3. EXAMPLE 12 with diffusion takes 60s to run on 100x100 with PG=AG=3.
– N0va
May 20 at 7:48
• Yeah, I tested without diffusion, the phenomenon is similar (EXAMPLE 11 becomes stiff and stops halfway, EXAMPLE 12 becomes unstable). I've also adjusted PrecisionGoal and AccuracyGoal, but they don't help, either. May 20 at 9:10

I am not sure if there are options to improve NDSolve. I tried stiffness method, but that did not help.

But since the inviscid Burgers' equation has analytical solution, then you can use that.

Where $$f(x)$$ above is your initial condition. Mathematica DSolve gives this also:

DSolve[{D[u[x, t], t] + u[x, t]* D[u[x, t], x] == 0,
u[x, 0] == 1/3 + 2 Sin[x]/3}, u[x, t], x, t]


gives

Solve[u[x, t] == 1/3 (1 + 2 Sin[x - t u[x, t]]), u[x, t]]


The problem is that the solution is implicit. The solution $$u$$ shows on both sides.

But you could generate sample points as follows. Here is at t=1.5

Clear["Global*"]
epsilon := 0.01;
f[z_] := 1/3 + 2 Sin[z]/3;
u[x_, t_ /; t > 0] := u[x, t] = f[x - u[x, t - epsilon]*t]
u[x_, t_] := u[x, 0] = f[x];
maxTime = 1.5;
data = Table[u[x, t], {t, 0, maxTime, epsilon}];
Plot[Evaluate[Last@data], {x, 0, 2*Pi}, PlotStyle -> Red,  GridLines -> Automatic, GridLinesStyle -> LightGray]


And at t=2

Since the trajectory of that characteristic is (from same link above, Wikipedia)

you can just use implicitplot? May be better and more accurate:

t = 1.5;
f[x_] := 1/3 + 2 Sin[x]/3;
ParametricPlot[{t*f[x] + x, f[x]}, {x, 0, 2*Pi}, AspectRatio -> 1,
GridLines -> Automatic, GridLinesStyle -> LightGray,
PlotStyle -> Red]


t = 2;


t = 3;

t = 2*Pi;

• Thanks for your answer! However I still have 2 questions🤔: 1、The new code still cannot work well for $t=2\pi$, maybe when $t$ is too large we just cannot use NDSolve? 2、(Most importantly) Your code results works well for $t=2$, but mine doesn't. Which difference (might) led to the improvement of the result? Can more differences to be made such that $t=2\pi$ actually CAN be solved by NDSolve? May 6 at 9:36
• @Ho-Oh Use the second approach (the implicit plot) that works for t=2*Pi. The first one was just a crude approach to solve the implicit equation using some epsilon value which I set to 0.01. I am sure there is a better way to numerically iterate over it than what I showed quickly. May 6 at 9:38
• @Ho-Oh CAN be solved by NDSolve? Possibly. I am not an expert in all the options NDSolve takes. Hopefully someone will know better on this. I tried few options I know about, and they did not help. May 6 at 10:14
• This answer is incomplete as it fails to mention that the presented implicit analytical solution is only valid for a limited amount of time. The multivalued solutions presented here for $t>1.5$ are not physical solutions to the Burger's equation. I discussed this in my answer. I am also not able to give a solution using NDSolve but I discuss a workaround of sorts.
– N0va
May 6 at 10:32
• the presented implicit analytical solution is only valid for a limited amount of time Sure. Until characteristics don't intersect. Ok. Then OP could use the solution only for the correct time where this solution is valid. May 6 at 10:39