Fourier Comparison of
Lax-Wendroff etc.

We now have a collection of solvers for the simple 1-D wave equation
(1)	(∂_t) u + a (∂_x) u = 0
We have three schemes described below with ν  = 
        a  Δ(t/Δ)x.

The Upwind scheme
(UW) 	(Subscript[U, i]^(n + 1))={({{ν (Subscript[U, i + 1]^n) + (1 - ν) 
                        (Subscript[U, i]^n), if ν < 0}, {ν (Subscript[U, i - 1]^n) + (1 - ν) 
                        (Subscript[U, i]^n), if ν > 0}})	
with CFL and stability condition |ν|≤1.

The Lax-Wendroff scheme
(LW)	(Subscript[U, i]^(n + 1))=\(\(ν\/2\) \((1 + ν)\)\ U\_\(i - 1\)\%n + \((1 - ν\^2)\) U\_i\%n - \(ν\/2\) \((1 - ν)\)\)(Subscript[U, i + 1]^n)
with CFL and stability condition |ν|≤1.

The Fat-Wendroff scheme
(FW)	(Subscript[U, i]^(n + 1))=\(\(ν\/24\) \((ν\^2 - 1)\)\ \((2 + ν)\) 
                      U\_\(i - 2\)\%n + \(ν\/6\) \((4 - ν\^2)\)\ \((1 + ν)\) U\_\(i - 1\)\%n\)
\(\(+\(1\/4\)\) \((1 - ν\^2)\) \((4 - ν\^2)\) U\_i\%n\)
+ (ν/6) ((ν^2) - 
                            4) (1 - ν)(Subscript[U, i + 1]^n) +\(\(ν\/24\) \((1 - ν\^2)\)\ \((2 - ν)\)\ U\_\(i + 2\)\%n\)
with CFL and condition |ν|≤2.

The plan is to compare these three schemes.  First we are going to look at their behavior on our square pulse and describe the effects of the schemes qualitatively.  Then we are going to perform our Fourier analysis to compare and contrast the effect of the schemes on the Discrete Fourier Modes (DFM).

Qualitative Analysis — Square Pulse

Here are the three different solutions to the square pulse problem along with the exact solution.

Initializations

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

Initial Condition

SquarePulse[ζ_] =  
        If[ 1 ≤ Abs[ζ] ≤ 2,  1, 
          0];
Δx = 0.1;  X =  20;
SquareInit  = 
        Table[ SquarePulse[x],  {x, 0, 
            X, Δx}];
ListPlot[
        SquareInit, 
PlotRange -> All];

Update Schemes

Clear[UWScheme, UIn, ν]
UWScheme[UIn_, ν_] :=  
      Module[
{UOut}, 
If[ν < 
            0, 
UIn - ν (RotateLeft[UIn]  - 
                  UIn), 
UIn - ν (UIn - 
                  RotateRight[
                    UIn])
]
]

Clear[LWScheme, UIn, ν]
LWScheme[UIn_, ν_] :=  
      Module[
{UOut}, 
ν/
              2 (1 + ν) RotateRight[UIn]  +  (1 - (ν^2)) 
            UIn  - ν/2 (1 - ν) 
            RotateLeft[UIn]
]

Clear[FWScheme, UIn, ν]
FWScheme[UIn_, ν_] :=  
      Module[
{UOut}, 
(1/24) ((ν^2) - 1) ν (2 + ν) 
            RotateRight[UIn, 
              2]
 + (1/6) (4 - (ν^2)) ν (1 + ν) RotateRight[UIn, 
              1]
 + (1/4) ((ν^2) - 
                4) ((ν^2) - 1) 
            UIn
 + (1/6) (4 - (ν^2)) ν (ν - 1) RotateLeft[UIn, 
              1]
 + (1/24) ((ν^2) - 
                1) ν (ν - 2) 
            RotateLeft[UIn, 2]
]

Runs

Here are some runs with a = 0.9, Δt = 0.1, and Δx = 0.1giving ν = 0.9.

Δt = 0.1;
a = 0.9;
ν = 
        a Δ(t/Δ)x;
NN = 100;
UWData  =  
        Nest[ UWScheme[#, ν] &, SquareInit, 
          NN];
LWData  =  
        Nest[ LWScheme[#, ν] &, SquareInit, 
          NN];
FWData  =  
        Nest[ FWScheme[#, ν] &, SquareInit, 
          NN];
ExactData  = 
        Table[ SquarePulse[x - NN a Δt],  {x, 0, 
            X, Δx}];
\(\(MultipleListPlot[UWData, LWData, FWData, 
        ExactData, 
PlotRange -> {{ν\ NN\  + 
                5, ν\ NN + 26}, 
            All}, 
PlotJoined -> 
          True, 
PlotLegend -> {"UW", "LW", "FW", "Exact"}];\)\)

Which scheme is the best and why?  I guess we should experiment a bit with a longer run.

Δt = 0.1;
a = 0.9;
ν = 
        a Δ(t/Δ)x;
NN = 600;
UWData  =  
        Nest[ UWScheme[#, ν] &, SquareInit, 
          NN];
LWData  =  
        Nest[ LWScheme[#, ν] &, SquareInit, 
          NN];
FWData  =  
        Nest[ FWScheme[#, ν] &, SquareInit, 
          NN];
ExactData  = 
        Table[ 
SquarePulse[
            Mod[x - NN a Δt, X]],  
{x, 0, 
            X, Δx}];
\(\(MultipleListPlot[UWData, LWData, FWData, 
        ExactData, 
PlotRange -> {
Mod[{ν\ NN\  - 10, ν\ NN + 26}, 
              Length[UWData]], 
All}, 
PlotJoined -> 
          True, 
PlotLegend -> {"UW", "LW", "FW", "Exact"}];\)\)

I am surprised that the FW scheme looks the best! Lets try it with a different a.

Δt = 0.1;
a = 0.1;
ν = 
        a Δ(t/Δ)x;
NN = 600;
UWData  =  
        Nest[ UWScheme[#, ν] &, SquareInit, 
          NN];
LWData  =  
        Nest[ LWScheme[#, ν] &, SquareInit, 
          NN];
FWData  =  
        Nest[ FWScheme[#, ν] &, SquareInit, 
          NN];
ExactData  = 
        Table[ 
SquarePulse[
            Mod[x - NN a Δt, X]],  
{x, 0, 
            X, Δx}];
\(\(MultipleListPlot[UWData, LWData, FWData, 
        ExactData, 
PlotRange -> {
Mod[{ν\ NN\  - 10, ν\ NN + 26}, 
              Length[UWData]], 
All}, 
PlotJoined -> 
          True, 
PlotLegend -> {"UW", "LW", "FW", "Exact"}];\)\)

With 1< | ν | <2 we know that the LW etc. schemes are unstable lets see what happens to the FW scheme!  Looks as though the FW scheme is also unstable for this range!

Δt = 0.1;
a = 1.3;
ν = 
        a Δ(t/Δ)x;
NN = 5;
FWData  =  
        Nest[ FWScheme[#, ν] &, SquareInit, 
          NN];
ExactData  = 
        Table[ 
SquarePulse[
            Mod[x - NN a Δt, X]],  
{x, 0, 
            X, Δx}];
\(\(MultipleListPlot[FWData, 
        ExactData, 
PlotRange -> {{0, 50}, 
            All}, 
PlotJoined -> 
          True, 
PlotLegend -> {"FW", "Exact"}];\)\)

Qualitative Analysis — Smooth Pulse

Here are the three different solutions to a smooth pulse problem along with the exact solution.

Initializations

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

Initial Condition

Plot[ 
      (E^(-(1/((((x - 1))^2) (((x - 2))^2))))), {x, -2, 
        4}, 
PlotRange -> All]

SmoothPulse[ζ_] =  
        (E^(-2 (((ζ - 1.5))^2)));
Δx = 0.1;  X =  20;
SmoothInit  = 
        Table[ SmoothPulse[x],  {x, 0, 
            X, Δx}];
ListPlot[
        SmoothInit, 
PlotRange -> All];

Runs

Here are some runs with a = 0.9, Δt = 0.1, and Δx = 0.1giving ν = 0.9.

Δt = 0.1;
a = 0.9;
ν = 
        a Δ(t/Δ)x;
NN = 100;
UWData  =  
        Nest[ UWScheme[#, ν] &, SmoothInit, 
          NN];
LWData  =  
        Nest[ LWScheme[#, ν] &, SmoothInit, 
          NN];
FWData  =  
        Nest[ FWScheme[#, ν] &, SmoothInit, 
          NN];
ExactData  = 
        Table[ SmoothPulse[x - NN a Δt],  {x, 0, 
            X, Δx}];
\(\(MultipleListPlot[UWData, LWData, FWData, 
        ExactData, 
PlotRange -> {{ν\ NN\  + 
                5, ν\ NN + 26}, 
            All}, 
PlotJoined -> 
          True, 
PlotLegend -> {"UW", "LW", "FW", "Exact"}];\)\)

Which scheme is the best and why?  I guess we should experiment a bit with a longer run.

Δt = 0.1;
a = 0.9;
ν = 
        a Δ(t/Δ)x;
NN = 600;
UWData  =  
        Nest[ UWScheme[#, ν] &, SmoothInit, 
          NN];
LWData  =  
        Nest[ LWScheme[#, ν] &, SmoothInit, 
          NN];
FWData  =  
        Nest[ FWScheme[#, ν] &, SmoothInit, 
          NN];
ExactData  = 
        Table[ 
SmoothPulse[
            Mod[x - NN a Δt, X]],  
{x, 0, 
            X, Δx}];
\(\(MultipleListPlot[UWData, LWData, FWData, 
        ExactData, 
PlotRange -> {
Mod[{ν\ NN\  - 10, ν\ NN + 26}, 
              Length[UWData]], 
All}, 
PlotJoined -> 
          True, 
PlotLegend -> {"UW", "LW", "FW", "Exact"}];\)\)

I am surprised that the FW scheme looks the best! Lets try it with a different a.

Δt = 0.1;
a = 0.1;
ν = 
        a Δ(t/Δ)x;
NN = 600;
UWData  =  
        Nest[ UWScheme[#, ν] &, SmoothInit, 
          NN];
LWData  =  
        Nest[ LWScheme[#, ν] &, SmoothInit, 
          NN];
FWData  =  
        Nest[ FWScheme[#, ν] &, SmoothInit, 
          NN];
ExactData  = 
        Table[ 
SmoothPulse[
            Mod[x - NN a Δt, X]],  
{x, 0, 
            X, Δx}];
\(\(MultipleListPlot[UWData, LWData, FWData, 
        ExactData, 
PlotRange -> {
Mod[{ν\ NN\  + 5, ν\ NN + 25}, 
              Length[UWData]], 
All}, 
PlotJoined -> 
          True, 
PlotLegend -> {"UW", "LW", "FW", "Exact"}];\)\)

Second Attempts

pic

A semi-obvious generalization of this is to use a quadratic interpolation between the three points 
	\(TraditionalForm\`U\_\(i - 1\)\%n, \ \(U\_\(i + 1\)\%n &\)\ U\_\(i + 1\)\%n\)
and move this up the characteristic.  The benefit of this is that it can deal with both signs of ν with one formula.  This idea gives the Lax-Wendroff scheme.
(7)	(Subscript[U, i]^(n + 1))=\(\(ν\/2\) \((1 + ν)\)\ U\_\(i - 1\)\%n + \((1 - ν\^2)\) U\_i\%n - \(ν\/2\) \((1 - ν)\)\)(Subscript[U, i + 1]^n)
which is very easily implemented.

Initial Conditions

SquarePulse[ζ_] =  
        If[ Abs[ζ] < 1,  1, 0];
Δx = 0.1;  X =  10;
SquareInit  = 
        Table[ SquarePulse[x],  {x, -X, 
            X - Δx, Δx}];
ListPlot[
        SquareInit, 
PlotRange -> All];

GaussPulse[ζ_] =  
        (E^(-(ζ^2)));
Δx = 0.1;  X =  10;
GaussInit  = 
        Table[ GaussPulse[x],  {x, -X, 
            X - Δx, Δx}];
ListPlot[
        GaussInit, 
PlotRange -> All];

Update Schemes LW and UW

The plan is to implement the update scheme (7).  For the moment we restrict attention to constant coefficients on a periodic domain.

Clear[SchemeLW, UIn, ν]
SchemeLW[UIn_, ν_] :=  
      Module[
{UOut}, 
ν/
              2 (1 + ν) RotateRight[UIn]  +  (1 - (ν^2)) 
            UIn  - ν/2 (1 - ν) 
            RotateLeft[UIn]
]

Here is a test!

a  =  1;
Δt = 0.04;
ν  =  
      a Δt/Δx
NewData  =  SquareInit;
Data  =  
        Table[
NewData  =  
            SchemeLW[
              NewData, ν], 
{400}];

    
ListPlot3D[ 
        Data, 
Mesh -> 
          False, 
PlotRange -> {0, 1}];

Lets look at the pictures as a function of time.

Table[
ListPlot[
Data[[
            i]], 
PlotRange -> {-0.3, 1.3}, 
PlotLabel -> 
            i, 
PlotJoined -> 
            True
], 
{i,  
          Length[Data]}];

It looks as though the square pulse grows some wiggles bt it stays roughly the correct height and width!  Of course, it will work better on a smoother intial profile.

a  =  1;
Δt = 0.04;
ν  =  
      a Δt/Δx
NewData  =  GaussInit;
Data  =  
        Table[
NewData  =  
            SchemeLW[
              NewData, ν], 
{400}];

    
ListPlot3D[ 
        Data, 
Mesh -> 
          False, 
PlotRange -> {0, 1}];

Table[
ListPlot[
Data[[
            i]], 
PlotRange -> {-0.3, 1.3}, 
PlotLabel -> 
            i, 
PlotJoined -> 
            True
], 
{i,  
          Length[Data]}];

We still need to obey the CFL condition.  Here is a short run violating it!

a  =  1;
Δt = 0.12;
ν  =  
      a Δt/Δx
NewData  =  GaussInit;
Data  =  
        Table[
NewData  =  
            SchemeLW[
              NewData, ν], 
{100}];

ListPlot3D[ 
        Data, 
Mesh -> 
          False, 
PlotRange -> {0, 1}];

It is pretty dramatic.

Table[
ListPlot[
Data[[
            i]], 
PlotRange -> {-0.3, 1.3}, 
PlotLabel -> 
            i, 
PlotJoined -> 
            True
], 
{i,  
          Length[Data]}];

Fourier Analysis

The Fourier Stability Analysis lets us understand this behavior. The idea is simple in the same way that any pulse f(x) on the region  0 < x < 2 π can be reprented by the Fourier series 
	f(
          x) ∼ (∑_(j = 0)^∞) 
                (Subscript[c, j]) (e^(i j x)) where \(TraditionalForm\`c\_j = ∫\_0\%\(2  π\)\(f(
              x)\) \(e\^\(\(-i\)\ j\ x\)\) dx/π\) and (Subscript[c, 0]) = (∫_0^(2  π))f(
              x) dx/(2  π).
Just in case you do not believe me here is a function and the sum of the first few terms of the Fourier Series (FS).  Since the function is not very "wiggly" the FS does not need all that many terms to match the function well.

Clear[c, j]
f[x_]  = 
      x (E^(Sin[4 x] + Sin[(x^2)/π]))
c[0] :=  c[0] = 
        NIntegrate[ 
            f[x] , {x, 0, 
              2 π}]/(2  π)
c[j_] := c[j] =  
        NIntegrate[ 
            f[x] (E^(-I j x)), {x, 0, 
              2 π}]/π
FS[n_]  := 
      Re[Sum[ c[j] (E^(I j x)), {j, 0, n}]]
Plot[{f[x], FS[12]}, {x, 0, 
          2 π}, 
PlotStyle -> {{Thickness[0.02], 
              Hue[0]}, {Thickness[0.01], 
              Hue[0.3]}}, 
PlotRange -> 
          All
];

Exactly the same thing happens for a discrete data set. The only difference is that all the integrals become sums and the normalizations (the 1/π) turns up slightly differently! The coefficients are
	(Subscript[c, s])=(1/(Sqrt[n]))(∑_(r = 1)^n)(Subscript[u, r]) 
                (e^(2  π  i(r - 1) (s - 1)/n))	for 1 ≤ s ≤ n.
The input data (Overscript[u, ⇀]) = {(Subscript[u, 1]), (Subscript[u, 2]), (Subscript[u, 3]), …, 
          (Subscript[u, n])} is reconstructed by
	(Subscript[u, r])=(1/(Sqrt[n]))(∑_(s = 1)^n)(Subscript[c, s]) 
                (e^(-2 π  i(r - 1) (s - 1)/n))	for 1 ≤ r ≤ n.
Just in case you do not believe me here is a check

Clear[c, j]
f[x_]  =  
        If[x > 1, (E^(Sin[x] + Sin[(x^2)/π])), 
          0];
u  =  
        Table[ f[x], {x, 0,  2  π, 0.1}];
ListPlot[u]
n  =  Length[u]
c  =  
        Table[
(1/(Sqrt[n])) (∑_(r = 1)^n) 
                  u[[
                  r]] (E^(2 π I (r - 1) (s - 1)/n)), 
{s, 1, 
            n}];
ListPlot[
        Abs[c], 
PlotRange -> All];

There is a slight difference here between the FS and this Discrete Fourier Transform (DFT).  IN the FS the modes with the higher indices have the most "wiggles" and there is no "wiggliest" mode.  In the DFT it is the modes with the indices in the middle of the range that have more "wiggles" and the central mode is the wiggliest:  if there is no central mode then the two on either side are tied for "wiggliest". The wiggliest mode has n/2 wiggles and a ± pattern.

We can reconstruct the original data set (Overscript[u, ⇀]) using the formula to get

uRecon  =  
        Table[
(1/(Sqrt[n])) (∑_(s = 1)^n) 
                  c[[
                  s]] (E^(-2 π I (r - 1) (s - 1)/n)), 
{r, 1, 
            n}];
PP = 28
PartialuRecon  =  
        Table[
(1/(Sqrt[n])) (((∑_(s = 1)^(PP + 1))c[[
                        s]] (E^(-2 π I (r - 1) (s - 1)/n))) + ((∑_(s = n - PP)^n) c[[
                        s]] (E^(-2 π I (r - 1) (s - 1)/n)))), 
{r, 1, n}];
ListPlot[Re[uRecon]];
ListPlot[Re[PartialuRecon]]

The upshot of this is that we can understand the effect of a numerical scheme for a constant-coefficient, linear equation, with periodic boundary conditions by looking at how it behaves on these Discrete Fourier Modes.  We can do it analytically just the way the book does but we can also do it on the computer.

First we need to make a DFM data set.  Here it is!

Remove[MultipleListPlot]
Needs["Graphics`MultipleListPlot`"]

DFM[n_, s_]  :=  
        Table[(E^(-2 π I (s - 1) (r - 1)/n)), {r, 1, n}];
nn = 24;
ss = 2;
Table[
α = 
        E^(I τ); 
Test = α (DFM[nn, 
                ss] + DFM[
                nn, -ss]); 
MultipleListPlot[ 
        Re[Test], Im[Test], 
        Abs[Test], 
PlotJoined -> 
          True], 
{τ, 0, 2  π, 0.1}]

Next we need to throw it into one of our update schemes.

Clear[ν]
Simplify[SchemeLW[DFM[nn, ss], ν]]
Simplify[SchemeLW[DFM[nn, ss], ν]]/DFM[nn, ss]

Clear[ν, nn, ss]
nn = 24;
ss = 2;
ν = 0.7
SchemeLW[DFM[nn, ss], ν]/DFM[nn, ss]

Here is a check on the Amplification factors for all of the modes for the LW scheme.

Clear[ν, nn, ss]
nn = 124;
ν = -0.9
Table[
{ss, α = 
          Simplify[SchemeLW[
                  DFM[nn, 
                    ss], ν][[1]]/DFM[nn, ss][[1]]], Abs[α], 1, 
        Arg[α]}, 
{ss, 1, nn}]

Here is the same analysis for Scheme2

?Scheme2

Clear[ν, nn, ss]
nn = 124;
ν = -0.9
Stab2 =Table[
{ss, α = 
            Simplify[Scheme2[
                    DFM[nn, 
                      ss], ν][[1]]/DFM[nn, ss][[1]]], Abs[α], 
          Arg[α]}, 
{ss, 1, 
          nn}]

How do we compare the two schemes.