Stability and Resolution Analysis of the Wavelet Collocation Upwind Schemes for Hyperbolic Conservation Laws

1. Introduction

High-speed flows governed by hyperbolic conservation laws often contain steep gradient regions or even discontinuities, such as shock waves and contact discontinuities. Pirozzoli [1] pointed out that ideal numerical methods for such problems should be of high accuracy and free from numerical dissipation in smooth regions of the flows, and they can capture steep gradient regions and discontinuities robustly without significant spurious oscillations. The numerical methods for high-speed flows can be classified into two classes, one dealing with smooth flows and the other with shock waves. The central finite difference schemes with null dissipation error and spectral-like resolution are the ideal candidates for computations in smooth regions [2]. However, it is well known that the central discretization for smooth flows will induce severe spurious oscillations and lead to numerical instability in the presence of the discontinuities, where extra techniques are required to distinguish shocks [1]. Although artificial viscosity added near the discontinuities may help to suppress oscillations and stabilize calculations, it is quite difficult to introduce the corresponding appropriate numerical viscosity to suppress these oscillations successfully [3]. Hence, a preferable choice is to design schemes that can achieve stable solutions in the inviscid limit.

The upwind method based on the philosophy that a stable numerical method should propagate the information along the directions of characteristic waves have been proven as a good choice [4,5]. The upwind schemes can introduce the appropriate implicit numerical viscosity to ensure the stability for solving hyperbolic conservation laws [6]. High-order upwind schemes have been widely investigated on numerical simulation of high-speed flows during last decades owing to their capability in obtaining satisfactory resolution and high-order accuracy [7]. The typical examples of these schemes include essentially non-oscillatory (ENO) schemes [8], weighted ENO (WENO) schemes [9], multimoment constrained finite volume (MCV) methods [10], discontinuous Galerkin (DG) schemes [11], and dynamical wavelet Galerkin schemes [12], etc.

Multiresolution analysis in wavelet theory provides an effective approach to capture the localized structures of functions both in physical and frequency domains, which is suitable for simulations of high-speed flows. Many wavelet-based numerical methods have been developed in computational fluid dynamics [13], which can be categorized into two types. The first one only embeds the wavelet multiresolution analysis into the traditional high-order schemes. The remainder is called a pure wavelet numerical method that uses the wavelets as a set of bases to discretize the PDEs. Vasilyev and the co-authors have conducted much pioneering work in wavelet-based numerical methods for computational fluid dynamics, which have been used successfully in solving Burgers Equation with viscosity [14], laminar flame–vortex interaction [15], homogeneous turbulence problems [16] and supersonic channel flow [17].

However, comparatively few works are dedicated for hyperbolic partial differential equations (PDEs). Wavelet methods coupled with the DG schemes [18], the finite difference WENO [19], and the FDM with artificial viscosity [6] are devised for the hyperbolic PDEs. These methods use the wavelets to detect the localized steep structures and implement the adaptive or reconstruction process, but omit the merit of adaptive wavelet approximation on arbitrary finite domains. In the pure wavelet method framework, Restrepo and Leaf [20] firstly applied the classic wavelet Galerkin schemes with uniform nodes to solve hyperbolic equations and concluded that spurious oscillations would spread further away from the shock, leading to numerical instability, which was also encountered in other classic Galerkin schemes [21]. Recently, a dynamical wavelet Galerkin scheme was proposed by Pereira et al. [12] that overcomes the above drawback successfully by the energy dissipation introduced by a non-smooth projection operator. The collocation methods are more efficient compared with the wavelet Galerkin ones. Although adaptive wavelet collocation methods are developed for solving parabolic problems successfully based on symmetrical interpolation wavelet basis [14,22,23], no attempts are carried out to use the wavelet collocation method to construct the upwind schemes for the hyperbolic problems. The effects of the scaling function characteristics on the wavelet schemes have also been unexposed until now.

Motivated by the idea of the upwind schemes, we have considered building asymmetrical interpolation wavelets to achieve the upwind property, and proposed high-order wavelet collocation upwind schemes based on the asymmetrical interpolation wavelets and corresponding function approximation theory in our recent work [24]. The numerical results show that the order of accuracy of the wavelet schemes agrees with the expected ones. Unfortunately, this only provides information on the asymptotic convergence rate to the exact solution, and nothing is known about the stability and resolution of the schemes. Dissipation and dispersion analysis can reveal more details of the schemes. Taylor expansion analysis [25] and Fourier analysis [26] are the two main methods in evaluating the dissipation and dispersion properties of the numerical schemes. For the former one, Taylor expansion of the equation is conducted to obtain the modified equation which represents the actual partial differential equation solved by numerical schemes, and a truncated version of the modified equation are applied to gain both dissipative and dispersive errors [25]. Fourier analysis method is on the basis of Fourier series and its derivative operator to examine the dissipative and dispersive properties of the equation, which can be implemented more conveniently and applied in more situations [26]. In addition, a numerical analysis method is proposed by Pirozzoli [27] to evaluate the dissipation and dispersion performance for the nonlinear high-order schemes for hyperbolic conservation laws.

The goal of the present paper is to explore the effects of characteristics of the scaling functions on wavelet collocation upwind schemes by conducting numerical tests and the Fourier analysis on several wavelet upwind schemes. We obtain the dissipation and dispersion properties, and analyze the stability and resolution properties of these schemes. Then we provide a fundamental rule to construct stable, high-order and high-resolution schemes and verify the outcome by conducting two typical numerical tests.

The remainder of the paper is organized as follows. In section 2, we briefly introduce the wavelet approximation theory and wavelet collocation upwind schemes for uniform node distribution, and define the parameters to describe the asymmetry of the scaling function. Numerical experiments for the linear scalar equation and Fourier analysis on the schemes are carried out in Section 3. Finally, the main conclusion is summarized in section 4.

2. Wavelet collocation upwind schemes

In the present work, we consider one-dimensional scalar conservation laws to analyze the stability of the wavelet collocation upwind schemes:

$$u_t+f{\left(u\right)}_x=0$$

(1)

2.1. Wavelet approximation theory

2.1.1. Preliminaries

The Hilbert space $L^2(\Omega )$ is defined as the collection of square integrable functions that satisfy ${{\displaystyle {\int }_{\Omega }\left|f(x)\right|}}^2\mathrm{d}x<\infty$, where Ω is any open subset of R [28]. The space is equipped with the inner product:

$$\left(f,g\right):={\displaystyle {\int }_{\Omega }f(x)g(x) }\mathrm{d}x<\infty \hspace{1em}\hspace{1em}\mathrm{for} \mathrm{all} f,g\in L^2(\Omega ).$$

(2)

The corresponding $L^2(\Omega )$-norm is defined by

$${‖f‖}_{L^2\left(\mathsf{\Omega }\right)}={\left({{\displaystyle {\int }_{\mathsf{\Omega }}\left|f(x)\right|}}^2\mathrm{d}x\right)}^{1/2}\hspace{1em}\hspace{1em}\mathrm{for} \mathrm{all} f\in L^2(\Omega ).$$

(3)

The definition of $L^{\infty }(\Omega )$-norm is expressed as

$${‖f‖}_{L^{\infty }\left(\mathsf{\Omega }\right)}=\underset{x\in \Omega }{\sup }\{f(x)\}.$$

(4)

For any integer m ≥ 1, the space of all functions which are m times continuously differentiable over Ω is denoted by ${\mathrm{C}}^m(\mathsf{\Omega })$.

2.1.2. Approximation of functions on a finite domain

For general practical problems, Ω is a finite domain. We apply the boundary extension technique in our previous study based on the Lagrange interpolation to remove local errors induced by a loss of information outside the domain [29,30]. For all functions in $L^2(\Omega )$, the modified wavelet approximation at the resolution level J can be written as

$$P_{J}f(x)={\displaystyle \sum _{k\in {ℝ}_J}f(x_k){\mathsf{\Phi }}_{J,k}(x)},$$

(5)

in which the modified wavelet basis is given by

$$\begin{array}{ll}{\mathsf{\Phi }}_{J,k}(x)& ={\phi }_{J,k}(x)+{\displaystyle \sum _{n\in {\overline{\mathsf{\Xi }}}_k}{\mathsf{\Pi }}_k^n\left(x_n\right){\phi }_{J,n}(x)}\\ & ={\displaystyle \sum _{n\in {\mathsf{\Xi }}_k}{\mathsf{\Pi }}_k^n\left(x_n\right){\phi }_{J,n}(x)},\end{array}$$

(6)

$${\mathsf{\Pi }}_k^n\left(x_n\right)={\displaystyle \prod _{\begin{array}{l}i=1\\ i\ne n\end{array}}^{\eta }\frac{x_n-{\tilde{x}}_i}{{\tilde{x}}_k-{\tilde{x}}_i}},$$

(7)

where ${\overline{\mathsf{\Xi }}}_k$ is the set of serial number n denoting the external point x_n such that f (x_n) is exactly all the Lagrange interpolations. In Equation (7) the set {${\tilde{x}}_i$, i = 1, 2, …, η } is the collection of selected Lagrange interpolation nodes inside the Ω for the external point x_n and ${\mathsf{\Xi }}_k={\overline{\mathsf{\Xi }}}_k\cup \left\{k\right\}$ with the coefficient ${\mathsf{\Pi }}_k^k\left(x_k\right)\equiv 1$. Next we give two wavelet approximation theorems without proof which can be proved by using the approach as presented in the references [29,31].

Theorem 2.1

For all f ∈$L^2(\mathsf{\Omega })\cap {\mathrm{C}}^m(\mathsf{\Omega })$, suppose that m is large enough and let$\eta =min\left\{{\eta }_k,k\in {\overline{\mathsf{\Xi }}}_k\right\}$, approximation errors in$L^2$-norm and$L^{\infty }$-norm can be estimated by

$${‖P_{J}f-f‖}_{L^2\left(\mathsf{\Omega }\right)}\le C_{1,L}{2}^{-J\lambda }$$

(8)

$${‖P_{J}f-f‖}_{L^{\infty }\left(\mathsf{\Omega }\right)}\le C_{1,\infty }{2}^{-J\lambda }$$

(9)

where λ = min {N,η}. Constants C_1,L and C_1,∞ are dependent on the regularity of the derivatives f ⁽ⁿ⁾(x) but independent of the resolution level J [29].

Theorem 2.2

For all f ∈${\mathrm{L}}^2(\Omega )\cap {\mathrm{C}}^m(\Omega )$, suppose that m is large enough and let$\eta =min\left\{{\eta }_k,k\in {\overline{\mathsf{\Xi }}}_k\right\}$, approximation errors of the first order derivative in${\mathrm{L}}^2$-norm and${\mathrm{L}}^{\infty }$-norm can be estimated by

$${‖\frac{\mathrm{d}{P}_{J}f}{\mathrm{d}x}-\frac{\mathrm{d}f}{\mathrm{d}x}‖}_{L^2\left(\mathsf{\Omega }\right)}\le C_{2,L}{2}^{-J(\lambda -1)},$$

(10)

$${‖\frac{\mathrm{d}{P}_{J}f}{\mathrm{d}x}-\frac{\mathrm{d}f}{\mathrm{d}x}‖}_{L^{\infty }\left(\mathsf{\Omega }\right)}\le C_{2,\infty }{2}^{-J(\lambda -1)},$$

(11)

where λ = min {N,η}. Constants C_{2, L} and C_{2, ∞} are dependent on the regularity of the derivatives f ⁽ⁿ⁾(x) but independent of the resolution level J [29].

2.2. Wavelet collocation upwind scheme

Inspired by the concept of the upwind scheme, we have proposed pure adaptive wavelet collocation upwind schemes to resolve hyperbolic conservation laws in our recent research, which have been used in solving 1D hyperbolic conservation laws successfully [24]. The upwind property is achieved by building a couple of asymmetrical wavelets. We can calculate the derivative based on the wavelet approximation:

$$u^{\prime }(x_{J,l})={\displaystyle \sum _{k\in {ℝ}_J}{u}_{J,k}}{{\mathsf{\Phi }}^{\prime }}_{J,k}(x_{J,l}).$$

(12)

Owing to the interpolation property of the scaling functions, we can discretize the nonlinear terms in the conservative equations directly [32]:

$$f\left(u(x)\right)={\displaystyle \sum _{k\in {ℝ}_J}f\left(u_{J,k}\right)}{\mathsf{\Phi }}_{J,k}(x),$$

(13)

where f satisfies a uniform Lipschitz condition of order α with respect to u, α ≥ 1. We observe that decomposition coefficients of $f\left(u\right)$ can be evaluated by computing value of $f\left(u_{J,k}\right)$ in a very high efficiency.

To concentrate on the stability and resolution analysis of the wavelet schemes, only the linear scalar equation with a periodic boundary condition and$f^{\prime }(u)\ge 0$ is considered. We conduct spatial discretization at the resolution level J by the wavelet collocation upwind schemes in the uniform node distribution framework and obtain the following semi-discretized system:

$$\frac{\mathrm{d}{u}_J(x_l,t)}{\mathrm{d}t}+{\displaystyle \sum _{k\in {ℝ}_J}f\left(u_J(x_k,t)\right)}{{\mathsf{\Phi }}^{\prime }}_{J,k,+}(x_l)=0,\hspace{1em}\hspace{1em}k,l\in \Omega ,$$

(14)

where ${\mathsf{\Phi }}_{J,k,+}(x)$ denotes the scaling function of the positive upwind wavelet. Then we can solve the above ordinary partial equations by the classic explicit fourth-order Runge–Kutta method for time integration.

2.3. Asymmtrical wavelets

The detailed procedures of constructing the asymmetrical wavelets can be found in the Appendix A. We establish a number of asymmetrical wavelets from N = 3 to N = 10 to examine the stability and resolution of these wavelet schemes. Here N describes the regularity of the scaling function and means that the corresponding scaling function can reproduce polynomials up to the degree (N – 1).

For convenience of discussion, we define two parameters to depict the asymmetrical properties. We take the scaling functions with N = 6 and N = 7 as examples. The stencils for the scaling functions are shown in Figure 1. The bias magnitude of the stencil is defined by the difference between the number of nodes on the left side of $x_{J+1,l}$ and that on the right side of $x_{J+1,l}$ as denoted by the BM. To evaluate the symmetry of the scaling function, we also use symmetry factor (SF) to evaluate the symmetry of the scaling function which can be calculated by

$$\mathrm{SF}=\frac{\mathrm{IL}}{\mathrm{IR}},$$

(15)

where IL is the length of the support interval on the left side of the zero point, and IR is the length of the support interval on the right side of the zero point. The scaling function is exactly symmetrical when SF = 1. We can provide an explicit formula to compute the SF base on the BM as follows:

$$\mathrm{SF}=\frac{N-\left(\mathrm{BM}+1\right)}{N+\left(\mathrm{BM}-1\right)}.$$

(16)

It can be observed from Equation (16) that the SF approaches to 1 as N tends to infinity, which reveals that the symmetry of the scaling function is improved as N increases. For a specified N, the SF decreases with an increment in BM. The BM and SF of the scaling functions used in the present paper for analyzing the stability and resolution of the wavelet upwind schemes are listed in Table 1. We note that the BM has the same parity with the N. The SF with the N_odd is greater than that with the N_even when the BM_odd and N_odd are both smaller than the BM_even and N_even by one, respectively. Here the subscript represents the parity of the N. Some scaling functions are shown in Figure 2 as examples.

3. Stability and resolution analysis of the wavelet upwind scheme

In this section, we conduct two benchmark tests by applying different schemes from N = 3 to N = 10. Two kinds of norms are defined to evaluate numerical errors as follows:

$${‖e‖}_{l_{\infty }}=\underset{k}{\max }\left\{ \left|u_k^e-u_k^n\right| \right\},$$

(17)

$${‖e‖}_{l_2}={\left({\displaystyle \sum _k{\left|u_k^e-u_k^n\right|}^2\mathsf{\Delta }x}\right)}^{\frac{1}{2}},$$

(18)

where $u_k^e$ is the exact solution and $u_k^n$ is the numerical one.

Numerical examples are performed to the following one-dimensional linear scalar equation:

$$u_t+a{u}_x=0,x\in \left[-1,1\right],$$

(19)

where a = 1.

3.1. Advection of a sine wave

At first, we verify the order of accuracy and explore the stability of the proposed schemes by refinement tests. A sine wave as the initial condition is presented as follows:

$$u(x,0)=\sin (\pi x)\hspace{1em}\hspace{1em}\mathrm{periodic}.$$

(20)

The periodic boundary condition is addressed based on the extension method for periodic functions proposed in our previous work [33].

The numerical errors are calculated at t = 2. BM_odd = 1 and BM_even = 2 are taken for the different schemes. Time steps are set small enough to remove the influence of time integration. Numerical errors and orders of accuracy for all schemes are shown in Table 2. It should be noted that the errors which are smaller than 5.0E-13 are removed from the table since they approach to the machine precision of the double type. Theorem 2.2 show that theoretical order of accuracy of the scheme is (N − 1) order. It can be seen that the orders of accuracy are consistent with the expected ones for the schemes that N are smaller than 8. Higher order schemes give the expected order of accuracy when the number of nodes is small. The N = 4 scheme is unstable even for this smooth wave evolution in our tests. Therefore the result is absent from Table 2.

To explore the effect of asymmetry of the scaling function on stability, we also compute the numerical results by respectively using the wavelet schemes with N = 7, N = 9, and different BM values as shown in Table 3. It can be observed that the scheme with N = 7 and BM = 3 is unstable because the error increases significantly as J goes up to 8. We conclude that the scheme with the same N will be unstable when the asymmetry of the scaling function increases.

3.2. Advection of a square wave

To validate the stability of the schemes to solve a jump discontinuity problem, we conduct the advection of a square wave which has the following initial distribution:

$$u(x,0)=\left\{\begin{array}{l}1\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}-0.4\le x\le 0.4,\\ 0\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\mathrm{otherwise},\end{array}\right.\hspace{1em}\hspace{1em}\mathrm{periodic}.$$

(21)

The Fourier series of the square wave is the infinite superposition of sine and cosine waves, and contains all frequency components. Numerical results at t = 8 obtained by different wavelet schemes at the resolution level J = 8 are illustrated in Figure 3. It can be seen that spurious oscillations pollute numerical solutions near the jump discontinuities for almost all schemes. Figure 3b shows that the oscillations obtained by the scheme with N = 6 is gradually amplified in the positive direction, which indicates that the scheme is unstable for discontinuity problems. For all the schemes with specified parity of N, the length of the oscillation interval decreases with an increment in N. We continue to compute the results at t = 32 as shown in Figure 4. The results for the scheme with N = 6 is divergent for long-term time integration, which is not depicted in the figure. Although this scheme is stable for the single sine wave advection, it behaves unstable for that of the waves with multiple or infinite frequencies. It can be observed that the schemes with N∈odd confine spurious oscillations in the thinner region, and the amplitude of the oscillations attenuates rapidly. In addition, N = 7, BM = 3 and N = 9, BM = 3 schemes are unstable for the current problem.

3.3. Dissipation and dispersion analysis

Although we can devise wavelet upwind schemes based on our proposed method, the above numerical tests suggest that some of the wavelet schemes are unstable even for the linear scalar advection problems. Therefore, it is crucial to investigate the dissipation and dispersion properties and provide an instruction to choose “good” wavelets for stable wavelet upwind schemes.

We examine the dissipation and dispersion properties of the schemes base on the Fourier analysis method [2, 26]. Suppose that $u(x_l)=e^{ik{x}_l}$, then we can obtain the numerical derivative:

$${\delta }_x{u}_l=\frac{\tilde{k}}{\mathsf{\Delta }x}{e}^{ik{x}_l},$$

(22)

where $\tilde{k}$ is the corrected wavenumber which is a complex variable and $\mathsf{\Delta }x=1/2^J$. Therefore, $\tilde{k}$ can be rewritten as $\tilde{k}=k_r+i{k}_i$. The exact derivative is given by

$$\frac{\mathrm{d}}{\mathrm{d}x}{u}_l=ik{e}^{ik{x}_l}.$$

(23)

Comparing the above relation, we can obtain that k_r = 0 and k_i = kΔx for ideal numerical schemes. When conducting Fourier analysis on Equation (19), we can calculate the following numerical solution:

$$u\left(x_l,t\right)=e^{-k_{r}at/\mathsf{\Delta }x}{e}^{ik\left(x_l-at{k}_i/\left(k\mathsf{\Delta }x\right)\right)},$$

(24)

where k_r is the dissipation coefficient and k_i is the dispersion coefficient. The exact solution for Equation (19) can be computed by

$$u\left(x_l,t\right)=e^{ik(x_l-at)}$$

(25)

Then we can give the physical meaning of k_r and k_i. k_r depicts the attenuation of the wave amplitude at t induced by the numerical error, and k_i describes the change in propagation speed of the wave caused by the numerical method. It can be noted that the wave amplitude will be amplified over time if k_r is a negative value. The corresponding scheme shows the instability when solving the advection problems. Therefore, k_r for stable numerical schemes should be non-negative. Here kΔx is defined as effective wavenumber and denoted by α. When k_i / α > 1, the numerical wave will run faster than the exact one. And the numerical wave will fall behind the exact wave for k_i / α < 1. We remark that only the spectral method has the ideal dispersion property, which means that k_i /α = 1.

We compuate the dissipation coefficient of the different wavelet schemes as shown in Figure 5. For specified parity of N, it can be seen that k_r decreases with the increase of N. k_r of the scheme with N∈odd is smaller than that with N∈even, indicating that the wavelet upwind schemes with N∈odd show better dissipation property. To seek the reason for the instability of several schemes, we analyze the k_r for all the above schemes. We find that the negative k_r exists for schemes with N∈even. To clarify this fact, the locally enlarged version of the Figure 5b is illustrated in Figure 6. It can be observed that k_r are all negative when α < 0.8. This suggests that the schemes with N∈even will show a negative diffusion phenomenon when J is large enough, and the error will be amplified over time. This actually induces the instability of the wavelet schemes. We also find that the local extreme value approaches to 0 as N increases for N∈even, which means that the negative diffusion process is weakened, and the stability of the schemes are improved. This explains that the numerical tests for the schemes with N∈even in subsection 3.1 and 3.2 are still stable. But for longer time integration and some α near the local extreme value, the results might be divergent.

Based on the above analysis, we can conclude that only N∈odd schemes provide the correct implicit viscosity for solving the hyperbolic conservation laws. However, numercial tests in subsection 3.1 show that the scheme with N = 7 and BM = 3 is still unstable. As has been discussed, the asymmetry of the scaling function is another factor influencing the stability of the wavelet schemes. We further evaluate the dissipation coefficients of the schemes with N = 7, BM = 3 and N = 9, BM = 3 as shown in Figure 7. It can be found that k_r is negative when α is smaller than the specified value. Therefore, the schemes with BM = 3 are unstable. Now we can achieve a basic instruction that N∈odd, BM = 1 schemes are with non-negative dissipation coefficients and stable for hyperbolic conservation laws.

Then, we will explore the dispersion property related to the resolution of the schemes. For a specified wavelet scheme, an effective node number is defined by the point per wavelength abbreviated as PPW, which can be computed by the following relation

$$PPW=\frac{\lambda }{\mathsf{\Delta }x}=\frac{2\pi }{k\mathsf{\Delta }x}=\frac{2\pi }{\alpha }.$$

(26)

The spectral method has the optimal resolution corresponding to α = π and PPW = 2. It can be seen from Equation (26) that PPW is inversely proportional to α. The length of the interval α that k_i / α is approximately equal to 1 depicts the ability of the scheme to trace the wave accurately. We choose the interval of α that satisfies $\left|1-k_i/\alpha \right|<5\%$ and $\left|1-k_i/\alpha \right|<2\%$ to measure the maximum α which reflects the resolution of the schemes directly. The dispersion coefficients against α are plotted in Figure 8. It can be observed that N < 5 schemes have a large dispersion error and a low resolution. The maximum α that k_i / α meets the tolerance relation increases with an increment in N for the specified parity. To clarify the resolution more clearly, we list the maximum α in a tolerance range in Table 4. It can be seen that the maximum α gradually tends to π as N increases, and the schemes with N∈even behave better in resolution when N > 6. On the basis of the above analysis, we can obtain that the wavelets are more applicable to design the high-order schemes, and the scheme with larger N has the higher accuracy and better resolution.

Next, we conduct a numerical test with a smooth initial distribution to verify the theoretical analysis of the dissipation and dispersion performance. The test is the advection of a sine wave described in subsection 3.1 with u(x, 0) = sin (5πx) as the initial condition. The PPW is approximately equal to 6 at the resolution level J = 4. The numerical results obtained by the schemes with N∈odd and BM = 1 at t = 2 are shown in Figure 9. It can be seen that the higher order schemes show less dissipative, approach to the exact solution more accurately and reveal higher resolution for the wave. Therefore, the numerical results are in accordance with the stability and resolution analysis.

Finally, a numerical test of the advection of a multi-scale function is devised to show the capability of the wavelet schemes in recognizing smooth and discontinuous solutions in the high-speed flows. The initial condition consists of step functions, saw-tooth function, sine waves in different frequencies, which is shown as

$$u(x,0)=\left\{\begin{array}{ll}0.5& -0.8\le x\le -0.6,\\ -1.25x+0.25& -0.6\le x\le -0.2,\\ 0.5& -0.2\le x\le -0.1,\\ 0.5(1+\sin (40\pi x))& -0.1\le x\le 0.1,\\ 0.5& 0.1\le x\le 0.2,\\ 0.5+0.2\sin (5\pi (x-0.2))& 0.2\le x\le 0.6,\\ 0.5& 0.6\le x\le 0.8,\\ 0& \mathrm{otherwise}\end{array}\right.$$

(27)

For this numerical test, we apply the adaptive wavelet upwind schemes proposed in our previous study [24]. The main idea of adaptive node generation is to recognize trouble nodes based on the wavelet coefficients that are larger than a threshold parameter ε = 1.0⁻⁵, and insert nodes in the adjacent zones near the trouble nodes. Moreover, an integration reconstruction method is designed based on the Lebesgue differentiation theorem to suppress the spurious oscillations. We choose the basic resolution level J₀ = 6, the maximum resolution level J_max = 12, the same adaptive and reconstruction parameters for different schemes. The numerical results obtained by the adaptive wavelet schemes with N∈odd and BM = 1 at t = 2 are compared with that of the classic fifth-order finite difference WENO scheme (WENO-5) proposed by Jiang and Shu [34] as illustrated in Figure 10. It can be found that the higher order wavelet scheme can capture the discontinuities without spurious oscillations and distinguish different scale structures accurately with less nodes, which also verifies the better resolution of the higher order scheme. For the WENO-5 scheme, a uniform node distribution with N₁ = 2048 is required to depict all the details of the solution. The nodes required in the adaptive wavelet upwind schemes are about half of the WENO-5 scheme, showing that the adaptive wavelet schemes with the integration reconstruction can capture discontinuities free from the numerical oscillations and distinguish complex solutions efficiently.

5. Conclusions

In the present paper, we construct a system of wavelet collocation upwind schemes and conduct linear advection tests and Fourier analysis to uncover the effects of the characteristics of scaling functions on stability and resolution of the constructed schemes. The convergence rates of the schemes are consistent with the expectation. The dissipation analysis suggests that one can apply scaling functions of the wavelets with N∈odd, BM = 1 to construct the high-order and stable upwind schemes for hyperbolic conservation laws. The higher order scheme has a better resolution which tends to the spectral resolution as N increases. Two typical numerical tests verify that the constructed wavelet collocation upwind schemes have the desired properties and can be used to solve high-speed flows with multi-scale smooth structures and discontinuities in an efficient way.