Last week I showed how the accuracy, stability and general properties of an ODE integrator might be studied with the aid of Mathematica. This week I will do the same for a partial differential equations solution. Again, I will provide the commands used in Mathematica to conduct the analysis reported at the end of the post.

It is good to start as simple as possible. That was the reason for retreading the whole ODE stability analysis last week. Now we can steadily go forward toward looking at something a bit harder partial differential equations, starting with a first-order method for a first-order hyperbolic equation, the linear advection equation,

$u_t + u_x=0$, where the subscript denotes differentiation with respect to the variable. This equation is about as simple as PDEs get, but it is notoriously difficult to solve numerically.

Before getting to the analysis we can state a few properties of the equation. The exact solution is outrageously simple, $u \left(x, t \right) = u(x-t,0)$. This means that the temporal solution is simply defined by the initial condition translated by the velocity (which is one in this case) and time. Nothing changes it simply moves in space. This is a very simple form of space-time self-similarity. If we are solving this equation numerically, any change in the waveform is an error. We can also note that the integral of the value is preserved (of course) making this a “conservation law”. Later when you’d like to solve harder problems this property is exceedingly important.

Now we can proceed to the analysis. The basic process is to replace the function with an analytical representation and similar to ODEs we use the complex exponential (Fourier transform), $\exp\left(\imath j \theta\right)$, where $j$ is the grid index of our discretized function, and $\theta$ is the angle parameterizing frequency of the waveform. The analysis then proceeds much as in the style as the ODE work from last week, one substitutes this function into the numerical scheme and works out the modification of the waveform by the numerical method. We then take this modification to be the symbol of the operator $A\left(\theta\right) = \left| A \right\| \exp\left(\imath\alpha\right)$. In this form we have divided the symbol into two effects its amplitude and its modulation of the waveform or phase. Finishing our conceptual toolbox is the expression of the exact solution as $u\left(x,0\right)\exp\left(-\imath t \theta\right)$.

We are now ready to apply the analysis technique to the scheme. We can start off with something horribly simple like first-order upwind. The numerical method is easy to write down as $u_j^{n+1}=u_j^n-\nu\left(u_{j+1/2}^n-u_{j-1/2}^n\right)$ where $\nu= \Delta t / \Delta x$ is the Courant or CFL number and $u_{j+1/2}^n = u_j^n$ is the upwind edge value. The CFL number is the similarity variable (dimensionless) of greatest important for numerical schemes for hyperbolic PDEs. Now we plug our Fourier function into the grid values in the scheme and evaluate for a single grid point $j=0$. Without showing the trivial algebraic steps this gives $A = 1 - \nu\left(1-\exp(-\imath \theta)\right)$. We can make the substitution of the trigonometric functions for the complex exponential, $\exp\left(-\imath \theta\right) = \cos\left(\theta\right) – \imath \sin\left(\theta\right)$.

Now it is time to use these relations to provide the properties of the numerical scheme. We will divide these effects into two categories, changes in the amplification of the function that will define stability, $\latex \left| A \right|$ and the phase error $\alpha$. The exact solution has amplitude of one, and a phase of $\nu \theta$. Once we have separated the symbol into its pieces we can then examine the formal truncation error of the method (as $\theta\rightarrow 0$ is equivalent to $\Delta x\rightarrow 0$) in a straightforward manner.

We can also expand these in a Taylor series to get a result for the truncation error. For the amplitude we get the following $\left|A\right\| \approx 1 -\frac{1}{2} \left(\nu-\nu^2 \right)\theta^2 + O\left(\theta^4\right)$. The phase error can be similarly treated, $\alpha \approx 1 + \frac{1}{6}\left(1-2\nu + \nu^2\right) + O\left(\theta^4\right)$. Please note that the phase error is actually one order higher than I’ve written because of its definition where I have divided

through by $\nu\theta$. The last bit of analysis we conduct is to make an estimate of the rate of convergence as a function of the mesh spacing and CFL number. Given the symbol we can compute the error $E=A - \exp\left(-\imath \nu\theta\right)$. We then compute the error with a refined grid by a factor of two and note that it must applied twice to get the solution to the same point in time. The error for the refined calculation is $E_{\frac{1}{2}} = A_{\frac{1}{2}} - \exp \left( - \frac{\imath \nu \theta}{2} \right)$, which is squared to account for two time steps being taken to get to the same simulation time, $E_{\frac{1}{2}}:=E_{\frac{1}{2}}^2$.Given these errors the local rate of convergence is simple, $n = \log\left(\left|E\right|/\left|E_\frac{1}{2}\right| \right)/log\left(2\right)$. We can then plot the function where we see that the convergence rate deviates significantly from one (the expected value) for finite values of $\theta$ and $\nu$.

We can now apply the same machinery to more complex schemes. Our first example is the time-space coupled version of Fromm’s scheme, which is a second-order method. Conducting the analysis is largely a function of writing the numerical scheme in Mathematica much in the same fashion we would use to write the method into a computer code.

The first version of Fromm’s scheme uses a combined space time differencing introduced by Lax-Wendroff implemented using a methodology similar to Richtmyer’s two-step scheme, which makes the steps clear. First, define a cell-centered slope $s_j^n =\frac{1}{2}\left( u_{j+1}^n - u_{j-1}^n\right)$ and then use this to define a edge-centered, time-centered value, $u_{j+1/2}^{n+1/2} = u_j^n + \frac{1}{2}\left(1 - \nu\right) s_j^n$. This choice has a “build-in” upwind bias. If the velocity in the equation were oriented oppositely, this choice would be $u_{j+1/2}^{n+1/2} = u_{j+1}^n - \frac{1}{2}\left(1 - \nu\right)s_{j+1}^n$ instead ($\nu<0$). Now we can write the update for the cell-centered variables as $u_j^{n+1} = u_j^n - \nu\left(u_{j+1/2}^{n+1/2} - u_{j-1/2}^{n+1/2}\right)$, substitute in the Fourier transform and apply all the same rules as for the first-order upwind method.

Just note that in the Mathematica the slope, and edge variables are defined as general functions of the mesh index $j$ and the substitution is accomplished without any pain. This property is essential for analyzing complicated methods that effectively have very large or complex stencils.

The results then follow as before. We can plot the amplitude and phase error easily and the first thing we should notice is the radical improvement over the first-order method, particularly the amplification error at large wavenumbers (i.e., the grid scale). We can go further and use the Taylor series expansion to express the formal accuracy for the amplification and phase error. The amplification error is two orders higher than upwind and is $\left| A \right| \approx 1$. The phase error is smaller than the upwind scheme, but the same order, $\alpha\approx 1$. This is the leading order error in Fromm’s scheme.

We can finish by plotting the convergence rate as a function of finite time step and wavenumber. Unlike the upwind scheme as the wavenumber approaches one the rate of convergence is larger than the formal order of accuracy.

The Mathematica commands used to conduct the analysis above:

(* 1st order 1-D *)

U[j_] := T[j]

U1[j_] := U[j] – v (U[j] – U[j – 1])

sym = 1/2 U[0] + 1/2 U1[0] – v/2 (U1[0] – U1[-1]);

T[j_] := Exp[I j t]

Series[sym – Exp[-I v t], {t, 0, 5}]

Simplify[sym]

Sym[v_, t_] := 1/2 E^(-2 I t) (-2 E^(I t) (-1 + v) v + v^2 + E^(2 I t) (2 – 2 v + v^2))

rg1 = Simplify[ComplexExpand[Re[sym]]];

ig1 = Simplify[ComplexExpand[Im[sym]]];

amp1 = Simplify[Sqrt[rg1^2 + ig1^2]];

phase1 = Simplify[ ArcTan[-ig1/rg1]/(v t)];

Series[amp1, {t, 0, 5}]

Series[phase1, {t, 0, 5}]

Plot3D[amp1, {t, 0, Pi}, {v, 0, 1}, AxesLabel -> {t, v}, LabelStyle -> Directive[18, Bold, Black]]

Export[“amp-1.jpg”, %]

ContourPlot[amp1, {t, 0, Pi}, {v, 0, 1}, PlotPoints -> 250, Contours -> {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.99, 1}, ContourShading -> False, Axes -> {True, True}, AxesLabel -> {t, v}, ContourLabels -> All, LabelStyle -> Directive[18, Bold, Black]]

Export[“amp-1-cont.jpg”, %]

Plot3D[phase1, {t, 0, Pi}, {v, 0, 1}, AxesLabel -> {t, v}, LabelStyle -> Directive[18, Bold, Black]]

Export[“phase-1.jpg”, %]

ContourPlot[phase1, {t, 0, Pi}, {v, 0, 1}, PlotPoints -> 250, Contours -> {-1, -0.5, -0.25, 0, 0.25, 0.5, 0.75, 0.9, 0.99, 1, 1.25, 1.5}, ContourShading -> False, Axes -> {True, True}, AxesLabel -> {t, v}, ContourLabels -> All, LabelStyle -> Directive[18, Bold, Black]]

Export[“phase-1-cont.jpg”, %]

err = Sym[v, t] – Exp[-I v t];

err2 = Sym[v, t/2]^2 – Exp[-I v t];

Plot3D[Log[Abs[err]/Abs[err2]]/Log[2], {t, 0.01, Pi/2}, {v, 0, 1}, AxesLabel -> {t, v}, LabelStyle -> Directive[18, Bold, Black]]

Export[“conv-rate-1.jpg”, %]

ContourPlot[Log[Abs[err]/Abs[err2]]/Log[2], {t, 0, Pi}, {v, 0, 1}, PlotPoints -> 250, Contours -> {0.5, 0.75, 0.9, 0.95, 0.99}, ContourShading -> False, Axes -> {True, True}, AxesLabel -> {t, v}, ContourLabels -> All, LabelStyle -> Directive[18, Bold, Black]]

Plot3D[Abs[sym/Exp[-I v t]], {t, 0, Pi}, {v, 0, 5}]

ContourPlot[ If[Abs[sym/Exp[-I v t]] <= 1, Abs[sym/Exp[-I v t]], 0], {t, 0, Pi}, {v, 0, 5}, PlotPoints -> 250, Contours -> {0.1, 0.25, 0.5, 0.75, 0.9, 1}, ContourShading -> False, Axes -> {True, True}, AxesLabel -> {t, v}, ContourLabels -> All, LabelStyle -> Directive[18, Bold, Black]]

errs = Sym[v, t/2]^2 – Sym[v, t];

errs2 = Sym[v, t/4]^4 – Sym[v, t/2]^2;

Plot3D[Log[Abs[errs]/Abs[errs2]]/Log[2], {t, 0.01, Pi/2}, {v, 0, 1}]

errt = Sym[v/2, t]^2 – Sym[v, t];

errt2 = Sym[v/4, t]^4 – Sym[v/2, t]^2;

Plot3D[Log[Abs[errt]/Abs[errt2]]/Log[2], {t, 0, Pi/2}, {v, 0, 1}]

(* classic fromm *)

U[j_] := T[j]

S[j_] := 1/2 (U[j + 1] – U[j – 1])

Ue[j_] := U[j] + 1/2 S[j] (1 – v)

sym2 = U[0] – v (Ue[0] – Ue[-1]);

T[j_] := Exp[I j t]

Series[sym2 – Exp[-I v t], {t, 0, 5}];

Simplify[Normal[%]];

Collect[Expand[Normal[%]], t]

Simplify[sym2]

Sym2[v_, t_] := 1/32 E^(-4 I t) (v^2 – 10 E^(I t) v^2 + E^(6 I t) v^2 + 2 E^(5 I t) v (-4 + 3 v) – 4 E^(3 I t) v (-10 + 7 v) + E^(2 I t) v (-8 + 31 v) – E^(4 I t) (-32 + 24 v + v^2))

rg2 = Simplify[ComplexExpand[Re[sym2]]];

ig2 = Simplify[ComplexExpand[Im[sym2]]];

amp2 = Simplify[Sqrt[rg2^2 + ig2^2]];

phase2 = Simplify[ ArcTan[-ig2/rg2]/(v t)];

Series[amp2, {t, 0, 5}]

Series[phase2, {t, 0, 5}]

Plot3D[amp2, {t, 0, Pi}, {v, 0, 1}, AxesLabel -> {t, v}, LabelStyle -> Directive[18, Bold, Black]]

Export[“amp-2.jpg”, %]

ContourPlot[amp2, {t, 0, Pi}, {v, 0, 1}, PlotPoints -> 250, Contours -> {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.99, 1}, ContourShading -> False, Axes -> {True, True}, AxesLabel -> {t, v}, ContourLabels -> All, LabelStyle -> Directive[18, Bold, Black]]

Export[“amp-2-cont.jpg”, %]

Plot3D[phase2, {t, 0, Pi}, {v, 0, 1}, AxesLabel -> {t, v}, LabelStyle -> Directive[18, Bold, Black]]

Export[“phase-2.jpg”, %]

ContourPlot[phase2, {t, 0, Pi}, {v, 0, 1}, PlotPoints -> 250, Contours -> {-1, -0.5, -0.25, 0, 0.25, 0.5, 0.75, 0.9, 0.99, 1, 1.25, 1.5}, ContourShading -> False, Axes -> {True, True}, AxesLabel -> {t, v}, ContourLabels -> All, LabelStyle -> Directive[18, Bold, Black]]

Export[“phase-2-cont.jpg”, %]

err = Sym2[v, t] – Exp[-I v t];

err2 = Sym2[v, t/2]^2 – Exp[-I v t];

Plot3D[Log[Abs[err]/Abs[err2]]/Log[2], {t, 0.01, Pi/2}, {v, 0, 1}, AxesLabel -> {t, v}, LabelStyle -> Directive[18, Bold, Black]]

Export[“conv-rate-2.jpg”, %]

ContourPlot[ Log[Abs[err]/Abs[err2]]/Log[2], {t, 0.01, Pi}, {v, 0.01, 1}, PlotPoints -> 250, Contours -> {1, 1.5, 1.75, 1.9, 2, 2.1, 2.2, 2.3, 2.4}, ContourShading -> False, Axes -> {True, True}, AxesLabel -> {t, v}, ContourLabels -> All, LabelStyle -> Directive[18, Bold, Black]]

Export[“conv-rate-2-cont.jpg”, %]

ContourPlot[ If[Abs[sym2/Exp[-I v t]] <= 1, Abs[sym2/Exp[-I v t]], 0], {t, 0, Pi}, {v, 0, 1}, PlotPoints -> 250, Contours -> {0.1, 0.25, 0.5, 0.75, 0.9, 1}, ContourShading -> False, Axes -> {True, True}, AxesLabel -> {t, v}, ContourLabels -> All, LabelStyle -> Directive[18, Bold, Black]]

errs = Sym2[v, t] – Sym2[v, t/2]^2;

errs2 = Sym2[v, t/4]^4 – Sym2[v, t/2]^2;

Plot3D[Log[Abs[errs]/Abs[errs2]]/Log[2], {t, 0.01, Pi/2}, {v, 0, 1}, AxesLabel -> {t, v}, LabelStyle -> Directive[18, Bold, Black]]

errt = Sym2[v, t] – Sym2[v/2, t]^2;

errt2 = Sym2[v/4, t]^4 – Sym2[v/2, t]^2;

Plot3D[Log[Abs[errt]/Abs[errt2]]/Log[2], {t, 0.0, Pi/2}, {v, 0, 0.1}]

(* 2nd order Fromm – RK *)

U[j_] := T[j]

S[j_] := 1/2 (U[j + 1] – U[j – 1])

Ue[j_] := U[j] + 1/2 S[j]

U1[j_] := U[j] – v (Ue[j] – Ue[j – 1])

S1[j_] := 1/2 (U1[j + 1] – U1[j – 1])

Ue1[j_] := U1[j] + 1/2 S1[j]

sym2 = 1/2 U[0] + 1/2 U1[0] – v/2 (Ue1[0] – Ue1[-1]);

T[j_] := Exp[I j t]

Series[sym2 – Exp[-I v t], {t, 0, 5}];

Simplify[Normal[%]];

Collect[Expand[Normal[%]], t]

Simplify[sym2]

Sym2[v_, t_] := 1/32 E^(-4 I t) (v^2 – 10 E^(I t) v^2 + E^(6 I t) v^2 + 2 E^(5 I t) v (-4 + 3 v) – 4 E^(3 I t) v (-10 + 7 v) + E^(2 I t) v (-8 + 31 v) – E^(4 I t) (-32 + 24 v + v^2))

rg2 = Simplify[ComplexExpand[Re[sym2]]];

ig2 = Simplify[ComplexExpand[Im[sym2]]];

amp2 = Simplify[Sqrt[rg2^2 + ig2^2]];

phase2 = Simplify[ ArcTan[-ig2/rg2]/(v t)];

Series[amp2, {t, 0, 5}]

Series[phase2, {t, 0, 5}]

Plot3D[amp2, {t, 0, Pi}, {v, 0, 1}, AxesLabel -> {t, v}, LabelStyle -> Directive[18, Bold, Black]]

Export[“amp-2rk.jpg”, %]

ContourPlot[amp2, {t, 0, Pi}, {v, 0, 1}, PlotPoints -> 250, Contours -> {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.99, 1}, ContourShading -> False, Axes -> {True, True}, AxesLabel -> {t, v}, ContourLabels -> All, LabelStyle -> Directive[18, Bold, Black]]

Export[“amp-2rk-cont.jpg”, %]

Plot3D[phase2, {t, 0, Pi}, {v, 0, 1}, AxesLabel -> {t, v}, LabelStyle -> Directive[18, Bold, Black]]

Export[“phase-2rk.jpg”, %]

ContourPlot[phase2, {t, 0, Pi}, {v, 0, 1}, PlotPoints -> 250, Contours -> {-1, -0.5, -0.25, 0, 0.25, 0.5, 0.75, 0.9, 0.99, 1, 1.25, 1.5}, ContourShading -> False, Axes -> {True, True}, AxesLabel -> {t, v}, ContourLabels -> All, LabelStyle -> Directive[18, Bold, Black]]

Export[“phase-2rk-cont.jpg”, %]

err = Sym2[v, t] – Exp[-I v t];

err2 = Sym2[v, t/2]^2 – Exp[-I v t];

Plot3D[Log[Abs[err]/Abs[err2]]/Log[2], {t, 0.01, Pi/2}, {v, 0, 1}, AxesLabel -> {t, v}, LabelStyle -> Directive[18, Bold, Black]]

Export[“conv-rate-2rk.jpg”, %]

ContourPlot[ Log[Abs[err]/Abs[err2]]/Log[2], {t, 0.01, Pi}, {v, 0.01, 1}, PlotPoints -> 250, Contours -> {1, 1.5, 1.75, 1.9, 2, 2.1, 2.2, 2.3, 2.4}, ContourShading -> False, Axes -> {True, True}, AxesLabel -> {t, v}, ContourLabels -> All, LabelStyle -> Directive[18, Bold, Black]]

Export[“conv-rate-2rk-cont.jpg”, %]

ContourPlot[ If[Abs[sym2/Exp[-I v t]] <= 1, Abs[sym2/Exp[-I v t]], 0], {t, 0, Pi}, {v, 0, 1.5}, PlotPoints -> 250, Contours -> {0.1, 0.25, 0.5, 0.75, 0.9, 1}, ContourShading -> False, Axes -> {True, True}, AxesLabel -> {t, v}, ContourLabels -> All, LabelStyle -> Directive[18, Bold, Black]]

errs = Sym2[v, t] – Sym2[v, t/2]^2;

errs2 = Sym2[v, t/4]^4 – Sym2[v, t/2]^2;

Plot3D[Log[Abs[errs]/Abs[errs2]]/Log[2], {t, 0.01, Pi/2}, {v, 0, 1}]

errt = Sym2[v, t] – Sym2[v/2, t]^2;

errt2 = Sym2[v/4, t]^4 – Sym2[v/2, t]^2;

Plot3D[Log[Abs[errt]/Abs[errt2]]/Log[2], {t, 0.0, Pi/2}, {v, 0, 0.1}]