Semi-Lagrangian scheme for solving hyperbolic equation of first order
Синь В., Вяткин А. В., Шайдуров В. В. Semi-Lagrangian scheme for solving hyperbolic equation of first order // Молодой ученый. 2013. №9. С. 6-13.
We consider numerical solving the hyperbolic equation
equipped with suitable initial and boundary conditions for known velocity coefficient
Among the successful numerical methods for solving this equation we mention such nonoscillatory conservative finite difference shemes as TVD (total variation diminishing), TVB (total variation bounded), ENO (essentially nonoscillatory), and CABARET (Compact Accurate Boundary Adjusting high REsolution Technique) shemes (see -  and the reference there). This approach uses the traditional difference approximation of temporal derivative and is liable to Courant — Friedrichs — Lewy (CFL) condition for ratio between steps in time and space.
The other approach focuses on the approximation of the whole operator of problem (1) as the partial derivative in some direction of the space This approach is developed under different names: methods of trajectories or modified characteristics and semi-Lagrangian one (see, for example, -  and the reference there).
In order to highlight the essential ingredients of suggested approach we shall operate with one-dimensional problem again, keeping in mind that we shall extend this method in subsequent papers. In this paper, we will first show how to consider the possible boundary conditions in contrast to the periodic conditions of the previous paper . Then we modify the presented algorithm for big and huge velocity and give a numerical example that demonstrates the stability independent of CFL condition.
2. The statement of problem and the main theorem
So, in the rectangle consider equation
Assume that coefficient is given at and is positive for simplicity sake.
Unknown function satisfies the initial condition
and boundary condition
Assuming the continuity of at we get at the coordinate origin. Moreover, for continuity of first partial derivatives in (1) the following equality must hold:
By simple calculations we can obtain the necessary condition for the continuity of the second partial derivatives at the point (0, 0), etc. We will not dwell on the question of the sufficiency of this conditions for the smoothness of the solution and at once we assume sufficient smoothness of the velocity and solution for further considerations.
Let us take two time lines with and two nodes
For both these nodes we construct the characteristics of equation (1) at segment [20, 21]. They satisfy the ordinary differential equation with different initial values:
Fig. 1. Trajectories in standard (first) situation
These characteristics define two trajectories in plane when The typical situation is considered in the previous paper  when each of these trajectories crosses line in some point We supposed that they are not mutually crossed and therefore For this case the following result was proved in .
Theorem 1. For smooth solution of equation (1) in the standard (first) situation (Fig. 1) we have equality
But the boundary condition (3) may produce two other situations. First, for small the trajectory may be interrupted at line at point because function is unknown for But trajectory continues up to line (see Fig. 2). And second, both trajectories and are interrupted at line at points and (see Fig. 3). We enumerate these situations from one to three.
Fig. 2. Second situation: trajectory C1 is interrupted and C2 does not
Fig. 3. Third situation: both trajectories C1 and C2 are interrupted
For last situations we prove two equalities.
Theorem 2. For smooth solution of equation (1) with boundary condition (3) in the second situation (Fig. 2) we get the equality
Proof. Define by the curvilinear pentagon bounded by lines And define by the corresponding parts of these lines, which form the boundary (see Fig. 4). Introduce also the external normal defined at each part of boundary except 5 vertices of pentagon.
Fig. 4. Integration along the boundary in second case
Now use formula by Gauss-Ostrogradskii [20, 21] in the following form:
where sing means scalar product. Since the boundary consists of five parts we calculate the integral over separately on each line:
Along the line the external normal equals Then
At arbitrary point the tangent vector is Therefore the external normal (that is orthogonal to it) equals
It implies equality
Along the line the external normal equals Then
For other two parts of the boundary we use the same way to calculate the integrals:
Combination of (7) — (9) and (11) — (14) implies (6):
The next result is justificated like previous theorem with some simplification. Therefore we give it without any proof.
Theorem 3. For smooth solution of equation (1) with boundary condition (3) in the third situation (Fig. 3) we get the equality
Note that this situation contains the case when and coincides with
So, consideration of the boundary conditions in the calculation expressions resulted in an additional calculation of the integrals along of known function In this sense, consideration of boundary conditions makes minor modifications into the algorithms discussed below, so in further considerations we return to the periodic case.
3. Piece-wise linear discrete approximation
So, let the condition of periodicity holds for the problem (1) — (2), namely functions are supposed periodical in with period and are smooth enough for further considerations.
Now we formulate some modification of numerical algorithm from  for solving problem (1) — (2). First, take integer and construct uniform mesh in with nodes and meshsize Than take integer and construct uniform mesh in with nodes and meshsize Then make the following cycle for supposing that the approximate solution is known yet at previous time-level for
1. With the help of values in these points and periodicity we construct the piecewise linear (periodical) interpolant
2. For each point construct approximate trajectories down to time-level for example, by Runge-Kutta method. They produce cross-points If goes outside segment we use periodicity of our data. One can see that we get values which do not coincide with exact values Let we solve equations (4) with the following accuracy:
where is small enough.
3. For each interval compute integral
by trapezoid quadrature formula separately at each nonempty subinterval where is linear.
4. Due to Theorem 1 it is supposed that
Finally we put
Thus, we complete our cycle which may be repeated up to last time-level
Condensed form of this algorithm in terms of piecewise linear periodical interpolants is written as follows:
So, we get approximate discrete solution at each time-level First we prove the conservation law in discrete form.
Let a discrete function is given, and we construct piecewise linear interpolant with period 1.
Theorem 4. For any initial condition the approximate solution (16) — (20) satisfies the equality:
Proof. The justification directly coincides with the proof of Theorem 2 in  with substitution instead of
It is interesting that this discrete conservation law is exactly valid for an approximate values
Now we prove a stability of algorithm (17) — (20) in the discrete norm analogous to that of space
Theorem 5. For any intermediate discrete function the solution of (17) — (20) satisfies the inequality:
Proof. Again the justification directly coincides with the proof of Theorem 3 in  with substitution instead of
Now evaluate an error of approximate solution in introduced discrete norm.
Theorem 4. For sufficiently smooth solution of problem (1) — (2) we have the following estimate for the constructed approximate solution:
with a constant independent of
Proof. We prove this inequality by induction in For this inequality is valid because of exact initial condition (2): Suppose that estimate (25) is valid for some and prove it for
So, at time-level we have decomposition
with a discrete function that satisfies the estimate
Because of Taylor series in of in the vicinity of point we get equality
Because of Theorem 1
Instead of let use its piecewise linear periodical interpolant Then
Thus, we get equality
For we use (21) and (26):
where values of are constructed by piecewise linear periodical interpolation.
Now let us evaluate the difference caused by error (17):
From the properties of piecewise linear interpolant it follows that
Now let subtract (31) from (30), multiply its modulus by use (32), and sum for all
Due to Theorem 3 last term in brackets is combined into Thus
Let put then this inequality is transformed with the help (27):
that is equivalent to (27).
We can see that at last time level we get inequality
In some sense we got a restriction on temporal meshsize to get convergence. For example, to get first order of convergence, it is enough to take
with any constant independent of But this restriction is not such strong for constant as CFL condition:
Moreover, it is opposite in meaning: here the greater the better accuracy.
Thus, this approach is convenient for the problems with huge velocity which come from a computational aerodynamics: we have computational stability on the base of Theorem 4 and conservation law on the base of Theorem 3.
Now discuss the choice of in (17). From (35) the better choice is i.e., For this purpose the standard Runge-Kutta method is acceptable that gives with appropriate accuracy and stability conditions [22, 23].
4. Numerical experiment
Let take and solve the equation (1) with initial condition
subject to periodicity. Then exact solution is
The result of implementing the presented algorithm is given in Table 1. Here in this experiment we set (as the most implemented ratio in computations). The first line shows a number n of mesh nodes and the middle line contains the values at last time level The last line demonstrates the order of accuracy.
Order of accuracy
Thus we indeed have at least the first order of accuracy on h when / h is fixed.
Thus, we continue presentation of the numerical approach  which is more convenient for huge velocity Here we show the treatment with boundary condition instead of periodical one, and then we examine theoretically and by numerical example the effect of the approximate solving the characteristics equations instead of exact process.
Again we have to note that the accuracy is the higher the less time steps done in the algorithm. But for the equations with nonzero right-hand side a small time step will be crucial for appropriate approximation.
1. Harten A., Osher S.: Uniformly high-order accurate non-oscillatory schemes // I. SIAM J. Numer. Anal. — 1987. — V. 24. — P. 279–309.
2. Harten A., Engquist B., Osher S., Chakravarthy S.: Uniformly high-order accurate non-oscillatory schemes, III // J. Comput. Phys. — 1987. — V. 71. — P. 231–303.
3. Osher S., Tadmor E.: On the convergence of different approximations to scalar conservation laws // Math. Comput. — 1988. — V. 50. — P. 19–51.
4. Sanders R.: A third-order accurate variation nonexpansive difference schem for single nonlinear conservation laws // Math. Comput. — 1988. — V. 51. — P. 535–558.
5. Чирков Д. В., Черный С. Г.: Сравнение точности и сходимости некоторых TVD-схем // Вычислительные технологии. — 2000. — Т. 5, № 5. — С. 86–107.
6. Cockburn B., Shu C.-W.: TVB Runge-Kutta projection discontinuous Galerkin finite element method for conservation laws II: general framework // Math. Comput. — 1988. — V. 52. — P. 411–435.
7. Shu C.-W.: TVB boundary treatment for numerical solution of conservation laws // Math. Comput. — 1987. — V. 49. — P. 123–134.
8. Shu C.-W.: Total-Variation-Diminishing time discretizations // SIAM J. Sci. Statist. Comput. — 1988. — V. 9. — P. 1073–1084.
9. Shu C.-W., Osher S.: Efficient implementation of essentially non-oscillatory shock-capturing schemes // J. Comput. Phys. — 1988. — V. 77. — P. 439–471.
10. Sweby P.: High-resolution schemes using flux limiters for hyperbolic conservation laws // SIAM J. Numer. Anal. — 1984. — V. 21. — P. 995–1011.
11. Головизнин В. М., Зайцев М. А., Карабасов С. А., Короткин И. А. Новые алгоритмы вычислительной гидродинамики для многопроцессорных вычислительных комплексов. — Москва: МГУ. — 2013.
12. Priestley A. A quasi-conservative version of the semi-Lagrangian advection scheme // Mon. Weather Rev. — 1993. — V. 121. — P. 621–629.
13. Scroggs J. S., Semazzi F. H. M. Aconservative semi-Lagrangian method for multidimensional fluid dynamics applications // Numer. Meth. Part. Diff. Eq. — 1995. — V. 11. — P. 445–452.
14. Phillips T. N., Williams A. J. Conservative semi-Lagrangian finite volume schemes // Numer. Meth. Part. Diff. Eq. — 2001. — V. 17. — P. 403–425.
15. Iske A. Conservative semi-Lagrangian advection on adaptive unstructured meshes // Numer. Meth. Part. Diff. Eq. — 2004. — V. 20. — P. 388–411.
16. Синь Вэнь, Вяткин А. В., Шайдуров В. В. Characteristics-like approach for solving hyperbolic equation of first order // Молодой ученый. — 2013. — № 3(50). — С. 5–12.
17. Chen H., Lin Q., Shaidurov V. V., Zhou J. Error estimates for triangular and tetrahedral finite elements in combination with trajectory approximation of first derivatives for advection-diffusion equations // Numerical Analysis and Applications. — 2011. — Vol. 4, № 4. — P. 345–362.
18. Чен Х., Лин К., Шайдуров В. В., Жоу Ю. Оценки ошибки для треугольных и тетраэдральных конечных элементов в комбинации с траекторной аппроксимацией первых производных для уравнений адвекции-диффузии // Сиб. журн. вычисл. матем. — 2011. — Т. 14, № 4. — C. 425–442.
19. Shaidurov V. V., Shchepanovskaya G. I., Yakubovich V. M. Numerical simulation of supersonic flows in a channel // Russ. Journ. of Numer. Anal. and Math. Modelling. — 2012. — Vol. 27, № 6. — P. 585–601.
20. Streeter V. L., Wylie E. B.: Fluid mechanics. — London: McGraw-Hill. — 1998.
21. Polyanin A. D.: Handbook of linear partial differential equations for engineers and scientists. — Boca Raton: Chapman & Hall/CRC Press. — 2002.
22. Новиков Е. А.: Явные методы для жестких систем. — Новосибирск: Наука. — 1997.
23. Hairer E., Wanner G., Nørsett S. P.: Solving Ordinary Differential Equations 1: Nonstiff Problems. — Berlin: Springer. — 1993.