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 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]), 
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"}];\)\)

\(\(MultipleListPlot[UWData, 
        ExactData, 
PlotRange -> {{ν\ NN\  + 
                5, ν\ NN + 26}, 
            All}, 
PlotJoined -> 
          True, 
PlotLegend -> {"UW", "Exact"}];\)\)

\(\(MultipleListPlot[LWData, 
        ExactData, 
PlotRange -> {{ν\ NN\  + 
                5, ν\ NN + 26}, 
            All}, 
PlotJoined -> 
          True, 
PlotLegend -> {"LW", "Exact"}];\)\)

\(\(MultipleListPlot[FWData, 
        ExactData, 
PlotRange -> {{ν\ NN\  + 
                5, ν\ NN + 26}, 
            All}, 
PlotJoined -> 
          True, 
PlotLegend -> {"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 = 10;
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

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"}];\)\)

\(\(MultipleListPlot[UWData, 
        ExactData, 
PlotRange -> {{ν\ NN\  + 
                5, ν\ NN + 26}, 
            All}, 
PlotJoined -> 
          True, 
PlotLegend -> {"UW", "Exact"}];\)\)

\(\(MultipleListPlot[LWData, 
        ExactData, 
PlotRange -> {{ν\ NN\  + 
                5, ν\ NN + 26}, 
            All}, 
PlotJoined -> 
          True, 
PlotLegend -> {"LW", "Exact"}];\)\)

\(\(MultipleListPlot[FWData, 
        ExactData, 
PlotRange -> {{ν\ NN\  + 
                5, ν\ NN + 26}, 
            All}, 
PlotJoined -> 
          True, 
PlotLegend -> {"FW", "Exact"}];\)\)

\(\(MultipleListPlot[
        LWData - 
          ExactData, 
PlotRange -> {{ν\ NN\  + 
                5, ν\ NN + 26}, 
            All}, 
PlotJoined -> 
          True, 
PlotLegend -> {"LW", "Exact"}];\)\)

\(\(MultipleListPlot[
        FWData - 
          ExactData, 
PlotRange -> {{ν\ NN\  + 
                5, ν\ NN + 26}, 
            All}, 
PlotJoined -> 
          True, 
PlotLegend -> {"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 = 6000;
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"}];\)\)

Fourier Analysis

Clear[ν, FourierAnalysis, Scheme]
FourierAnalysis[ Scheme_, ν_] :=  
      Module[
{n, DFM, s, 
          r, λ}, 
n = 
          145; 
DFM[s_] = 
          Table[(E^(-2 π I (s - 1) (r - 1)/n)), {r, 1, n}]; 
λ[
            s_] = Scheme[
                DFM[s], ν][[1]]/DFM[
                s][[1]]; 
Show[
GraphicsArray[{
ListPlot[
Table[{π (s - 1)/n, 
                    Abs[λ[s]]}, {s, 1, 
                    Ceiling[
                      n/2]}], 
PlotRange -> {0.8, 
                    1. }, 
DisplayFunction -> 
                  Identity, 
Ticks -> {
Table[ θ, {θ, 
                        0, π/2, π/
                          8}], 
Automatic
}
], 
ListPlot[
Table[{π (s - 
                            1)/n, ((Arg[λ[s]] )/(2 ν  π (s - 1)/n))}, {s, 2, 
                    Ceiling[n/2]}], 
PlotRange -> 
                  Automatic, 
PlotStyle -> 
                  Hue[Random[]], 
DisplayFunction -> 
                  Identity, 
Ticks -> {
Table[ θ, {θ, 
                        0, π/2, π/
                          8}], 
Automatic
}
]
}], 
DisplayFunction -> $DisplayFunction
]
]

ν = 0.1;
Show[
FourierAnalysis[
          UWScheme, ν][[1, 
        2]], 
FourierAnalysis[
          LWScheme, ν][[1, 
        2]], 
FourierAnalysis[
          FWScheme, ν][[1, 
        2]], 
DisplayFunction -> $DisplayFunction
]