It is well known that lower order numerical methods, such as various monotone methods, behave well near discontinuities, while high order numerical methods, such as the Lax-Wendroff method, work well in smooth regions (see [14, 15] for example). Monotone methods are first-order accurate due to Godunov's theorem. They do not produce non-physical phenomena such as smearing of the solution or spurious oscillations. On the other hand, high-order methods yield non-physical phenomena when shocks are presented.

The total variation diminishing (TVD) technique has long been used as one of the various techniques for avoiding the spurious or non-physical oscillations exhibited by high-order schemes. It requires one to modify a high-order method such as the Lax-Wendroff method (see [2, 13, 15, 19] and etc., for example) using flux- or slope- limiter techniques.

The main goal of this paper is to introduce a correction technique to the upwind method, another known method for eliminating undesired non-physical phenomena, to solve a simple scalar linear hyperbolic equation without the help of any limiter techniques. The first step towards realising this goal is to split the hyperbolic equation u_t + a u_x=0 by introducing a new flux variable v= a u_x (see [10]). As a result, we have a system of ordinary differential equations consisting of v= a u_x and u_t + v =0. Our conservative method then consists of two-steps: the first-step is to set up a correction scheme for v=au_x and the second-step is a scheme for u_t+v=0. Once we deal with the two equations separately and then combine them, we have a new conservative scheme which is composed of the upwind scheme plus a correction term.

The correction term behaves like a high-resolution correction (see p.163 in [15]) while the upwind method of the form Ujm−ν2(Ujm−Uj−1m) obtains a low-resolution behavior. As a result, we obtain a stable scheme which is second-order accurate for u. It has been also been frequently observed [4, 11, 12, 13, 14, 15, 17, 18, 19] that a high resolution method using a flux-limiter has a restricted TVD region, and so the limiter functions must be chosen within this TVD region. Our proposed method does not use any limiter functions, so this is not a concern. Moreover it is not only second-order accurate, but also has the almost l₁ contracting property.

Our other goal is to computationally compare the proposed algorithm with the well-known upwind, Lax-Wendroff and flux-limiter methods, and to show the second-order local truncation error of the proposed algorithm. The order of convergence of the upwind correction scheme in comparible to that of min-mode type schemes reported in [8] and [16] and that of the fully-discrete high-resolution schemes with van Leer's flux limiter reported in [7], for example.

This paper is layed out as follows. In Section 2, we will present how the algorithm can be derived. In Section 3, the almost l₁ contracting property and truncation error estimates are shown. In Section 4, numerous examples are presented by comparing the present method with other well-known methods using flux-limiters. In the concluding section, we mention and discuss a possibility for expanding the present method to a nonlinear scalar hyperbolic equations.

We consider the following simple linear hyperbolic equation

with initial data

By introducing a new variable v= a u_x, (2.1) can be written as

Lemma 2.1. Assume that the solution of (2.1) is sufficiently smooth. Then the system of ordinary differential equations (2.3) is equivalent to

Proof. From (2.3), it follows that

which implies (2.4). On the other hand, assuming (2.4) we have

Hence, we have v - a u_x = c. Then, the initial condition shows v = a u_x. This implies (2.3).

Now, let us discuss the discretization for which we assume that the domain (−∞,∞)×[0,T] has the discrete mesh points (x_j, t_m) given by

where h and τ denote a mesh width and a time step, respectively.

We will present a conservative method for u_t = - v by applying 3^rd order Runge-Kutta method which leads us

where, for an approximation of v(x,tm+1/2), we consider v(x,tm+1/2) as a convex combination of v(x, t_m) and v(x,t_m+1) such that

Later, a parameter γ in (2.6) will be taken to adjust the convergence of the proposed algorithm. With this, (2.5) can be written as

For the approximation of the values of the solution u(x,t) at (x_j, t_m), we use the notation U^m_j, but for the conservation laws we approximate the cell average of u(x,t_m), rather than u(x_j, t_m), i.e.,

Also we use the notation V^m_j to approximate the cell average of v(x, t_m), i.e.,

Taking the cell-average of the both sides of (2.7), we have

For a conservative scheme for v = a u_x, we employ the Simpson's numerical quadrature

to approximate the integral of v = u_x on time interval [tm,tm+1]

As a result, it follows that

in which we approximate v(x,tm+τ2) with the linear combination of v(x, t_m) and v(x,t_m+1) given by

Here, the positive constants δ, which differ from γ in (2.6), will be chosen later. Using (2.13), we can rewrite (2.12) as

Taking the cell-average for (2.14) on Ij:=[xj−1/2,xj+1/2] (see [14, 15], for example)

leads to the conservative numerical method for v = u_x. Let us set

Then, (2.15) becomes

Hence, using (2.10) for u_t = - v and (2.17), the new conservative method with

can be written as

where the parameters δ and γ will be chosen with proper accuracy in Section 3. Actually, the two constants

will be chosen so that (2.20) and the correction term V^m_j from (2.19) lead to the following algorithm.

It will be shown in Section 3 that the local truncation error of (2.22) is of second-order and that of (2.23) is of first-order. Due to Lemma 2.1, it may be suggested to allow the correction term V^m_j in (2.22) updated by (2.23) as the exact solution to v_t + a v_x = 0. Since the exact solution for the equation vt+avx=0 is v(x,t)=v0(x−at), we get the revised algorithm as following:

Algorithm 2.1. For the linear problem ut+aux=0, with Uj0=1h∫xj−1/2xj+1/2 u 0(x)dx and Vj0=ah[u0(xj+1/2)−u0(xj−1/2)] as initial data, iterates for m=0,1,⋯, and j=0,1,⋯

where the integer index k is chosen to be

Let us denote (2.19) by Vm+1=N(Um,Vm) and (2.20) by Um+1=N(Um,Vm,Vm+1) for the accuracy of convergence where V^m and U^m are replaced with the exact solution u(x,t) and v(x,t) at (x_j, t_m) respectively in the method (2.19) and (2.20). The local truncation errors E(vm):=1τ[N(um,vm)−vm+1] and E(um):=1τ[N(um,vm,vm+1)−um+1] for (2.19) and (2.20), respectively, can be shown using the Taylor expansion of u and v at the point (x_j, t_m).

Assuming that its exact solution is smooth enough, the following relations are hold:

In fact, since v = au_x and vt=−a2uxx.vtt=a3uxxx, we have

and, using ut+v=0,a2uxx=utt=−vt and vtt=a3uxxx=−uttt,

we have

Hence, if one takes the constants δ and γ in (3.1) and (3.2) as

then it follows that for a fixed aτh=:ν

With parameters from (3.3), the algorithm (2.19) and (2.20) becomes

whose combination leads to (2.22) and (2.23).

Summarizing the above arguments, we have

Theorem 3.1. Suppose that the solution is smooth enough. Then the local truncation errors for {Um} and {Vm} generated by the algorithm (3.5) and (3.6) for the problem u_t + au_x =0 are

Hence, (3.5) and (3.6) are first and second-order accuracy respectively.

Let us denote (2.24) as Um+1=N(Um,V0). The local truncation error E(um):=1τ[N(um,v0)−um+1] will be analyzed in the following theorem. For the initial V^0_k, note that the exact solution to vt+avx=0 is v(x,t)=v0(x−at)=au′0(x−at) and that, due to the relation of the index k and j in algorithm (2.24), we have

Theorem 3.2. Suppose that the solution is smooth enough. Then the local truncation error for {Um} generated by Algorithm (2.24) for the problem ut+aux=0 is

Proof. First note that the correction term Vk0 is Vmj in algorithm (2.24) due to (3.7). Therefore, using ut+v=0, aux−v=0, utt=a2uxx and uttt−a3uxxx=0, it follows that

This completes the proof.

We note that using (3.8) one may have a modified equation for ut+aux=0 (see [1, 15, 20] for example). Since v= au_x and utt=a2uxx, the modified equation turns out to be

which is compared to the modified equation ut+aux=−16ah2(1−ν2)uxxx for Lax-Wendroff scheme (see [2, 3, 15] for example).

Theorem 3.3. Let 0≤ν≤1. The sequence {Um} generated by Algorithm for the problem ut+aux=0 satisfies

where ‖ ⋅ ‖1 denotes the l₁ norm of a vector W:=(w1,w2,⋯,wk).

Proof. From the algorithm

it follows that

Hence, one has the conclusion by taking summation..

Note that, according to the above theorems, the new method has almost l₁ contracting property on ‖Um‖1 with second-order accuracy.

We will take a typical linear hyperbolic equation ut+aux=0 (a=1) with several initial data. The numerical solution from the proposed Algorithm will be compared by upwind method, Lax-Wendroff method and several flux methods.

Example 4.1. The first initial condition is given by

By comparing the numerical solutions of the proposed upwind correction scheme (UC) with the classical upwind method and the Lax-Wendroff method, we see the effects of the correction term Vk0 with the weight ν(h−τ). While the numerical behavior of the UC scheme without the correction term are same as those of upwind method with ν², the numerical behavior do reveal little spurious oscillations or smearing of the solution (see Figure 4.1). Next, we compare the numerical solutions of the UC scheme with the Godunov's method for several flux limiters such as minmod, superbee, van Leer and mc. As shown in Figure 4.2, one may notice that the proposed scheme yields much better approximations to exact solution than the Godunov's method yields.

Figure 1. The numerical solutions for the linear equation for the initial data (4.1) at t=1, 2, 3, and 5 with the step size h=1/100.

Figure 2. The numerical solutions for the linear equation for the initial data (4.1) at t = 1, 2, 3, and 5 with the step size h = 1=100.

Example 4.2. The second initial condition will be chosen as the smooth

Comparing the numerical solutions of the UC scheme with the classical upwind method and the Lax-Wendroff method, one may see that the smooth exact solution can be almost exactly approximated comparing to second-order Lax-Wendroff method even if both two methods keep the same second-order accuracy.

Example 4.3. The third initial condition is given by

In this example, the initial data has only one singularity. The numerical solutions of the UC scheme are compared with the Godunov's method for several flux limiters.

The UC scheme shows a better approximation comparing to those of the Godunov's method( see [2, 15] , etc.). Throughout numerous demonstrated examples, one can verify that such better approximations can be obtained. The reason is that the developed method (UC) is of second-order accurate with almost l₁ contracting property.

The newly developed method which works for linear scalar hyperbolic equations uses the upwind algorithm plus a correction term whose weight is chosen as ν(ha−τ). This particular weight is derived by combining the first-order correction algorithm for v = u_x and the second-order algorithm for u_t = -v. It is shown that the proposed upwind correction method has second-order accurate local truncation error with almost l₁ contracting property while a high-order method using flux-limiter function preserves the TVD property.

As done for a linear hyperbolic equation, one may apply the developed algorithm to a nonlinear scalar hyperbolic equation ut+fx(u)=0 with given initial data u(x,0)=u0(x) following known techniques (see [14, 15] for example ). However, in an attempt to expand the developed techniques in this paper to a nonlinear scalar equation directly, one will likely have to find a nice weight to make the computations feasible. Such issues will be discussed in a upcoming paper.

Figure 3. The numerical solutions for the linear equation for the initial data (4.2) at t = 1, 3, 6, and 10 with the step size h = 1=20.

Figure 4. The numerical solutions for the linear equation for the initial data (4.3) at tt = 1, 3, 5, and 9 with the step size h = 1=50.