Acta Polytechnica https://doi.org/10.14311/AP.2021.61.0049 Acta Polytechnica 61(SI):49–58, 2021 © 2021 The Author(s). Licensed under a CC-BY 4.0 licence Published by the Czech Technical University in Prague MODIFIED EQUATION FOR A CLASS OF EXPLICIT AND IMPLICIT SCHEMES SOLVING ONE-DIMENSIONAL ADVECTION PROBLEM Tomáš Bodnára, b, ∗, Philippe Frauniéc, Karel Kozela a Czech Technical University in Prague, Faculty of Mechanical Engineering, Karlovo Náměstí 13, 121 35 Prague, Czech Republic b Czech Academy of Sciences, Institute of Mathematics, Žitná 25, 115 67 Prague, Czech Republic c Université de Toulon, Mediterranean Institute of Oceanography - MIO, BP 20132 F-83957 La Garde cedex, France ∗ corresponding author: Tomas.Bodnar@fs.cvut.cz Abstract. This paper presents the general modified equation for a family of finite-difference schemes solving one-dimensional advection equation. The whole family of explicit and implicit schemes working at two time-levels and having three point spatial support is considered. Some of the classical schemes (upwind, Lax-Friedrichs, Lax-Wendroff) are discussed as examples, showing the possible implications arising from the modified equation to the properties of the considered numerical methods. Keywords: Modified equation, finite difference, advection equation. 1. Introduction Numerical solution of differential equations became a standard tool in many disciplines of theoretical science as well as in applied sciences and engineering. Numerical and computational methods brought the possibility to solve non-trivial problems described by ordinary and partial differential equations for which it was impossible to obtain an analytical solution by standard mathematical methods. A wide range of numerical methods was developed during the years for specific problems. Typically the physical problem is first described mathematically, so the mathematical model is created. Then this mathematical model is solved numerically. The obtained numerical solution is an approximation to the solution of the mathematical model, which itself is just an approximation of the physical problem. We keep aside the error introduced by the inaccuracies and approximations made when developing the math- ematical model from the physical problem. Here the focus is on the difference between the discrete numer- ical solution of the mathematical model and its exact (analytical) solution. The numerical solution strongly depends on the method used to obtain it. So the properties and quality of the numerical solution may differ from the solution of the original mathematical model. The aim of this paper is to show that the numerical solution of certain problems (equations) is much closer to the solution of some modified equation, rather than to the solution of the original problem. Many impor- tant properties of the numerical solution can than be easily seen (and a-priori expected) from the behavior of the known solution of the modified problem. This rather general principle will be demonstrated on a finite-difference approximation of the solution of one-dimensional advection equation. The structure of the paper is as follows. First the advection, diffusion and dispersion equations are in- troduced and discussed. Then several explicit schemes for advection equation are presented and analyzed at the discrete level. Modified equation is first derived for upwind scheme and then extended for a general class of explicit and implicit schemes. Finally the modified equations are discussed from the point of view of numerical diffusion and dispersion. 2. Model problems In order to be able to explain the properties of modi- fied equations and the underlying numerical schemes three model problems are introduced. They describe the physical phenomena of advection, diffusion and dispersion. All these problems can be mathemati- cally formulated using a linear evolutionary partial differential equation. In all cases the unknown func- tion u(x,t) is the sought subject to the initial data u(x,t = 0) = η(x). For a sketch of an example of initial data see Fig. 1. The initial value problem can be thus solved an- alytically using the Fourier transform method. The Fourier transform (in space) needed to obtain the analytical solutions is defined by: û(ξ,t) = 1 √ 2π ∫ +∞ −∞ u(x,t)e−iξxdx , while the inverse is u(x,t) = 1 √ 2π ∫ +∞ −∞ û(ξ,t)eiξxdξ . 49 https://doi.org/10.14311/AP.2021.61.0049 https://creativecommons.org/licenses/by/4.0/ https://www.cvut.cz/en T. Bodnár, Ph. Fraunié, K. Kozel Acta Polytechnica u x u(x,t )= (x)η0 Figure 1. Initial condition. Using this transformation, the analytical solutions of model linear problems (discussed further) can be derived. The details can be found for example in the appendix of the book [1]. 2.1. Advection equation The advection problem can be described by a first order PDE of the form: ut + aux = 0 In this paper we consider the advection velocity a > 0, but for a ≤ 0 the solution can be obtained as well. This advection equation is supplemented by the initial data u(x,t = 0) = η(x) . Applying the Fourier transform we can get the expres- sion for the Fourier image of the solution û(ξ,t) = e−iξatη̂(ξ) . Using the inverse transform, the solution can be writ- ten in the form u(x,t) = 1 √ 2π ∫ +∞ −∞ e−iξatη̂(ξ)eiξxdξ and finally u(x,t) = 1 √ 2π ∫ +∞ −∞ η̂(ξ)eiξ(x−at)dξ It’s not difficult to see that the solution corresponds in fact to the initial data η(x) shifted along the x-axis at velocity a, i.e. u(x,t) = η(x−at) . An illustration of such a situation can be seen in Fig. 2. When the advection equation is solved numerically the discrete solution differs from the exact one, de- pending on the numerical method used. The numerical solution often exhibits some non-physical oscillations (due to dispersion) or smeared solution gradients (due to diffusion). This is why the diffusion and dispersion equations are briefly presented hereafter. 0t=t +2 tδt=t + t0 δt=t 0 u x a Figure 2. Advection equation solution evolution. 0t=t δ0t=t +2 t t=t + t0 δu x Figure 3. Diffusion equation solution evolution. 2.2. Diffusion equation The diffusion equation contains second spatial deriva- tive, multiplied by a diffusion coefficient b. ut + buxx = 0 Also this linear PDE can easily be solved analytically using the Fourier transform to obtain û(ξ,t) = e−i 2ξ2btη̂(ξ) and after the inverse transform the solution u(x,t) = 1 √ 2π ∫ +∞ −∞ eξ 2btη̂(ξ)eiξxdξ . A straightforward comparison with the solution of the advection equation reveals that now the individual Fourier modes found in the initial data η(x) do not change their position. They only change their ampli- tude by the factor eξ 2bt, which means the exponential decay for b < 0 and growth for b > 0. This is why only the decaying case with b < 0 is usually physi- cally acceptable leading to an asymptotically stable evolution in time. It should also be noted that the decay depends on the square of wave number ξ and thus the rapidly oscillating modes decay much faster. An illustration of the diffusion process and diffusion equation solution is shown in Fig. 3. 2.3. Dispersion equation The linear dispersion equation can be written as ut + cuxxx = 0 where c is the dispersion coefficient. The Fourier image of the solution is û(ξ,t) = e−i 3ξ3ctη̂(ξ) 50 vol. 61 Special Issue/2021 Modified equation for explicit and implicit schemes from which follows the solution u(x,t) = 1 √ 2π ∫ +∞ −∞ eiξ 3ctη̂(ξ)eiξxdξ that can be rearranged to u(x,t) = 1 √ 2π ∫ +∞ −∞ η̂(ξ)eiξ(x+cξ 2t)dξ This form is very similar to the solution of the advec- tion equation u(x,t) = 1 √ 2π ∫ +∞ −∞ η̂(ξ)eiξ(x−at)dξ . Comparing the corresponding expressions in the ex- ponential, it can be seen that during the dispersive evolution the Fourier modes are also shifted along the x-axis, but each mode is shifted at a different velocity depending on the wave number as −cξ2. The minus sign indicates that for c > 0 the modes propagate against the sense of the x-axis. In general it means that the initial data are decomposed into individual modes and each of them propagates at a different speed, proportional to ξ2. 2.4. Advection-Diffusion-Dispersion equation The above described model problems involving the ad- vection, diffusion and dispersion can be combined into advection-diffusion or advection-dispersion equations. The solutions of these combined equations can also be derived analytically. They will share the proper- ties and behavior of both original equations solutions. These combinations will be important later in this paper while discussing the properties of the modified equations, where often the diffusive and dispersive term appear due to the discretization of advection equation. See the discussions in the sections 4.1–5. 3. Discrete analysis The problem of numerical solution of advection equa- tion is one of the classical topics in numerical math- ematics. There exists a wide range of numerical schemes to discretize this problem. Here we only show a few classical methods as an illustration. Further we will consider advection in one spatial dimension ut + aux = 0 a > 0 , where u(x,t) is the sought solution and uni ≈ u(xi, tn) is its discrete numerical approximation at point xi = i∆x and time instant tn = n∆t, i,n being integer in- dices in spatial and temporal coordinates respectively. The numerical approximation of the solution can be constructed e.g. using the finite-difference dis- cretization, replacing the derivatives in the original equation by the corresponding divided differences for- mulas developed from the Taylor expansions. The x t i−1 i+1i n+1 n Figure 4. Initial condition. final formulas for few classical schemes are shown in the Table 1. Some details concerning the process of derivation of some of these schemes will be shown in the following sections. More details can be found e.g. in the classical textbooks [2], [3] or [4]. 3.1. Up-wind & Down-wind decomposition A family of simple explicit schemes working on a three-point spatial stencil and two time levels (see Fig. 4) can be derived using forward and backward difference approximations of the spatial derivative. When a general (convex) combination of forward and backward differences is used, with weighting factor α, the explicit finite difference approximation will take the form un+1i −u n i ∆t + a [ α ( uni+1 −u n i ∆x ) + +(1 −α) ( uni −u n i−1 ∆x )] = 0 This family of schemes can be formally written using the forward (down-wind) difference {D} and backward (up-wind) difference {U} as un+1i −u n i ∆t + a [ α{D} + (1 −α){U} ] = 0 . It’s evident, that the classical upwind and downwind schemes can be recovered for special choices of α = 0 and α = 1 respectively. The simple central scheme can be obtained for a symmetric choice represented by α = 1/2. A little bit less apparent, but still easy to verify is the representation of Lax-Friedrichs and Lax-Wendroff schemes, that also belong to this family. Values of the upwind/downwind blending coefficient α for all these schemes are summarized in the Tab. 2. 3.2. Central & Upwind decomposition In the same way as every explicit scheme with a three- point stencil can be written in the form of convex combination of forward and backward differences, it can also be formally rewritten as a weighted sum of central and upwind difference approximations. The central approximation {C} can simply be expressed as an average of upwind and downwind (backward and forward) differences as {C} = {D} + {U} 2 =⇒ {D} = 2{C}−{U} . 51 T. Bodnár, Ph. Fraunié, K. Kozel Acta Polytechnica Scheme Formula Up-wind [U] un+1i = u n i − a∆t ∆x (u n i −u n i−1) Down-wind [D] un+1i = u n i − a∆t ∆x (u n i+1 −u n i ) Central [C] un+1i = u n i − a∆t 2∆x(u n i+1 −u n i−1) Lax-Friedrichs [LF] un+1i = 1 2 (u n i+1 + u n i−1) − a∆t 2∆x(u n i+1 −u n i−1) Lax-Wendroff [LW] un+1i = u n i − a∆t 2∆x(u n i+1 −u n i−1) + a2∆t2 2∆x2 (u n i+1 − 2u n i + u n i−1) Table 1. Examples of some simple classical discretizations for advection equation. Scheme Coefficient α [U] 0 [D] 1 [C] 12 [LF] 12 − ∆x 2a∆t [LW] 12 − a∆t 2∆x Table 2. Coefficient of up-wind/down-wind decom- position of classical schemes. The family of explicit schemes can thus be rewrit- ten using the central difference approximation {C} and additional upwinding {U}, where the blending parameter now has the value 2α. un+1i −u n i ∆t + a [2α{C} + (1 − 2α){U}] = 0 3.3. Central & Viscous decomposition By taking a divided difference of forward and back- ward differences (approximations of ux), the central approximation of the second derivative uxx can be constructed. This is often used in approximation of diffusive (or viscous terms) containing second deriva- tives. This viscous-like term {V} can be formally written as {V} = {D}−{U} ∆x =⇒ {U} = {C}− ∆x 2 {V} . Using this definition of {V}, the whole family of schemes can be written as a combination of the central approximation part and a viscous (diffusive) part. un+1i −u n i ∆t + a [ {C}− (1 − 2α) a∆x 2 {V} ] = 0 3.3.1. Numerical viscosity By a simple rearrangement of this formula, it’s easy to see that all the schemes from the considered explicit family can be written in the form, where only the central approximation of the first derivative is kept on left, while the remaining viscous-like part is moved to the right-hand side un+1i −u n i ∆t + a{C} = (1 − 2α) a∆x 2 {V} , Scheme Coefficient α � µ [U] 0 1 a∆x2 [D] 1 −1 −a∆x2 [C] 12 0 0 [LF] 12 − ∆x 2a∆t 1 γ ∆x2 2∆t [LW] 12 − a∆t 2∆x γ a2∆t 2 Table 3. Numerical viscosity coefficients for some classical schemes. or shortly un+1i −u n i ∆t + a{C} = µ{V} . It’s not difficult to see that this discrete equation corresponds to a finite-difference approximation of the advection-diffusion equation rather than just to the original advection equation. ut + aux = 0 −→ ut + aux = µuxx The extra term on the right hand side corresponds to the numerical diffusion (viscosity), where the coeffi- cient µ depends on the method being used. µ = (1 − 2α)︸ ︷︷ ︸ � a∆x 2 Values of this numerical viscosity coefficient for few classical explicit methods are listed in the table 3. In this context, each scheme within this explicit fam- ily can be considered as a simple central scheme with different amounts of added numerical viscosity (or up- winding, alternatively). When taking into account the definition of the non-dimensional Courant-Friedrichs- Lewy parameter γ = a∆t∆x , which is positive (and bounded by stability conditions), the non-dimensional viscosity coefficient � can be defined. The value of the numerical viscosity coefficient de- termines the essential behavior and properties of the specific numerical method. Just by changing the pa- rameter µ, the schemes can become more diffusive or dispersive, stable or unstable and also their order of accuracy will change. These properties for some schemes are summarized in the table 4. 52 vol. 61 Special Issue/2021 Modified equation for explicit and implicit schemes Scheme � µ Accuracy Stability Behavior [D] −1 −a∆x2 < 0 (∆t) 1/(∆x)1 Unstable Dispersive [C] 0 0 (∆t)1/(∆x)2 Unstable Dispersive [LW] γ < 1 a 2∆t 2 (∆t) 2/(∆x)2 Stable Dispersive [U] 1 a∆x2 (∆t) 1/(∆x)1 Stable Diffusive [LF] 1 γ > 1 ∆x 2 2∆t (∆t) 1/(∆x)1 Stable Diffusive Table 4. Properties and behavior of selected explicit numerical schemes. The schemes in the table 4 are sorted from top to bottom according to increasing numerical viscosity. It can be seen (as expected) that when the numerical viscosity coefficient is negative, the scheme is unstable. By increasing its value, the scheme becomes stable and also can improve its accuracy. However by further increase of the numerical viscosity, the behavior of the method becomes more diffusive and the formal order of accuracy drops down again. In summary, the identification of the amount of numerical viscosity embedded in a numerical scheme is the key point in understanding its behavior. Here it was done at the discrete level, by identifying the discrete diffusive term in the scheme. Similarly even more detailed analysis can be performed at the con- tinuous level, leading to so called modified equation. 4. Modified equation The modified equation approach is well known, often used to assess the order of accuracy for finite-difference schemes. When doing the Taylor expansions to de- velop the finite-difference approximations, it is possi- ble to go beyond just finding the order of the leading term in the truncated Taylor series. It’s possible to find analytically the form of the leading term, add it to (keep it in) the original equation and study the behavior of the modified equation that includes this extra term introduced by the discretization of the problem. The properties of the modified problem so- lution will be close to the properties of the numerical solution of the original problem. 4.1. Modified equation for upwind scheme When solving the the advection equation ut + aux = 0 a > 0 the Up-wind scheme can be written as un+1i −u n i ∆t + a [ uni −u n i−1 ∆x ] = 0 , using the explicit Euler discretization in time and backward difference in space (i.e. upwind when the advection velocity a > 0). Considering a sufficiently smooth interpolant u(xi, tn) = uni , the Taylor expansions can be used to derive the corresponding approximation formulas. u(xi, tn+1)−u(xi, tn) ∆t +a [ u(xi, tn)−u(xi−1, tn) ∆x ] = 0 The terms u(xi, tn+1) and u(xi, tn+1) appearing in this formula are then obtained from (truncated) the Taylor series u(xi, tn+1) = u(xi, tn) + ∆tut(xi, tn)+ + ∆t2 2 utt(xi, tn) + ∆t3 6 uttt(xi, tn) + O(∆t4) u(xi−1, tn) = u(xi, tn) − ∆xux(xi, tn)+ + ∆x2 2 uxx(xi, tn) − ∆x3 6 uxxx(xi, tn) + O(∆x4) . Using these values, the corresponding difference ap- proximations can be obtained for temporal and spatial derivatives. u(xi, tn+1) −u(xi, tn) ∆t = ut(xi, tn)+ + ∆t 2 utt(xi, tn) + ∆t2 6 uttt(xi, tn) + O(∆t3) u(xi, tn) −u(xi−1, tn) ∆x = ux(xi, tn)− − ∆x 2 uxx(xi, tn) + ∆x2 6 uxxx(xi, tn) + O(∆x3) When these difference approximations are put together as in the numerical scheme, the original advection equation appears on the right hand side, together with some extra terms, representing the leading order terms in the remainder of the truncated Taylor series. u(xi, tn+1)−u(xi, tn) ∆t +a [ u(xi, tn)−u(xi−1, tn) ∆x ] = = ut(xi, tn) + aux(xi, tn) + . . . . . . . . . = 0 It means that although the difference scheme approx- imates the advection equation, some extra (higher order) terms appear on the right hand side due to discretization. ut(xi, tn) + aux(xi, tn) = − ∆t 2 utt(xi, tn) − ∆t2 6 uttt(xi, tn) + O(∆t3)+ + a∆x 2 uxx(xi, tn) − a∆x2 6 uxxx(xi, tn) + O(∆x3) 53 T. Bodnár, Ph. Fraunié, K. Kozel Acta Polytechnica Partial derivatives with respect to time and also space, both appear on the right hand side. The derivatives with respect to time can be converted to spatial deriva- tives using the original advection equation (or corre- sponding scheme Taylor’s expansions): ut + aux = 0 =⇒ ∂ · ∂t = −a ∂ · ∂x utt = a2uxx & uttt = −a3uxxx This leads to ut(xi, tn) + aux(xi, tn) = − a2∆t 2 uxx(xi, tn) + a3∆t2 6 uxxx(xi, tn) + O(∆t3)+ + a∆x 2 uxx(xi, tn) − a∆x2 6 uxxx(xi, tn) + O(∆x3) . When only the leading order term is kept, the modified equation takes the form: ut(xi, tn) + aux(xi, tn) = = ( a∆x 2 − a2∆t 2 ) uxx(xi, tn) + O(∆t2; ∆x2) . In the approximate version, the extra term containing second (spatial) derivative uxx appears on the right hand side. ut(xi, tn) + aux(xi, tn) .= a∆x 2 (1 −γ)uxx(xi, tn) It means that while solving the advection equation by the Up-wind scheme, we obtain rather the solution of the advection-diffusion equation with the diffusion coefficient corresponding to the numerical viscosity µ ut + aux = 0 −→ ut + aux = a∆x 2 (1 −γ)︸ ︷︷ ︸ µ uxx In more detail, instead of the first order approxima- tion of the advection equation we have obtained a second order approximation of the advection-diffusion equation. 1st order approximation of original equation ut + aux = 0 + O(∆t; ∆x) 2nd order approximation of modified equation ut + aux = a∆x 2 (1 −γ)uxx + O(∆t2; ∆x2) So the numerical solution of the original (advection) problem is in fact much closer to the solution of the modified (advection-diffusion) problem, with all the consequences it may have on the solution behavior. From the form of the modified equation we can see the order of accuracy of the approximation of the original problem, the diffusive (or possibly dispersive) x t i−1 i+1i n+1 n α 4 α α α1 3 2 Figure 5. Computational stencil for general implicit and explicit schemes. character of the modified equation, representing the behavior of the numerical solution. It’s also possible to see e.g. how the numerical viscosity coefficient µ scales with time-step ∆t and spatial step ∆x. From the requirement of positivity of the diffusion coefficient µ > 0 it’s possible to estimate the stability of the underlying numerical scheme, i.e. the limitations for the CFL parameter γ (evidently γ < 1 is required for the Up-wind scheme). 4.2. General modified equation The procedure of deriving the modified equation can be applied to all explicit schemes discussed so far. In fact it can be applied to a much larger family including also implicit schemes. Further we will work with a family of explicit and implicit schemes, where the approximation of the spatial derivative is obtained as a linear combination of the forward and backward differences at time levels n and n + 1. The whole family of such schemes can be written in the form: un+1i −u n i ∆t + + a [ α1 ( uni+1 −u n i ∆x ) + α2 ( uni −u n i−1 ∆x )] + +a [ α3 ( un+1i+1 −u n+1 i ∆x ) +α4 ( un+1i −u n+1 i−1 ∆x )] = 0 . The computational stencil for such family of schemes is shown in the Fig. 5, introducing also the weights α1 – α4 used for blending the individual differences in the linear combination. Due to consistency the condition α1 + α2 + α3 + α4 = 1 should be verified, moreover obviously |α3| + |α4| = 0 for explicit schemes, while the opposite |α3| + |α4| 6= 0 characterizes the implicit schemes. For simplicity, the coefficients of the explicit part are marked in blue, while the implicit part is green. The coefficients α1 – α4 for some examples of clas- sical explicit and implicit schemes are shown in the Tab. 5. The last two rows in the table correspond to the Wendroff scheme denoted by [W] and Crank- Nicolson scheme denoted by [CN]. These two schemes are examples of implicit schemes, using the difference approximations at both time levels n and n + 1. For 54 vol. 61 Special Issue/2021 Modified equation for explicit and implicit schemes Scheme α1 α2 α3 α4 [U] 0 1 0 0 [D] 1 0 0 0 [C] 12 1 2 0 0 [LF] 12 − ∆x 2a∆t 1 2 + ∆x 2a∆t 0 0 [LW] 12 − a∆t 2∆x 1 2 + a∆t 2∆x 0 0 [W] 12 0 0 1 2 [CN] 14 1 4 1 4 1 4 Table 5. Coefficients of some explicit and implicit schemes with three-point stencil. Scheme Modified equation [U] ut + aux = (1 −γ)a∆x2 uxx [D] ut + aux = −(1 + γ)a∆x2 uxx [C] ut + aux = −γa∆x2 uxx [LF] ut + aux = ( 1γ −γ) a∆x 2 uxx [LW] ut + aux = −(1 −γ2)a∆x 2 6 uxxx [W] ut + aux = −(2 + 3γ + γ2)a∆x 2 12 uxxx [CN] ut + aux = −(2 + γ2)a∆x 2 12 uxxx Table 6. Modified equations for some explicit and implicit schemes. each of these schemes it’s possible to derive the modi- fied equation similarly as in the case of the Up-wind scheme. These modified equations are listed in the Tab. 6 Using the same procedure a general modified equa- tion (with up to 3rd order terms) for the whole family of considered schemes can be obtained in the form: ut + aux = �2uxx + �3uxxx . This is formally an advection-diffusion-dispersion equation with the numerical diffusion coefficient �2 and dispersion coefficient �3. These coefficients can be derived analytically to the form of functions depending on the blending weights α1 – α4. �2 = − a∆x 2 { (α1 −α2) + (α3 −α4)+ + γ [ (α1 + α2)2 − (α3 + α4)2 ]} �3 = − a∆x2 6 { 1 + 3γ [ (α21 −α 2 2) − (α 2 3 −α 2 4) ] + + 2γ2 [ (α1 + α2)3 + (α3 + α4)3 ]} 5. Numerical diffusion and dispersion The general modified equation for the whole ex- plicit/implicit family of schemes can be rewritten in 0 0.5 1 1.5 2 0.8 1 1.2 1.4 1.6 1.8 2 2.2 x−coordinate u ut + aux = (1 −γ) a∆x 2 uxx Figure 6. Up-wind [U] scheme solution and modified equation. the form ut + aux = �2uxx + �3uxxx = = �̃2 a∆x 2 uxx + �̃3 a∆x2 6 uxxx where the non-dimensional coefficients �̃2 and �̃3 are functions of the CFL parameter γ. This dependence of numerical diffusion and dispersion coefficients is obvious from the Tab. 7. As already noted in the case of the discrete analysis of classical schemes (end of section 3), the amount of added numerical viscosity (diffusion) is responsible for the essential properties and behavior of each scheme. 5.1. Diffusive schemes The schemes, for which the leading order term on the right-hand side of the modified equation is the diffusive term �̃2 a∆x2 uxx are of the first order of accuracy. If the sign of the numerical viscosity coefficient �̃2 is positive, the scheme is stable and it’s behavior is diffusive. The negative sign of �̃2 indicates that the scheme is unstable. The typical behavior of first order schemes is shown in the Fig. 6 for Up-wind and in Fig. 7 for Lax- Friedrichs scheme. The same initial value problem is solved for piece-wise constant initial data. The analytical solution corresponds to the "shifted" original data. The discrete numerical solution obtained using each scheme is marked by circles in the corresponding plots. While the sign of �̃2 determines the stability of the numerical method, the size of �̃2 is responsible for the amount of added numerical diffusion. The higher the coefficient, the more diffused (smeared) the solution will be. By comparing the Lax-Friedrichs and the upwind schemes modified equations it is clear that for given CFL parameter γ the numerical viscosity coefficient �̃2 = (1 γ −γ ) of the [LF] scheme is always greater than the �̃2 = (1−γ) of the [U] scheme. So the Lax-Friedrichs will be more diffusive than the Up-wind scheme. This is also well visible in the corresponding numerical solutions in the Figs. 6 and 7. 55 T. Bodnár, Ph. Fraunié, K. Kozel Acta Polytechnica Scheme �̃2 �̃3 Behavior [U] 1 −γ · · · �2 > 0 Diffusive [D] −(1 + γ) · · · �2 < 0 Unstable [C] −γ · · · �2 < 0 Unstable [LF] 1 γ −γ · · · �2 > 0 Diffusive [LW] 0 −(1 −γ2) �3 6= 0 Dispersive [W] 0 −12 (2 + 3γ + γ 2) �3 6= 0 Dispersive [CN] 0 −12 (2 + γ 2) �3 6= 0 Dispersive Table 7. Coefficients of numerical diffusion and dispersion for selected schemes. 0 0.5 1 1.5 2 0.8 1 1.2 1.4 1.6 1.8 2 2.2 x−coordinate u ut + aux = (1 γ −γ )a∆x 2 uxx Figure 7. Lax-Friedrichs [LF] scheme solution and modified equation. 5.2. Dispersive schemes For some schemes the first-order diffusive term van- ishes and the dominant role in the modified equation is played by the dispersive term �̃3 a∆x 2 6 uxxx. Due to this, such schemes are of a second order and their behavior is dispersive. The numerical solution be- haves much more like the solution of the advection- dispersion equation. It means that the advected initial data are decomposed into individual Fourier modes and each of them propagates at a different velocity, depending on the corresponding wave-number. This kind of behavior can be observed for Lax-Wendroff and Wendroff schemes in the Fig. 8 and 9 respectively. Again, for a given CFL parameter γ the numerical dis- persion coefficient �̃3 = −(1 −γ2) of the [LW] scheme is smaller (in the absolute value) than the coefficient �̃3 = −(2 + 3γ +γ2)/2 for the [W] scheme. This results in a less dispersive (less oscillatory) solution obtained using the Lax-Wendroff scheme. 6. Limitations and extensions In this paper the presentation was limited to the finite-difference approximation of linear advection, us- ing schemes with three-point stencil operating at two time levels. Most of these limiting assumptions can be removed or relaxed. Larger stencil and more time levels can be used to construct the scheme, it will only make the analysis of the scheme and its modified 0 0.5 1 1.5 2 0.8 1 1.2 1.4 1.6 1.8 2 2.2 x−coordinate u ut + aux = −(1 −γ2) a∆x2 6 uxxx Figure 8. Lax-Wendroff [LW] scheme solution and modified equation. 0 0.5 1 1.5 2 0.8 1 1.2 1.4 1.6 1.8 2 2.2 x−coordinate u ut + aux = −(2 + 3γ + γ2) a∆x2 12 uxxx Figure 9. Wendroff [W] scheme solution and modified equation. equation derivation more difficult (time consuming), however it can easily be done. The only important limitation is that the schemes should be linear in the sense that their coefficients can’t depend on the solu- tion. The non-linear schemes case will be significantly more complicated, although some kind of local lin- earization (frozen coefficients) would probably give at least some basic information. Multidimensional schemes can be analyzed in a very similar way, at least on regular, non-distorted grids. For finite-volume schemes such analysis can be per- formed on arbitrary grids, based on the expression and splitting of the numerical flux as a sum of the simple 56 vol. 61 Special Issue/2021 Modified equation for explicit and implicit schemes average central flux and a dissipative stabilizing part [5, 6]. Based on the extra dissipation, the properties of the scheme can be shown. 7. Remark on applications Based on a detailed knowledge of the diffu- sive/dispersive behavior of numerical schemes, suit- able strategy can be chosen to improve their properties, namely the stability and resolution. The classical way is to modify the embedded numerical viscosity of the scheme. It can be done for example by the following approaches: (1.) Modify the coefficients – As it was shown in sec- tion 3, each scheme can be written as a sum of the central (non-diffusive) part and additional in- ternal dissipative part (or sum of down-wind and up-wind parts alternatively). The coefficient of in- ternal numerical viscosity embedded in the scheme can be modified and adjusted by varying e.g. the blending parameter α and balancing the amount of forward and backward differences in the approxi- mation. For modified Lax-Friedrichs scheme with reduced numerical viscosity see e.g. [7]. (2.) Build a composite scheme – Two schemes are used, one with higher accuracy (typically disper- sive) and one with lower accuracy (typically dif- fusive). During the time-stepping process, several steps are performed using the higher order (but usu- ally oscillatory) scheme, followed by a "smoothing" step performed using a lower order diffusive scheme. The ratio of high/low order steps can be tuned to optimize the performance of the combined compos- ite method. An example of this technique using a combination of Lax-Wendroff and Lax-Friedrichs scheme can be found in [8]. (3.) Use artificial viscosity – The splitting of numeri- cal methods into the central and viscous (diffusive) part offers the possibility to use always the (accu- rate but unstable) central scheme, followed by a separate stabilizing step applying an artificial vis- cosity term. In this way, the numerical viscosity is separated and no longer relies on the form embed- ded (hidden) in the numerical discretization. The separate, stand alone, numerical viscosity can be designed and tuned for the specific problem and best performance. Typically the second and fourth order damping terms were used proportional to sec- ond and fourth order derivatives �2uxx and �4uxxxx. More efficient non-linear numerical viscosities can be constructed easily by considering variable, so- lution dependent, numerical viscosity coefficients �2(u), �4(u). When these coefficients are kept small on the smooth solution and only locally increased where the solution has high gradients or is oscilla- tory, a good resolution properties can be preserved while the stability and robustness of the method is improved. See e.g. [9] for artificial viscosity exam- ple and application or [10] or alternative filtering techniques. 8. Conclusions and Remarks The discrete analysis based on a decomposition of the scheme into a central and diffusive (or upwind) part was shown for a family of explicit schemes. A general modified equation was developed for a large family of explicit and implicit schemes. This led to a discussion concerning the diffusive and dispersive properties of numerical schemes. Most of the partial conclusions were already formulated in the above sections. Let’s just point out again some key points. The numerical solution of the advection equation is “much closer” to the solution of advection-diffusion or advection-dispersion equation, depending which is the leading order term in the discretization error. The behavior and quality of the numerical solution heavily depends on the coefficients of the modified equation. Their knowledge can help to assess a-priori the behavior of a numerical method. The detailed knowledge of the structure of the lead- ing order terms on the right-hand side of the modified equation can be used to construct the “high(er) reso- lution” numerical methods. Acknowledgements T. Bodnár is greatful for the support provided by the European Regional Development Fund-Project “Center for Advanced Applied Science” No.CZ.02.1.01/0.0/0.0/16 019/0000778 and partly by the Czech Science Foundation under the grant No. P201-19-04243S. References [1] R. J. LeVeque. Finite difference methods for ordinary and partial differential equations : steady-state and time-dependent problems. SIAM, 2007. [2] C. Hirsch. Numerical computation of internal and external flows, vol. 1,2. John Willey & Sons, 1988. [3] J. D. Anderson. Computational Techniques for Fluid Dynamics, vol. 1-2 of Springer Series in Computational Physics. Springer-Verlag Berlin Heidelberg, 2nd edn., 1991. [4] C. A. J. Fletcher. Computational Fluid Dynamics - The Basics with Applications. McGraw-Hill, 1995. [5] R. J. LeVeque. Numerical Methods for Conservation Laws. Lectures in Mathematics. Birkhäuser Verlag, 1990. [6] R. J. LeVeque. Finite Volume Methods for Hyperbolic Problems. Cambridge Texts in Applied Mathematics. Cambridge University Press, 2002. [7] R. Dvořák, K. Kozel. Mathematical modeling in aerodynamics (in Czech). Vydavatelství ČVUT, 1996. [8] R. Liska, B. Wendroff. Composite schemes for conservation laws. SIAM Journal on Numerical Analysis 35(6):2250–2271, 1998. doi:10.1137/S0036142996310976. [9] T. Bodnár, L. Beneš, K. Kozel. Numerical simulation of flow over barriers in complex terrain. Il Nuovo Cimento C 31(5–6):619–632, 2008. doi:10.1393/ncc/i2008-10323-4. 57 http://dx.doi.org/10.1137/S0036142996310976 http://dx.doi.org/10.1393/ncc/i2008-10323-4 T. Bodnár, Ph. Fraunié, K. Kozel Acta Polytechnica [10] A. Sequeira, T. Bodnár. On the filtering of spurious oscillations in the numerical simulations of convection dominated problems. Vietnam Journal of Mathematics 47:851–864, 2019. doi:10.1007/s10013-019-00369-z. 58 http://dx.doi.org/10.1007/s10013-019-00369-z Acta Polytechnica 61(SI):49–58, 2021 1 Introduction 2 Model problems 2.1 Advection equation 2.2 Diffusion equation 2.3 Dispersion equation 2.4 Advection-Diffusion-Dispersion equation 3 Discrete analysis 3.1 Up-wind & Down-wind decomposition 3.2 Central & Upwind decomposition 3.3 Central & Viscous decomposition 3.3.1 Numerical viscosity 4 Modified equation 4.1 Modified equation for upwind scheme 4.2 General modified equation 5 Numerical diffusion and dispersion 5.1 Diffusive schemes 5.2 Dispersive schemes 6 Limitations and extensions 7 Remark on applications 8 Conclusions and Remarks Acknowledgements References