Needs["Graphics`MultipleListPlot`"]
Needs["Graphics`Legend`"]

Lax-Wendroff and Upwind 
(Subscript[u, t]) + (Subscript[(f(u)), x]) = 0

An alternate derivation of the LW scheme for the conservation form goes as follows
	u(x, t + Δt) = 
        u(x, t) + (Subscript[u, t])(x, t) Δt  + (Subscript[u, tt])(x, 
              t) Δ((t^2)/2) + …
Since the equation says
	(Subscript[u, t]) = -(Subscript[(f(u)), x])
we can differentiate with respect to t to get
	\(TraditionalForm\`u\_tt = \(\(-\(f(u)\)\_xt\) = \(\(-\((f' \((u)\) u\_t)\)\_x\) = \((f' \((u)\) \(f(u)\)\_x)\)\_x\)\)\)
we can substitute to get
	\(TraditionalForm\`u(x, t + Δt) = 
        u(x, t) - \(\(f(
                  u(x, t))\)\_x\) Δt\  + \(\((f' \((u)\) \(f(u)\)\_x)\)\_x\) \(Δt\^2\/2\) …\)
Dropping the Δ(t^3) and higher terms gives the approximation 
	\(TraditionalForm\`u(x, t + Δt) ≈ 
        u(x, t) - \(\(f(
                  u(x, t))\)\_x\) Δt\  + \(\((f' \((u)\) \(f(u)\)\_x)\)\_x\) Δt\^2\/2\)
which gives a numerical scheme as soon as we pick an approximation for the x derivatives.

We get the Lax-Wendroff scheme if we choose the centered difference of length 2Δx for the first derivative in the Δt term and centered differences of length Δx for the two derivatives in the Δ(t^2) term.  This gives
	\(TraditionalForm\`𝒰\_i\%\(n + 1\) = 𝒰\_i\%n - \ Δt \( f(𝒰\_\(i + 1\)\%n) - f(𝒰\_\(i - 1\)\%n)\)\/\(2\ Δx\)\)
\(TraditionalForm\`	\(+\ \(Δt\^2\/\(\(2\)\(\ \ \)\)\)\) \(\(\((f' \((𝒰\_\(i + 1/2\)\%n)\) \(f(𝒰\_\(i + 1\)\%n) - f(𝒰\_i\%n)\)\/Δx - 
                  f' \((𝒰\_\(i - 1/2\)\%n)\) \(f(𝒰\_i\%n) - f(𝒰\_\(i - 1\)\%n)\)\/Δx)\)\/Δx\)\(.\)\)\)
The only remaining problem is the need to evaluate f' at the "±1/2" points:  this is customarily solved by using the approximations 
	\(TraditionalForm\`f' \((𝒰\_\(i + 1/2\)\%n)\) ≈ \ \(f' \((𝒰\_\(i + 1\)\%n)\) + f(𝒰\_i\%n)\)\/2\)
	\(TraditionalForm\`f' \((𝒰\_\(i - 1/2\)\%n)\) ≈ \ \(f' \((𝒰\_\(i - 1\)\%n)\) + f(𝒰\_i\%n)\)\/2\)
It is simple to check that this reduces to the earlier LW formula when f(u) = a u.

If we choose the centered difference of length 2Δx for the first derivative in the Δt term and centered differences of length 2Δx for the two derivatives in the Δ(t^2) term we get a different scheme
	\(TraditionalForm\`𝒰\_i\%\(n + 1\) = 𝒰\_i\%n - \ Δt \( f(𝒰\_\(i + 1\)\%n) - f(𝒰\_\(i - 1\)\%n)\)\/\(2\ Δx\)\)
\(TraditionalForm\`	\(+\ \(Δt\^2\/\(\(2\)\(\ \ \)\)\)\) \(\(\((f' \((𝒰\_\(i + 1\)\%n)\) \(f(𝒰\_\(i + 2\)\%n) - f(𝒰\_i\%n)\)\/\(2  Δx\) - 
                  f' \((𝒰\_\(i - 1\)\%n)\) \(f(𝒰\_i\%n) - f(𝒰\_\(i - 2\)\%n)\)\/\(2  Δx\))\)\/\(2  Δx\)\)\(.\)\)\)
HW A:  Perform a linear linear stability analysis on this scheme.

We get an Upwind scheme if we stop after Δt term and use appropriate one sided derivatives to approximate the remaining first derivative. If we write 
	 \(TraditionalForm\`ℱ\_\(i ± 1/2\)\%n = \ \(f(𝒰\_\(i ± 1\)\%n)\  - f(𝒰\_i\%n)\)\/\(𝒰\_\(i ± 1\)\%n - 𝒰\_i\%n\)\) 
as approximations to the derivatives we get the upwind scheme
	(Subscript[𝒰, i]^(n + 1))=(Subscript[𝒰, i]^n)- Δ(t/Δ)x(({{(Subscript[ℱ, i + 1/2]^n)(
                        f((Subscript[𝒰, i + 1]^n)) - 
                          f((Subscript[𝒰, i]^n))), \(if\ ℱ\_\(i + 1/2\)\%n > 0\)}, {(Subscript[ℱ, i - 1/2]^n)(
                        f((Subscript[𝒰, i]^n)) - 
                          f((Subscript[𝒰, i - 1]^n))), \(if\ ℱ\_\(i - 1/2\)\%n > 0\)}})
This is not completely specified since it is clearly possible for both (or neither) of the conditions to hold.

Lets try them on Burger's Equation
	(Subscript[u, t]) + u (Subscript[u, x]) = 0
which is 
	(Subscript[u, t]) +  (Subscript[(f(u)), x]) = 0
where f(u) = (1/2) (u^2).
As an initial condition we will take a nice smooth pulse and see what happens.

Code for ConsLW and ConsUW

Here is the LW scheme.

ConsLWScheme[UIn_, f_, μ_]  :=  
      Module[
{fUIn}, 
fUIn = 
          f[UIn]; 
fDer  = ((( RotateLeft[fUIn] - fUIn))/(RotateLeft[UIn] - UIn)); 
fDerPlusHalf  =  
          1/2 (RotateLeft[fDer] + 
                fDer); 
fDerMinusHalf  =  
          1/2 (fDer + 
                RotateRight[
                  fDer]); 
UIn  - (1/2) μ (RotateLeft[fUIn] - 
                RotateRight[
                  fUIn]) + 
(1/2) (μ^2) (fDerPlusHalf (RotateLeft[fUIn] - fUIn) - 
                fDerMinusHalf (fUIn  -  
                      RotateRight[fUIn]))
]

Here is the UW scheme.

Clear[ConsUWScheme, UIn, f, μ]

    
ConsUWScheme[UIn_, f_, μ_]  :=  
      Module[
{fUIn, i}, 
fUIn = 
          f[UIn]; 
fDer  = ((( RotateLeft[fUIn] - fUIn))/(RotateLeft[UIn] - UIn)); 
 (* Print[
              fDer]; *) 
LDers  =  
          RotateLeft[fUIn] - fUIn; 
RDers  =  
          fUIn - RotateRight[fUIn]; 
UWDers  = 
          Table[
If[ 
              fDer[[i]] > 0,  
              fDer[[i]] 
                RDers[[
                  i]], -fDer[[
                    i]] LDers[[
                  i]]], 
{i, 
              Length[LDers]}]; 
UIn  -  μ UWDers
]

Here is the "other" scheme.

Clear[ConsOtherScheme, UIn, f, μ]

    
ConsOtherScheme[UIn_, f_, μ_]  :=  
      Module[
{fUIn, i}, 
fUIn = 
          f[UIn]; 
fDer  = ((( RotateLeft[fUIn] - fUIn))/(RotateLeft[UIn] - UIn)); 
UIn  - (1/2) μ (RotateLeft[fUIn] - 
                RotateRight[
                  fUIn]) + 
(1/8) (μ^2) (RotateLeft[fDer] (RotateLeft[fUIn, 2] - fUIn) - 
                RotateRight[
                    fDer] (fUIn  -  
                      RotateRight[fUIn, 
                        2]))
]

Testing the Code

Setting up some smooth initial conditions.

Δx = 0.091;
IC  =  
        Table[ (E^(-(x^2))), {x,  -6, 
            6, Δx}];
ListPlot[IC];

Running the LW, UW, and other code on the linear case as a check.

Δt = 0.05;
μ = Δt/Δx
f[u_] = u;
LWData  =  IC;
UWData  =  IC;
NoOfSteps = 200;
LWData  = 
        Nest[ ConsLWScheme[#, f, μ] &, IC, 
          NoOfSteps];
UWData  = 
        Nest[ ConsUWScheme[#, f, μ] &, IC, 
          NoOfSteps];
OtherData  = 
        Nest[ ConsOtherScheme[#, f, μ] &, IC, 
          NoOfSteps];
\(\(MultipleListPlot[\ LWData, UWData, IC, 
        OtherData +  .1, 
PlotLabel -> 
          ToString[
StringForm["Linear, n=``", 
              NoOfSteps]
], 
PlotJoined -> True, 
PlotRange -> {\(-0.1\), 
            1}, 
PlotLegend -> {"LW", "UW", "IC", "Other"}
];\)\)

Looks as though they work OK for the linear problem.  Lets check them on Burger's Eq.

Δt = 0.0125;
μ = Δt/Δx
f[u_] = 1/2 (u^2);
LWData  =  IC;
UWData  =  IC;
NoOfSteps = 200;
LWData  = 
        Nest[ ConsLWScheme[#, f, μ] &, IC, 
          NoOfSteps];
UWData  = 
        Nest[ ConsUWScheme[#, f, μ] &, IC, 
          NoOfSteps];
OtherData  = 
        Nest[ ConsOtherScheme[#, f, μ] &, IC, 
          NoOfSteps];
\(\(MultipleListPlot[\ LWData, UWData, IC, 
        OtherData, 
PlotLabel -> 
          ToString[
StringForm["Linear, n=``", 
              NoOfSteps]
], 
PlotJoined -> True, 
PlotRange -> 
          All, 
PlotRange -> {{80, 100}, 
            All}, 
Frame -> 
          True, 
PlotLegend -> {"LW", "UW", "IC", "Other"}
];\)\)

Δt = 0.0125;
μ = Δt/Δx
f[u_] = 1/2 (u^3);
LWData  =  IC;
UWData  =  IC;
OtherData  =  IC;
NoOfSteps = 200;
\(Table[
LWData\  = \ 
        ConsLWScheme[LWData, f, μ]; 
UWData\  = 
        ConsUWScheme[UWData, f, μ]; 
OtherData\  = 
        ConsOtherScheme[OtherData, 
          f, μ]; 
MultipleListPlot[\ LWData, UWData, 
        OtherData, 
PlotLabel -> 
          ToString[
StringForm["Linear, n=``", 
              n]
], 
PlotJoined -> 
          True, 
PlotRange -> {{80, 100}, 
            All}, 
Frame -> 
          True, 
PlotLegend -> {"LW", "UW", "Other"}
], 
{n, 1, 
        NoOfSteps}]\)