Abstract
The fifth-order accurate Weighted Essentially Nonoscillatory space discretization developed by Jiang and Shu (WENO-JS) is studied theoretically. An exact Nonlinear Spectral Method (NSM) is developed based on an innovative yet simple methodology. The NSM explains the behaviour of nonsmooth solutions because it is valid for arbitrary modified wave numbers (MWN). The NSM clarifies the effects of the time integration methods and the Courant number. The mode isolation assumption, extensively used to analyse WENO-JS, is elucidated, and several novel findings are presented. The improved performance of the combination of WENO-JS with the forward Euler time integration method, compared to the Linear Fifth-Order Upwind discretization, is illustrated. The overdamping of the combination of WENO-JS with the popular third-order total variation diminishing Runge-Kutta method is discovered. Thus, the NSM covers several gaps in the current literature.
1. Introduction
Generally, high-order accurate discretizations suffer from spurious oscillations near high gradients due to Gibbs oscillations [1]. The fifth-order accurate Weighted Essentially Nonoscillatory method was developed by Jiang and Shu (WENO-JS) to achieve high-order accuracy while avoiding Gibbs oscillations [2]. The design of WENO-JS is based on a clever nonlinear modification of the fifth-order Linear Upwind (LUW5) discretization. The relatively wide five-point stencil of LUW5 can be subdivided into three small substencils. In WENO-JS, each smooth substencil is assigned a higher weight and vice versa. Hence, the constant weights adopted by LUW5 are replaced by nonlinear rational weights in WENO-JS.
A wide variety of WENO-JS improvements were developed, and a sample of the recent developments was provided in [3]. On the contrary, theoretical studies of WENO-JS are relatively rare. Spectral analysis of WENO-JS is a theoretical methodology that quantifies the solution in terms of dissipation and dispersion [4–7]. The spectral analysis methodologies used to analyse WENO-JS can be classified into two main categories: Geometric Stability Analysis (GSA) and Approximate Dispersion Relation (ADR).
The pioneering study of [5] adopted GSA to study the Linear Stability of WENO-JS and concluded that the forward Euler time integration method is linearly unstable. The results of [5] were revised by [6] who illustrated how linear stability can be satisfied because the discretization spectrum is equal to a finite set. However, the GSA provided by [5, 6] adopted linearized approximations of WENO-JS. For instance, many results were derived by [5] using Taylor expansions in terms of the modified wave number (MWN). Hence, these results were limited to small MWN that describe smooth solutions. This is a major disadvantage because the main goal of WENO-JS is modelling non-smooth solutions coexistent with high MWN. Moreover, the GSA of [6] was mainly based on LUW5. Hence, both [5, 6] did not identify the discrepancy between WENO-JS and LUW5 at high MWN.
The ADR is an alternative methodology based on discrete Fourier transform (DFT) and was originally developed by [4]. Within ADR, an initial condition defined on an arbitrary grid is advanced numerically in time. Next, DFT is applied to the solution to obtain the spectrum. The ADR algorithm accounts for nonlinearity and is valid for arbitrary MWN. The ADR was used to analyse the cost efficiency of WENO-JS [4], and a wide application was presented in [8]. The ADR was further developed by [7] and applied to analyse synthetic turbulent energy spectra [9, 10]. The design of the ADR algorithm requires adopting an arbitrarily small time step to minimize the effects of time integration. There are no definite criteria to select the numerical grid parameters within the ADR. In addition, neither the effects of the Courant number (CN) nor those of the numerical integration method were provided in [4, 7–10].
A cornerstone assumption that was adopted in the past GSA and ADR studies was ignoring the higher-order modes induced due to WENO-JS nonlinearity. However, theoretical validation of this assumption was not provided yet. This assumption will be denoted by mode isolation in this work. An exact methodology is proposed in this study to address the challenges mentioned above. A novel Nonlinear Spectral Method (NSM) is introduced to analyse WENO-JS. The benefits of the application of spectral methods to nonlinear problems are well explained in [11]. The development of the novel NSM is based on a continuous presentation of the discrete solution. This subtle and key idea does not compromise accuracy and proves to be effective within the analysis. Unlike [5, 6], NSM includes the full nonlinear rational weights of WENO-JS, and the results are valid for the full range of MWN. The NSM does not adopt any arbitrary grids, and the effects of CN and various time integration methods are described. A quantitative analysis of the mode isolation assumption is provided. To the author’s best knowledge, these contributions are provided for the first time.
2. Numerical Discretization
The numerical algorithms studied in this work will be detailed. Consider the scalar conservation law given by subject to the suitable boundary and initial conditions. Here, stands for a conserved quantity, and denotes the corresponding flux. Space and time domains are defined as and , respectively. To obtain a numerical solution of Equation (1) using the finite difference method, the space domain is discretized using the grid given by
Here, stands for the grid spacing, and stands for the total number of grid subdivisions. The numerical solution is defined as
2.1. Fifth-Order Space Discretization
The fifth-order accurate space discretization algorithm is explained in numerous sources and will be presented briefly. The reader may refer to [1] for more details. The main task is to obtain a high-order finite difference discretization of the flux space derivative . To simplify notations, let stand for . The space derivative of at will be written as
The attention is shifted towards the calculation of . This quantity is commonly termed as the numerical flux at point . Generally, is calculated based on the five-point stencil shown in Figure 1. This five-point stencil includes three smaller substencils as also illustrated in Figure 1. The three substencils are given by
Various choices exist to calculate . Consider the choice of that yields a third-order approximation given by
Here, the superscript corresponds to the substencil selected to calculate .
The linear combination of that is used to calculate is given by
The values of the constants are provided explicitly below for each .
The three choices for can be combined to provide a fifth-order accurate flux using the weighted sum given by
The selection of the weights will be detailed in the next subsections.
2.1.1. Fifth-Order Linear Upwind Discretization
Generally, are interpolated using a high-order polynomial. The following constants are defined [2, 12]:
The same values of are provided by [1] in the opposite order. The reason for this slight difference is that [1] listed the formulas to calculate the flux at point instead of . For the fifth-order Linear Upwind (LUW5) method, are given by
Hence, the fifth-order accurate linear numerical flux is given by
The algorithm of LUW5 is linear because constant values are assigned to , and accurate results may be obtained in smooth regions. However, the high-order polynomial used in LUW5 may produce spurious oscillations near high gradients due to the Gibbs phenomenon.
2.1.2. Fifth-Order Weighted Essentially Nonoscillatory Discretization
The fifth-order WENO-JS method is designed to overcome the disadvantages of LUW5 [1], and are designed to be higher for smooth and vice versa. Similar to [5, 6], the current analysis will focus on WENO-JS. Smoothness is quantified using the smoothness indicators () based on a sub-stencil . The norm squared was used to calculate in [2] using the following formula:
Here, is the derivative of the polynomial used to interpolate at . This definition of in Equation (13) yields the following expressions:
The modified weights are calculated as
Here, is an arbitrary small number introduced to avoid division by zero. The normalized weights are given by
Hence, for WENO-JS, the weights are given by and WENO-JS flux is given by
Substituting by (given by Equation (12)) or (given by Equation (18)) in Equation (4) results in LUW5 or WENO-JS space discretizations, respectively.
2.2. Numerical Time Integration
The operator is introduced to simplify notations. Referring to Equation (4), is defined as
The operator encapsulates the space discretization algorithm. The time integration algorithm is based on the method of lines [1, 4, 5]. A semidiscrete approximation of Equation (1) is given by
To perform the numerical integration of Equation (20), the time domain is discretized using given by
The time step is defined as . Numerical integration can be done using a variety of standard methods. The single-stage, first-order Euler forward solution of Equation (20) is written as follows [1]:
Here, and are the solution vectors at and , respectively. The symbol is introduced. And Equation (22) is rewritten as follows:
For a wave speed of unity, is commonly termed the Courant number (CN). The first-order Euler forward method can be regarded as a member of the set of Runge-Kutta (RK) methods [13]. Hence, the method defined by Equation (23) will be denoted by RK1. Higher-order multistage time integration algorithms can be implemented using the intermediate solution vectors and . The second-order explicit total variation diminishing Runge-Kutta (TVDRK2) solution of Equation (20) is given by
The third-order explicit total variation diminishing Runge-Kutta (TVDRK3) solution of Equation (20) is given by
The combination of WENO-JS space discretization with Equations (23), (24), and (25) will be denoted by WENO-RK1, WENO-TVDRK2, and WENO-TVDRK3, respectively. Similarly the combinations of LUW5 with Equations (23), (24), and (25) will be denoted by LUW5-RK1, LUW5-TVDRK2, and LUW5-TVDRK3, respectively.
3. Analysis
Similar to [4–6, 9], the analysis will focus on the linear flux . Hence, Equation (1) is rewritten as
The space periodic boundary condition is adopted. Consider the initial condition . The exact solution is obtained by the method of characteristics as follows:
Hence, the exact solution corresponds to a right-travelling wave whose form is identical to the initial condition and speed equal to 1.
3.1. Spectral Analysis of Numerical Solutions
Any periodic intermediate solution is exactly presented as a summation of Fourier modes.
The total number of modes is equal to the grid size . The modified wave number (MWN) termed as is given by
The operation defined by Equation (28) is termed the discrete Fourier transform (DFT). An extended explanation of DFT is provided in Appendix A. Generally, the range of is . However, for any real-valued periodic function, the total number of unique modes is and the rest of the modes are obtained by Hermitian symmetry [14]. Hence, the effective range of is reduced to be as adopted by [4, 8]. The Hermitian symmetric property of DFT of real functions is derived in Appendix B.
The numerical solution of any numerical algorithm can be written as
The spectral parameters and stand for the magnification factor and the phase shift, respectively. Generally, and are expressed as
The effect of the numerical algorithm on each mode is fully described in terms of and . A value of or indicates whether mode is magnified or damped, respectively. The significance of will be clarified as follows. Rewrite Equation (30) in terms of , and factor out .
Each mode in can be regarded as a wave, that travels at a speed . Hence, is obtained by shifting by a distance . Hence, Equation (33) is rewritten in the form of the travelling wave solution:
Comparing the arguments of the function in Equations (33) and (34), we get
Hence, is proportional to .
3.1.1. Mode Isolation
Consider a monochromatic intermediate solution that is described by a single mode [4] and defined as
For this monochromatic case, let us denote all wave parameters with subscript (e.g., ) as the principal mode parameters. It will be useful to write explicitly the exact solution for this case using Equation (27). First, Equation (36) is rewritten as
Next, the exact solution is obtained by replacing with and using the definition of Equation (2).
Finally, Equation (38) is rewritten in terms of the spectral parameters, introduced in Equation (30).
Upon comparison of Equations (38) and (39), the spectral parameters of the exact solution are given by
Let us consider the numerical solution . Based on Equation (36), the expression given by Equation (30) is unaltered. However, the expressions given by Equations (31) and (32) are simplified to be
As the accuracy of the numerical discretization is improved, the expressions provided by Equations (42) and (43) approach those provided by Equations (40) and (41), respectively.
So far, no limitations have been imposed on the presented analysis. However, it will be useful to introduce an assumption that was a cornerstone in [4–6]. This assumption will be denoted by the Mode Isolation or MI for abbreviation. This assumption is stated as follows.
Assumption 1. Let be a monochromatic intermediate discrete solution described by Equation (36). Then, resulting from Equation (23) is well approximated as
In other words, when a monochromatic is considered, the numerical solution can be equated to a single mode. Hence, all the summation terms in Equation (30) are equated to zero, except that corresponding to the principal mode. The validity of Assumption 1 will be thoroughly clarified in the next sections.
3.2. Linear von Neumann Analysis of LUW5
In this subsection, the Linear von Neumann (VN) analysis will be presented by application to the LUW5-RK1 algorithm. The VN analysis is generally used to study finite domain problems with periodic boundary conditions [15]. Adopting Equation (11) yields a linear scheme, and the operator is given by [6]
Alternatively, Equation (45) is rewritten in the following compact form.
Here, are constants whose values are listed explicitly in Equation (45). The numerical solution of LUW5-RK1 based on Equation (23) can be written as
The coefficients of on the right-hand side of Equation (47) are constants. Hence, Equation (47) is linear in terms of , and the analysis can be simplified as follows. The VN analysis considers a discrete Fourier Mode [16] that corresponds to the intermediate solution given by Equation (36). Although one mode is considered out of the multiple modes in Equation (28), the analysis is still general [17]. Also, adopting the complex notation is allowed, and Equation (36) can be rewritten in terms of as follows:
The expression of Equation (48) is substituted into Equation (47), and we get
To simplify Equation (49), the term is factored out to obtain
Here, , and the second equality is obtained by substitution by Equation (48). The spectral parameters are calculated as
Hence, the numerical solution of LUW5-RK1 is written as
The following points should be emphasized regarding the application of the VN method to LUW5-RK1 and other linear discretizations.
The expressions for and provided by Equations (51) and (52) are exactly valid for arbitrary values of and . (i)As illustrated by Equation (50), Assumption 1 is exactly satisfied by LUW5-RK1 and all linear discretizations(ii)Referring to Equations (51) and (52), the spectral parameters and are independent of and . Hence, Equations (42) and (43) are simplified to be
3.2.1. Extension to High-Order Time Integration
The analysis based on VN can be extended to high-order multistage time integration methods. It should be noted that the algorithms described in Equations (24) and (25) adopt Equation (23) as a building block. Hence, will be used frequently. Since LUW5 is linear, the following simplification of expressions including is justified:
The symbol in Equation (56) stands for any expression that does not depend explicitly on . The analysis will be based on adopting given by Equation (48). To analyse LUW5-TVDRK2, Equation (24) is rewritten as
Hence, for LUW5-TVDRK2, the spectral parameters are defined as
Similarly, LUW5-TVDRK3 is analysed using Equation (25). The derivation follows the same ideas listed in Equation (57). To avoid distraction, the derivation is provided in a relatively brief format.
Hence, for LUW5-TVDRK3, the spectral parameters are defined as
3.3. Nonlinear Spectral Analysis of WENO-JS
It should be noted that the numerical solution of Equation (26) using WENO-JS is nonlinear. The reason is that used in WENO-JS are not constants and are based on rational expressions in terms of . Hence, the steps explained in Section 3.2 cannot be applied directly to WENO-JS. Still, the nonlinear features of WENO-JS could be analysed using a slight modification of the VN method. This modification is based on a continuous presentation of WENO-JS. This presentation’s main benefit is clarifying two main solution properties: periodicity and homogeneity.
As explained in Section 3.2, the VN method is based on analysing Fourier modes for periodic finite domains. It will be shortly illustrated how the solution periodicity will be used to develop a nonlinear analysis that enjoys several advantages of the VN method. For instance, the grid parameters and will be eliminated and the analysis will depend only on and . The homogeneity of WENO-JS weights and flux will eliminate the dependence of the numerical solution on if the MI is adopted. It could be shown that the periodicity and the homogeneity properties are satisfied by several other nonlinear numerical algorithms. However, only WENO-JS will be considered in the current work.
3.3.1. Continuous Presentation of WENO-JS
The algorithm of WENO-JS is designed originally to deal with discrete quantities. However, WENO-JS algorithm steps can be applied to continuous quantities as follows. The continuous version of is obtained by substituting with the linear flux into Equation (14) and replacing by , where . Hence, based on Equation (14) the continuous version of is given by
The superscript in stands for continuous. The reader should note the following: (i)Although was skipped in Equation (61), time dependence is still included. In fact . This symbolic simplification is justified since operates only in space(ii) is defined for a continuous range of . This is in contrast to defined for the discrete points in defined by Equation (5)(iii) is equal to a second-order homogeneous polynomial [18] in terms of . Hence, can be regarded as a composite homogeneous function
Based on the discussion above, the following equation is introduced:
Here, stands for the composite function, and stands for the polynomials defined in the right-hand side of Equation (61). Since is a second-order homogeneous function, the following equation is satisfied:
Lemma 1. Let be a space periodic function with period . Then, is also space periodic with period .
Proof. By definition,
Substituting by into Equation (62), we get
The second equality is obtained based on Equation (64). Hence, is also periodic with period .
To proceed further, Equation (16) is rewritten as
Following the same direction, a continuous version of is defined as
A composite function form can be also used to describe as follows. Substituting Equation (62) into Equation (67), we get
Corollary 2. Let be a space periodic function with period . Then, is also space periodic with period .
Corollary 2 can be deduced using the same steps used to prove Lemma 1.
The value of should be proportional to the machine precision [12]. Hence, will be neglected and assigned to zero within the current exact analysis. Consequently, the expression for is written based on Equation (68) as
The second equality in Equation (69) is based on Equation (63). Hence, is a zero-order homogeneous function as illustrated by Equation (69) that. A continuous version of as defined by Equation (7) is given by
Similarly, a continuous version of defined by Equation (18) is given by
The flux can be regarded as the composite function given by
To clarify homogeneity, the flux is applied to as follows:
The second equality is based on Equations (68) and (70). Hence, is a first-order homogeneous function.
Corollary 3. Let be a space periodic function with period . Then, is also space periodic with period .
Corollary 3 can be deduced using the same steps used to prove Lemma 1 and Corollary 2.
Based on , a continuous operator is defined as
Corollary 4. Let be a space periodic function with period . Then, based on WENO-JS is also space periodic with period .
The proof is straightforward and is based on Corollary 3.
It can be illustrated that is a first-order homogeneous composite function as follows:
The first equality is justified by Equation (73).
The main theorem will be presented through the analysis of the relatively simple WENO-RK1 algorithm.
Theorem 5. Consider the initial value problem of Equation (26), subject to space periodic boundary conditions. Let the initial condition be a space periodic function with period . Then, the numerical solution based on WENO-JS space discretization and RK1 time integration is also space periodic with period .
Proof. The proof proceeds by induction. Let be a space periodic function with period . Applying Equation (23), we get The periodicity of is clarified by Corollary 4, and is periodic by assumption. Hence, the whole right-hand side of Equation (76) is periodic with period . The proof is straightforward for the base case since is periodic with period by assumption.
Corollary 6. Let be a space periodic function with period . Consider the numerical solution based on WENO-JS space discretization and RK1 time integration. If is continuous, then it can be expressed using the following Fourier series:
The coefficients , , , and are the Fourier expansion coefficients defined as follows:
Proof. By definition, Hence, from Theorem 5, Consequently, can be expressed as the infinite trigonometric series provided by Equation (77). The series expression is provided by Equations (3.7, 3.10) in [19], and provided by Equation (77). The series coefficients , , , and are provided by Equations (3.5, 3.6, 3.9) in [19] and given by Equations (78)–(81). The condition of continuity is included to assure the series convergence. For more details, the reader may refer to Theorem 4.3 in [19].
3.3.2. Fourier Analysis of a Monochromatic Solution
Consider a monochromatic with a single period and MWN . As illustrated in Theorem 5, the subsequent solutions will be also exactly periodic with the same period . The variable is defined as a continuous version of index to be
Both variables and are interchangeable. Hence, is defined as
Also, the numerical solution defined in Equation (76) can be written as
Using the continuous presentation, Equation (37) is rewritten in terms of instead of as
Also, Equation (37) is rewritten in terms of as
The periodicity of and are defined as
The periods are defined as
The main analysis step is substituting Equation (88) into Equation (76). Based on Corollary 6, can be expressed as a series expansion of Equation (77). In principle, the coefficients , , , and can be obtained from Equations (78), (79), (80), and (81). However, equivalent formulas can be written in terms of as follows:
Although and refer to the same quantity, the latter is used in Fourier analysis. The specific advantage of as defined in Equation (88) is that is eliminated. The solution behaviour is fully clarified considering only , instead of considering the whole domain .
Although in Equation (87) is monochromatic, the numerical solution includes multiple modes as described in Equation (77). In other words, a monochromatic initial condition induces high-order modes. Assumption 1 will be revisited within the continuous presentation provided. As implied by MI, zero and high order modes can be neglected and Equation (77) can be rewritten as follows:
The spectral parameters are defined based on Fourier expansion coefficients as follows:
Using and , Equation (96) is rewritten as
The parameter describes a shift in the waveform of . This shift does not affect neither nor [19]. In addition, and depend linearly on as illustrated by Equation (75). Consequently, the values and can be adopted without loss of generality. Hence, and can be expressed using Equations (54) and (55), respectively.
The analysis described in Equations (97) and (98) will be denoted by the Nonlinear Spectral Method, abbreviated as NSM. Upon introduction of MI, the focus is limited to calculating only the principal mode spectral parameters, namely, and . Still, it should be emphasized that the error introduced by MI does not affect the accuracy of NSM. In fact the principal mode parameters and are still recovered exactly. The following features of the NSM should be noted: (i)The results of the NSM are valid for any and similar to VN, and the computational effects are eliminated. Hence, the limitations of the ADR developed by [4] are avoided(ii)The expressions resulting from plugging Equation (88) into Equation (76) are mightily lengthy. Hence, this work uses the computer algebra system MAXIMA [20](iii)The integrands of Equations (93) and (94) include rational expressions in terms of trigonometric functions. This is attributed to the form of the nonlinear weights of WENO-JS. Still, the definite integrals can be calculated up to arbitrary accuracy using numerical quadrature. The reader could refer to [21] or other standard textbooks for the details of numerical integration using the trapezoidal approximation(iv)The nonlinear nature of the WENO-JS algorithm should be well captured by the present analysis. The rational expressions of are adopted, unlike [6] where were approximated as constants(v)Past theoretical studies [5, 22] adopted the assumption of small MWN () because Taylor expansion was used. Such limitations are avoided, and the NSM is valid for arbitrary values of
3.3.3. Mode Isolation Error
The error introduced by Assumption 1 will be quantified. Consider the well-known Parseval identity:
It is straightforward to verify that
Hence, the percentage error normalized with respect to is defined as
The error is proportional to the magnitude of zero and higher-order modes truncated by . The reader is reminded that MI is exactly satisfied by LUW5. Hence, a considerable discrepancy may exist between the results of WENO-JS and LUW5 due to , even though the spectra of both methods may be similar in some cases.
3.3.4. Extension to High Order Time Integration
The algorithm described for WENO-RK1 can be applied to WENO-TVDRK2 and WENO-TVDRK3 with slight modifications. The continuous versions of Equations (24) and (25) are defined as follows:
The algorithm described in Section 3.3.2 is modified by replacing Equation (76) with Equations (103) and (104) to analyse WENO-TVDRK2 and WENO-TVDRK3, respectively.
3.4. Geometric Stability Analysis
The spectral methods developed by [5, 6] termed the Geometric Stability Analysis (GSA) will be briefly presented. The spectral properties of the spatial discretization operator are defined in terms of given by
Here, is given by Equation (48) because MI is a cornerstone in GSA [5, 6]. The continuous spectrum is given by
Generally, numerical time integration methods are characterized in terms of given by
The amplification factor for a numerical scheme is given by
The relation between and the spectral parameters in the current work is given by
For any numerical time integration method, the Linear Stability Domain and its boundary are given by
The Geometric Linear Stability Limit (GLSL) is defined as the largest such that
It was illustrated by [5] that the GLSL exists for neither WENO-RK1 nor WENO-RK2, and both were concluded to be linearly unstable. However, the grid given by Equation (2) was used in [6] to define a discrete spectrum that consists of the set of eigenvalues given by
Hence, Equation (111) should be replaced by
It was illustrated by [6] that Equation (113) can be satisfied by both WENO-RK1 and WENO-RK2. Given a specific value of , define the intersection point between and to be . The nonzero eigenvalue of that is closest to the origin is equated to . Hence, the linear stability condition can be given by
For LUW5-RK1, the following formula of that is valid for small MWN was derived by [6]
The function plays a key role in the GSA. However, linearised procedures were used to calculate in [5, 6]. For instance, the GSA of [6] was mainly based on Equation (50) which describes LUW5 instead of WENO-JS. The GSA provided by [5] was mainly based on Taylor expansion and was limited to small . The proposed NSM can be used to calculate exactly as follows. Define and as
Hence, is given by
3.5. Analysis Based on Numerical Experiments
The results provided above can be further verified using numerical experiments. The methodology follows closely the procedure outlined by [4]. Also, complete derivations are provided in Appendix A. For given values of , , and , a grid is adopted based on Equation (2). A value of is selected, where . A monochromatic intermediate condition is calculated based on Equation (87). The solution is obtained by advancing one time step using one of Equations (23), (24), and (25). Finally, Fast Fourier transform (FFT) procedure is applied to to get the spectrum , where . According to the MI assumption, only is considered and other modes are neglected. Hence, the magnification factor and the phase shift are calculated as
The spectral parameters are calculated using the FFT of the solution obtained using a specific grid size . Hence, the procedure will be sensitive to the machine’s precision. A unity domain size was adopted for this paper’s numerical experiments.
4. Results
4.1. Significance of the Nonlinear Analysis
4.1.1. Linear and Nonlinear Weights
To clarify the significance of the NSM, the weights of both LUW5 () and WENO-JS () were plotted. The curves are provided in Figure 2 for four values of , and , covering small and large MWN. A low MWN corresponds to smooth solutions at which the performance of LUW5 and WENO-JS should be close. Hence, for , the difference between and is small. As the value of MWN increases the solution smoothness decreases and the influence of nonlinear WENO-JS weights should be stronger. For instance, at , we get based on Equation (91), and a complete wave is located within 4 grid spacings . Consequently, for large , the difference between and cannot be ignored. This can be easily observed in Figure 2 for .
4.1.2. Mode Isolation Error Results
The error induced by Assumption 1 for WENO-RK1, WENO-TVDRK2, and WENO-TVDRK3 is plotted in Figure 3 for three values of , , and . This range of covers both small and moderate values. Several features can be observed:
The value of vanishes as . This corresponds to the trivial case of a constant , where . (i)On the other hand, also vanishes as . As already illustrated, this corresponds to the maximum sampled by the grid. Hence, the solution is described exactly using a single mode and the MI is exactly satisfied(ii)The value of is almost negligible for the interval . For this range of , modes extend over more than grid points and the solution will be relatively smooth. Hence, WENO-JS accuracy should be relatively higher and MI is well satisfied. This is the opposite for , where modes extend over less than grid points and exhibits strong variation(iii)As increases, increases and MI becomes less satisfied. This is expected since the numerical accuracy is generally inversely proportional to (iv)The value of is generally similar and less than for the three methods WENO-RK1, WENO-TVDRK2, and WENO-TVDRK3 at . However, at , the value of is considerably higher for WENO-RK1 than both WENO-TVDRK2 and WENO-TVDRK3
Based on the above discussion, it is expected that the discrepancy between WENO-JS and LUW5 results should be high for solutions including modes of and using .
4.2. Numerical Solution Spectra
The spectra of WENO-JS are plotted for and , and for the three time integration methods RK1, TVDRK2, and TVDRK3 in Figure 4. Hence, the effects of both and the time integration methods are clarified. For each value of , three spectra are provided based on the following: the NSM explained in Section 3.3.2, FFT numerical experiments explained in Section 3.5, and the exact solution spectrum given by Equations (40) and (41). The results of both FFT and NSM are almost identical, except for some minor discrepancies. Hence, the validity of NSM at the full range of is illustrated. The deviation between the spectra of WENO-JS and the exact solution is proportional to . This inverse proportionality between accuracy and is expected. Most of the spectra of WENO-JS are close to the spectra of the exact solution at . However, of WENO-RK1 is considerably high for at , as shown in Figure 4(a).
Due to the current limited scope, the GSA will be applied to RK1 only. The results of the GSA of WENO-RK1 and LUW5-RK1 are shown in Figure 5. In Figure 5(a), the red unit circle represents of RK1 [13], and the blue and black curves represent of LUW5 and WENO-JS, respectively. Both spectra are calculated exactly using Equations (50) and (117) for LUW5 and WENO-JS, respectively. The full range of the MWN is included in Figure 5(a). For small in the neighbourhood of the imaginary axis, WENO-JS and LUW5 have very close . However, of LUW5 is closer to the imaginary axis than that of WENO-JS. Also, a strong discrepancy can be observed for large away from the imaginary axis.
(a)
(b)
To explain how Equation (113) is satisfied, adopting is plotted for both WENO-JS and LUW5 in Figure 5(b). Similar results can be obtained for other small values of . The region of Figure 5(b) is very close to the imaginary axis and corresponds to low MWN. Two very close values of for LUW5-RK1 are obtained based on Figure 5(b) and Equation (115) to be and , respectively. Hence, the relation with the GSA of [6] is illustrated. Also, it is shown in Figure 5(b) that for WENO-RK1. Hence, the intersection between and is higher for WENO-RK1 than that for LUW5-RK1 at . Also, the grid spacing is restricted to and for WENO-RK1 and LUW5-RK1, respectively. Hence, WENO-RK1 admits a wider range of compared to LUW5-RK1.
In practice, WENO-RK1, LUW5-RK1, WENO-RK2, and LUW5-RK2 require prohibitively small [6]. Hence, the next results will focus on WENO-TVDRK3 because it is considered the most popular combination [23].
The spectra are plotted for WENO-TVDRK3 and LUW5-TVDRK3 using moderate and high values of in Figure 6. Generally, the magnitudes of and of WENO-TVDRK3 are considerably lower than those of LUW5-TVDRK3 for . The low values of and correspond to more dissipation and deceleration of high MWN modes. This is expected because WENO-JS is designed to filter out high MWN modes. Specifically, the weighting algorithm of WENO-JS is focused on the solution behaviour within the five-point stencil. The width of this stencil is , and it is covered by the MWN modes of . Hence, the effect of WENO-JS should be stronger on the MWN modes of . The spectra of WENO-TVDRK3 and LUW5-TVDRK3 show somehow similar trends for . However, the spectra differ considerably at higher . Specifically, two small peaks can be observed for WENO-TVDRK3 at in the range as shown in Figure 6(e). These two peaks merge into a single spike at located at the high MWN as shown in Figure 6(g). This spike contradicts the main requirement of WENO-JS, which is damping high MWN modes. In addition, the effect of this spike will extend to other modes due to the high at (Figure 3). Hence, the convergence of WENO-TVDRK3 solutions at high will deteriorate.
4.3. Numerical Experiments
The significance of the NSM results is clarified using numerical experiments. All the numerical experiments are done using , domain size , grid spacing , and . The solutions are provided at a relatively large time to elucidate any long-term effects. Three experiments are performed adopting the initial conditions given by
The numerical and exact solutions based on Equations (119), (120), and (121) are provided in Figures 7–9, respectively. The results of LUW5-TVDRK3 and WENO-TVDRK3 are provided in the upper and lower parts of Figures 7–9, respectively. The solutions and the magnification factors are provided in the left and right parts of Figures 7–9, respectively. Based on Equation (27), the solutions plotted at correspond to traversing outside and inside the domain and returning to the original position times. Hence, the exact and the initial conditions are coincident and plotted using the same dashed line.
The initial condition of Equation (119) includes a sharp discontinuity and is selected to clarify the performance of WENO-JS for nonsmooth solutions. For the exact solution, approximates a Sinc function with multiple lobes. This behaviour is expected and is explained in [14] and other standard textbooks.
Considering LUW5-TVDRK3, approximates the first three lobes in the range as shown in Figure 7(b). However, almost vanishes for high MWN in the range and the spectrum is almost truncated. This truncation leads to an oscillatory behaviour and Gibbs oscillations can be observed in the solution plotted in Figure 7(a). This truncated spectrum is expected since was close to in the range for LUW5-TVDRK3 spectra presented in Figure 6. The behaviours of LUW5-TVDRK3 solutions are almost similar for .
For WENO-TVDRK3, the solution behaviour is strongly affected by . Considering low and moderate CN (), multiple lobes are included in and cover the whole MWN range as shown in Figure 7(d). This distribution of approximates the product of two Sinc functions that describes a relatively smeared but non-oscillatory solution. Hence, WENO-TVDRK3 solutions at approximate the exact solution avoiding Gibbs oscillations as shown in Figure 7(c). Inverse proportionality between and is observed in Figure 7(c). However, does not vanish at high MWN in the range . This behaviour is attributed to the nonlinear nature of WENO-JS, where low MWN modes contribute to high MWN modes as illustrated in Section 3.3.2 and quantified by .
The distribution of for WENO-TVDRK3 considerably deteriorates at the high CN (). The Sinc distribution disappears and a strong peak is observed at as shown in Figure 7(d). This peak is expected and corresponds to the spike that is clarified in Figure 6(g). As already mentioned in Section 4.2, this spike has an adverse effect on the solution. The spike tends to magnify high MWN modes leading to nonsmooth and irregular solutions. On the other hand, the WENO-JS mechanism tends to smooth the solution and remove the irregularity. As a result, the solution of WENO-TVDRK3 for shown in Figure 7(c) is overdamped and infested by small yet distinguishable spurious oscillations.
The detection of this overdamping phenomenon is generally difficult using numerical experiments. For instance, consider the solutions based on Equation (120) plotted in Figure 8. This initial condition, given by Equation (120), includes only one low MWN mode (). Both WENO-TVDRK3 and UW5-TVDRK3 solutions are almost indistinguishable, and the effect of disappears. The reason for this is that overdamping occurs due to the influence of high MWN modes. A numerical experiment was done by [5] at using a low MWN initial condition. Hence, overdamping was hidden, and it was erroneously concluded that the performance of WENO-TVDRK3 is generally satisfactory at this high CN.
Finally, the sensitivity of WENO-TVDRK3 to high wave numbers is illustrated by considering the initial condition of Equation (121). This condition consists of a high amplitude low MWN mode (), and a low amplitude high MWN mode (). The amplitude of the high MWN mode () is almost negligible compared to that of the low MWN mode. The influence of mode is negligible for UW5-TVDRK3 solution plotted in Figure 9(a), and it is effectively filtered in UW5-TVDRK3 spectra plotted in Figure 9(b). However, the solution of WENO-TVDRK3 suffers overdamping at as illustrated in Figure 9(c). A small yet distinguishable peak located at can be observed in the WENO-TVDRK3 spectrum for in Figure 9(d). This peak corresponds to the spike that is clarified in Figure 6(g).
5. Conclusions
The nonlinear Weighted Essentially Nonoscillatory method developed by Jiang and Shu (WENO-JS) is analysed to clarify its spectral properties. An exact Nonlinear Spectral Method (NSM) is developed based on the Fourier expansion of a continuous presentation of WENO-JS solutions. The nonlinear rational formula of WENO-JS weights is fully included. Unlike other theoretical methods, the NSM is valid for arbitrary modified wave numbers (MWN). The mode isolation (MI) assumption, which was used extensively in the literature, is quantified for the first time, and its impact on the performance of WENO-JS is explained. The MI error increases for high MWN and moderate Courant numbers (CN), and the discrepancy between the results of WENO-JS and those of the fifth-order Linear Upwind (LUW5) discretization cannot be neglected.
The NSM can reproduce the spectra of numerical experiments while being independent of any grid parameters. The effects of both the CN and the numerical time integration methods are elucidated. At low CN, the performance of the combination of WENO-JS with the forward Euler time integration method is better than that of LUW5. At high CN, the combination of WENO-JS with the third-order variation diminishing Runge-Kutta method suffers extra damping compared to LUW5. This peculiar phenomenon was overlooked in past studies because it occurs due to a spike located at the high MWN spectra, and its detection via numerical experiments is difficult. The results of the current study should be extended to include other time integration algorithms and nonlinear discretizations.
Appendix
A Discrete Fourier Transform
Consider the periodic function of period . Samples of are taken at points , where . The samples are obtained at equally spaced points .
Hence, the sequence is obtained. The sequence can be written as the following sum [14]:
Here, . The coefficients are termed as the discrete Fourier transform or DFT as an abbreviation. The operator is used to define DFT. This is defined as
Similarly, the inverse operation of DFT can be defined using the operator notation:
For a better presentation, the wave number is defined as
And the modified wave number (MWN) is defined as
Hence, Equations (A.2) and (A.3) can be rewritten as
B DFT of a Monochromatic Wave
Consider the real monochromatic data series, with an amplitude , a wave number , and an arbitrary phase shift .
The discrete Fourier transform is obtained as
Consider the case :
Here, . The terms are constants w.r.t. the summation. Using the geometric series sum formula , it is concluded that since . Hence,
Consider the case :
Similarly, is defined as . Following the same results of the geometric series, it is concluded that vanishes. Hence,
Following the same steps, it can be easily deduced that that for . These results can be verified by substitution by Equations (A.12) and (A.14) into Equation (A.7).
Finally, referring to Equations (A.12) and (A.14), it is deduced that
The conclusions are summarized in the following points: (i)The amplitude of mode is recovered from the DFT coefficient as (ii)The phase shift of mode is recovered from the DFT coefficient as (iii)The range of as defined in Equation (A.6) is . However, the effective range of is reduced to due to the property of Hermitian Symmetry, as illustrated by Equations (A.16) and (A.17)(iv)The results above were derived for a monochromatic signal. However, due to the linearity of DFT [14], the results still hold for any periodic signal. Consider the general case:
The DFT can be obtained as shown in Equation (A.3). Hence, the amplitude and phase shift of mode can be recovered back from the DFT as follows:
As mentioned in [6], can be assumed even without loss of generality. Based on the above discussion, a discrete mesh yields only a finite number of distinct modes. These modes can be represented using , or where .
Data Availability
Data will be provided upon reasonable request.
Conflicts of Interest
The author declares no conflicts of interest. The research was performed as a part of the employment of the author at Nile University.
Acknowledgments
The author is grateful for the administrative support and encouragement provided by Nile University. Also, the author would like to thank his mentors and colleagues for their suggestions during the manuscript preparation.