General Routines

\(discretize::usage = 
	
      "Returns a discretization of the interval [a,b] into nCells cells.  The number of nodes returned is nCells+1."; 


    discretize[a_, b_, nCells_] := 
      Module[{Δx}, 
		
        Δx = \((b - a)\)\/nCells; 
		
        Return[Table[N[a + Δx\ i], {i, 0, nCells}]]; 
		]; 
    \)

\(initialCondition::usage = 
      "Returns a set of points {xi,yi} that represent the discretized initial conditions of a 1D finite difference problem."; 


    initialCondition[a_, b_, nCells_, iCFunction_] := 
      Module[{xIC, yIC}, 
		xIC = discretize[a, b, nCells]; 
		
        yIC = Map[iCFunction, xIC]; 
		Return[Thread[{xIC, yIC}]]; 

        		]; \)

To obtain the initial conditions we start the scheme with the solution vector at t_0 and t_{-1} the same i.e. both icy.

\(\(advectionOneDimension[xL_, \ xR_, \ iCFunction_, \ nCells_, \ 
        scheme_, \ σ_, nTimeSteps_]\  := 
      Module[{ic, \ icx, \ icy}, 
		
        ic\  = \ initialCondition[xL, \ xR, \ nCells, \ iCFunction]; 
		
        icx\  = \ Map[First, \ ic]; 
		icy\  = \ Map[Last, \ ic]; 

        					
        yAtnTimeSteps = 
          FixedPoint[scheme[#, σ]&, {icy, icy}, nTimeSteps]; 
		\ 
        Return[Thread[{icx, \ yAtnTimeSteps[\([1]\)]}]]; 
		]; \)\)

Page 339

Problem 8.6-8.8

Exact solution:

\(\(uExact[a_, \ k_, \ t_]\  := 
\ \ \ \ 
      Plot[
		
        If[0 < x - a\ t < 1, 
		Sin[k\ π\ \((x\  - \ a\ t)\)], 
			
          0], \ {x, \ \(-0.5\), \ 2.5}, 
\ \ \ \ \ \ 
        PlotStyle\  -> \ RGBColor[0, \ 1, \ 0], 
\ \ \ \ \ \ 
        DisplayFunction\  -> \ Identity]; \)\)

\(plotTimes = Table[1\/10\ i, \ {i, \ 0, \ 10}]; \ \ 
     (*{0, \ 0.1, \ 1, \ 2}; *) \ 
\)\(temp\  = 
      Table[uExact[1, \ 1, \ plotTimes[\([i]\)]], 
        \ {i, \ 1, Length[plotTimes]}]; 

    Show[temp, \ DisplayFunction\  -> \ $DisplayFunction]; \)

Here we investigate how to get an second condition needed for the leapfrog method.  I tried to invert the difference scheme.  This seamed to work but I think that many terms are needed to have this transformation work.

\(ic = initialCondition[0, 1, 50, Sin[π\ #]&]; 
icy = Map[Last, ic]; 
    

f1[v_List, σ_] := 
      1\/\((1 + σ)\) - 
        \(σ\/\((1 + σ)\)\^2\) RotateRight[v] + 
        \(σ\^2\/\((1 + σ)\)\^3\) RotateRight[v, 2]; 

    f2[v_List, σ_] := 
      1 - σ + σ\^2 + 
        \((\(-σ\) + 2\ σ\^2)\) RotateRight[v] + 
        \(σ\^2\) RotateRight[v, 2]; 

lp0 = ListPlot[icy]; 

    lp1 = ListPlot[f1[icy, 0.5]]; 
lp2 = ListPlot[f2[icy, 0.5]]; 

    Show[{lp0, lp2}]; \)

Another option would be just to use the same vector for the other condition needed for the leap frog method.

Numerical schemes:

\(\(upwind[{v1_List, v0_List}, σ_] := 
      Module[{}, 
		
        \(Return[

          			{\((1 - σ)\)\ v1 + σ\ RotateRight[v1], v1}

          				]; \)
		]; \)\)

\(\(laxFriedrichs[{v1_List, v0_List}, \ σ_]\  := 
      Module[{}, 
		
        \(Return[

          			{\(1\/2\) 
              \((\((1\  + \ σ)\)\ RotateRight[v1]\  + \ 
                  \((1\  - \ σ)\)\ RotateLeft[v1])\), v1}
			]; 
        \)
		]; \)\)

\(\(laxWendroff[{v1_List, v0_List}, σ_] := 
      Module[{}, 
		
        \(Return[

          			{\(σ\/2\) \((1 + σ)\) RotateRight[v1] + 
              \((1 - σ\^2)\) v1 + 
              \(σ\/2\) \((\(-1\) + σ)\) RotateLeft[v1], v1}

          		]; \)
		]; \)\)

\(\(leapFrog[{v1_List, v0_List}, σ_] := 
      Module[{}, 
		
        \(Return[

          			{v0 - σ\ RotateLeft[v1] + σ\ RotateRight[v1], 
            v1}
		]; \)
		]; \)\)

\(\(Exercise[xL_, xR_, k_, scheme_, \ nCells_, \ σ_, \ nTimeSteps_]
         := \ Module[{a\  = \ 1, \ 
          Δx\  = \ \((xR\  - \ xL)\)\/nCells, \ ue, \ au, \ 
          lpa}, \ 
ue\  = \ 
          uExact[a, \ k, σ \( Δx\/a\) nTimeSteps]; \ 
        
ua\  = \ 
          advectionOneDimension[xL, \ xR, 
            Function[x, If[0 < x < 1, Sin[k\ π\ x], 0]], \ nCells, \ 
            scheme, \ σ, \ nTimeSteps]; \ 
        
lpa\  = \ 
          ListPlot[ua, DisplayFunction\  -> \ Identity]; 
		
		
        Show[{ue, \ lpa}, \ 
\ \ \ \ \ \ \ \ \ \ 
          DisplayFunction\  -> \ $DisplayFunction
			]; 
		]; \)\)

Exercise[-.5, 2.5, 1, laxFriedrichs, 150, 0.25, 10]Exercise[-.5, 2.5, 10, laxFriedrichs, 150, 0.25, 10]

Exercise[-.5, 2.5, 1, upwind, 150, 0.25, 10]Exercise[-.5, 2.5, 10, upwind, 150, 0.25, 10]

Exercise[-.5, 2.5, 1, leapFrog, 150, 0.25, 10]Exercise[-.5, 2.5, 10, leapFrog, 150, 0.25, 10]

Exercise[-.5, 2.5, 1, laxWendroff, 150, 0.25, 10]Exercise[-.5, 2.5, 10, laxWendroff, 150, 0.25, 10]

Page 340

Problem 8.12

\(ul = 2; 
ur = 0; 
sl = 1; 
u0[x_] := If[x < sl, ul, ur]; \)

\(\(uExact[a_, t_]\  := 
\ \ \ \ 
      Plot[u0[x - a\ t], \ {x, \ 2.0, \ 4.0}, 
		
        PlotStyle\  -> \ RGBColor[0, \ 1, \ 0], 
		
        DisplayFunction\  -> \ Identity]; \)\)

Numerical schemes:

To implement the boundary conditions we add an extra element to the ends of the array, compute with this element and then remove the two end elements.

\(\(upwind[{v1_List, v0_List}, σ_] := 
      Module[{vPadded, vNextPadded, vNext}, 
		
        vPadded = Join[{ul}, v1, {ur}]; 
		
        vNextPadded = 
          \((1 - σ)\)\ vPadded + σ\ RotateRight[vPadded]; 
		
        vNext = Take[vNextPadded, {2, Length[vNextPadded] - 1}]; 
		
        Return[{vNext, v1}]; 
		]; \)\)

\(\(laxFriedrichs[{v1_List, v0_List}, \ σ_]\  := 
      Module[{v1Padded, vNextPadded, vNext}, 
		
        v1Padded = Join[{ul}, v1, {ur}]; 
		
        vNextPadded = 
          \(1\/2\) 
            \((\((1\  + \ σ)\)\ RotateRight[v1Padded]\  + \ 
                \((1\  - \ σ)\)\ RotateLeft[v1Padded])\); 
		
        vNext = Take[vNextPadded, {2, Length[vNextPadded] - 1}]; 
		
        Return[{vNext, v1}]; 
		]; \)\)

\(\(laxWendroff[{v1_List, v0_List}, σ_] := 
      Module[{v1Padded, vNextPadded, vNext}, 
		
        v1Padded = Join[{ul}, v1, {ur}]; 
		
        vNextPadded = 
          \(σ\/2\) \((1 + σ)\) RotateRight[v1Padded] + 
            \((1 - σ\^2)\) v1Padded + 
            \(σ\/2\) \((\(-1\) + σ)\) RotateLeft[v1Padded]; 

        		vNext = Take[vNextPadded, {2, Length[vNextPadded] - 1}]; 
		
        Return[{vNext, v1}]; 
		]; \)\)

\(\(leapFrog[{v1_List, v0_List}, σ_] := 
      Module[{v0Padded, v1Padded, vNextPadded, v1Next, v0Next}, 
		
        v0Padded = Join[{ul}, v0, {ur}]; 
		
        v1Padded = Join[{ul}, v1, {ur}]; 
		
        vNextPadded = 
          v0Padded - σ\ RotateLeft[v1Padded] + 
            σ\ RotateRight[v1Padded]; 
		
        v1Next = Take[vNextPadded, {2, Length[vNextPadded] - 1}]; 
		
        v0Next = Take[v1Padded, {2, Length[v1Padded] - 1}]; 
		
        Return[{v1Next, v0Next}]; 
		]; \)\)

\(\(Exercise[xL_, xR_, scheme_, \ nCells_, \ σ_, \ nTimeSteps_] := 
      \ Module[{a\  = \ 1, \ 
          Δx\  = \ \((xR\  - \ xL)\)\/nCells, \ ue, \ au, \ 
          lpa}, \ 
ue\  = \ 
          uExact[a, σ \( Δx\/a\) nTimeSteps]; \ 
        
ua\  = \ 
          advectionOneDimension[xL, \ xR, Function[x, u0[x]], \ nCells, \ 
            scheme, \ σ, \ nTimeSteps]; \ 
        
lpa\  = \ 
          ListPlot[ua, DisplayFunction\  -> \ Identity]; 
		
		
        Show[{ue, lpa}, \ 
\ \ 			PlotRange -> {{2.0, 4.0}, All}, 

          \ \ \ \ \ \ \ \ \ \ DisplayFunction\  -> \ $DisplayFunction

          			]; 
		]; \)\)

Exercise[0, 4, upwind, 80, .8, 50]

Exercise[0, 4, laxFriedrichs, 80, .8, 50]

Exercise[0, 4, laxWendroff , 80, .8, 50]

Exercise[0, 4, leapFrog , 80, .8, 50]

Problem 8.13

\(pL = 1.5; 
pR = 3.5; \)

u0[x_] := If[0 < x < 1, Sin[2 π n x], 0];

\(n = 2; 
Plot[u0[x], {x, \(-2\), 2}]; 
n = 4; 

    Plot[u0[x], {x, \(-2\), 2}]; \)

\(\(uExact[a_, t_]\  := 
\ \ \ \ 
      Plot[u0[x - a\ t], \ {x, \ pL, \ pR}, 
		
        PlotStyle\  -> \ RGBColor[0, \ 1, \ 0], 
		
        DisplayFunction\  -> \ Identity]; \)\)

uExact[1, 2.1] // Show[#, DisplayFunction -> $DisplayFunction]&;

Numerical schemes:

\(\(upwind[{v1_List, v0_List}, σ_] := 
      Module[{}, 
		
        \(Return[{\((1 - σ)\)\ v1 + σ\ RotateRight[v1], v1}]; 
        \)
		]; \)\)

\(\(laxFriedrichs[{v1_List, v0_List}, \ σ_]\  := 
      Module[{}, 
		
        \(Return[{
            \(1\/2\) 
              \((\((1\  + \ σ)\)\ RotateRight[v1]\  + \ 
                  \((1\  - \ σ)\)\ RotateLeft[v1])\), v1}]; \)
		]; 
    \)\)

\(\(laxWendroff[{v1_List, v0_List}, σ_] := 
      Module[{}, 
		
        \(Return[{
            \(σ\/2\) \((1 + σ)\) RotateRight[v1] + 
              \((1 - σ\^2)\) v1 + 
              \(σ\/2\) \((\(-1\) + σ)\) RotateLeft[v1], v1}]; 
        \)
		]; \)\)

\(\(leapFrog[{v1_List, v0_List}, σ_] := 
      Module[{}, 
		
        \(Return[{v0 - σ\ RotateLeft[v1] + σ\ RotateRight[v1], 
            v1}]; \)
		]; \)\)

\(\(Exercise[xL_, xR_, scheme_, \ nCells_, \ σ_, \ nTimeSteps_] := 
      \ Module[{a\  = \ 1, \ 
          Δx\  = \ \((xR\  - \ xL)\)\/nCells, \ ue, \ au, \ 
          lpa}, \ 
ue\  = \ 
          uExact[a, σ \( Δx\/a\) nTimeSteps]; \ 
        
ua\  = \ 
          advectionOneDimension[xL, \ xR, Function[x, u0[x]], \ nCells, \ 
            scheme, \ σ, \ nTimeSteps]; \ 
        
lpa\  = \ 
          ListPlot[ua, DisplayFunction\  -> \ Identity]; 
		
		
        Show[{ue, lpa}, \ 
\ \ 			PlotRange -> {{pL, pR}, All}, 

          \ \ \ \ \ \ \ \ \ \ DisplayFunction\  -> \ $DisplayFunction

          			]; 
		]; \)\)

\(n = 2; 
Exercise[0, 3.5, upwind, 70, .8, 50]\)Exercise[0, 3.5, laxFriedrichs, 70, .8, 50]Exercise[0, 3.5, laxWendroff , 70, .8, 50]Exercise[0, 3.5, leapFrog , 80, .8, 50]

\(n = 4; 
Exercise[0, 3.5, upwind, 70, .8, 50]\)Exercise[0, 3.5, laxFriedrichs, 70, .8, 50]Exercise[0, 3.5, laxWendroff , 70, .8, 50]Exercise[0, 3.5, leapFrog , 80, .8, 50]

\(n = 2; 
CFL = 0.2; 
Exercise[0, 3.5, upwind, 70, CFL, 50]\)Exercise[0, 3.5, laxFriedrichs, 70, CFL, 50]Exercise[0, 3.5, laxWendroff , 70, CFL, 50]Exercise[0, 3.5, leapFrog , 80, CFL, 50]

Problem 8.16