Hyperbolic partial differential equation, numerical methods

Methods for solving hyperbolic partial differential equations using numerical algorithms.

Various mathematical models frequently lead to hyperbolic partial differential equations. Only very infrequently such equations can be exactly solved by analytic methods. The most widely used methods are numerical methods. They find extensive use in solving problems of the mechanics of a continuous medium, in particular for the equations of gas dynamics (cf. Gas dynamics, numerical methods of), which are quasi-linear.

Numerical methods for solving hyperbolic partial differential equations may be subdivided into two groups: 1) methods involving an explicit separation of the singularities of the solution; 2) indirect computation methods, in which the singularities are not directly separated but are obtained in the course of the computation procedure as domains with sharp changes in the solutions.

The first group includes, for instance, the method of characteristics, which is only used for solving hyperbolic partial differential equations; it has found extensive use in solving the problems of gas dynamics.

The methods in the second group yield non-singular difference schemes (cf. Difference scheme). Let there be given, for example, the hyperbolic equation

$$ \tag{1 } \frac{\partial w }{\partial t } = A \frac{\partial w }{\partial x } , $$

where $ A $ is an $ ( m \times m ) $- matrix with $ m $ distinct real eigen values and $ w = w( x, t) $ is a vector function with $ m $ components. The matrix $ A $ may either be a function in $ x $, $ t $( in such a case (1) is a linear hyperbolic system), or may depend on $ w = w ( x, t ) $( a quasi-linear system). In the latter case, let the system of equations (1) be reduced to divergence form:

$$ \frac{\partial w }{\partial t } = \ \frac{\partial F }{\partial x } + \psi , $$

where $ F $ is a vector function in $ w $, $ x $, $ t $ such that $ A = dF/dw $ and $ \psi $ is a vector function in $ w $, $ x $, $ t $. In the most important case $ A $, $ F $ and $ \psi $ depend only on $ w $. For (1) one may pose the Cauchy problem

$$ w ( x, 0) = w _ {0} ( x) $$

with suitable boundary conditions.

As a rule, the basis on which finite-difference schemes are constructed is an approximation by the integral conservation law corresponding to the differential equation, using certain quadrature formulas on the contour of integration of the differential cell. If the solutions are smooth, approximation by the integral conservation law is equivalent to direct approximation of the corresponding differential equation. Finite-difference schemes must satisfy the approximation and stability requirements. These requirements are mutually independent and, in a sense, contradict one another. In divergence systems of differential equations it is the divergence (or the conservation) condition of the finite-difference scheme which is fundamental. Moreover, the finite-difference schemes must also satisfy a number of other necessary conditions, like dissipativeness, economy, etc. A two-level explicit finite-difference scheme for a linear equation of type (1) has the form

$$ w _ {j} ^ {n + 1 } = \Lambda w _ {j} ^ {n} , $$

where $ \Lambda $ is a finite operator, i.e. can be represented in the form

$$ \Lambda = \sum _ { \alpha = - q _ {1} } ^ { {q _ 2 } } B _ \alpha T _ {1} ^ \alpha , $$

where $ B _ \alpha $ are $ ( m \times m) $- matrices with coefficients that depend on $ \tau $, $ h $, $ x $, $ t $; $ t = n \tau $, $ x = jh $; $ \tau $, $ h $ are the steps of the finite-difference grid along the axis $ t $ and $ x $, respectively; the numbers $ q _ {1} $ and $ q _ {2} $ do not depend on $ \tau $, $ h $; $ \kappa = \tau / h $; and $ T _ {1} $ is a displacement operator with respect to $ x $.

The consistency conditions lead to the relationships

$$ \sum _ { \alpha = - q _ {1} } ^ { {q _ 2 } } B _ \alpha = I, $$

$$ \sum _ { \alpha = - q _ {1} } ^ { {q _ 2 } } \alpha B _ \alpha = \kappa A, $$

where $ I $ is the identity matrix.

The implicit finite-difference scheme may be written as follows:

$$ \Lambda _ {1} w _ {j} ^ {n + 1 } = \Lambda _ {0} w _ {j} ^ {n} , $$

where $ \Lambda _ {1} $ and $ \Lambda _ {0} $ are finite operators,

$$ \Lambda _ {k} = \ \sum _ { \alpha = - q _ {1} } ^ { {q _ 2 } } B _ \alpha ^ {k} T _ {1} ^ \alpha ,\ \ k = 0, 1, $$

where $ B _ \alpha ^ {k} $ are $ ( m \times m ) $- matrices which depend on $ \tau $, $ h $, $ x $, $ t $, and the operator $ \Lambda _ {1} $ contains at least two non-zero matrices $ B _ \alpha ^ {1} $. The operator $ \Lambda _ {1} $ is assumed to be invertible, but its inverse need not be finite.

The various schemes may be subdivided into two classes by their approximation properties: the conditionally approximating schemes and the absolutely-approximating schemes. The conditionally-approximating finite-difference schemes approximate the initial differential equation for $ \tau $, $ h $ tending to zero with a certain interdependence between them: $ \tau = \phi ( h) $. Absolutely-approximating finite-difference schemes approximate the initial differential equation as $ \tau $, $ h $ tend to zero in accordance with an arbitrary law.

In the case of conditional approximation the finite-difference equation may approximate various differential equations for various conditions of limit transition. For instance, for the equation

$$ \tag{2 } \frac{\partial u }{\partial t } + a \frac{\partial u }{\partial x } = 0,\ \ a = \textrm{ const } > 0, $$

one may consider, e.g., the two finite-difference schemes

$$ \tag{3 } \left . \begin{array}{c} { \frac{u _ {j} ^ {n + 1 } - \overline{u}\; {} _ {j} ^ {n} } \tau } + a { \frac{u _ {j + 1 } ^ {n} - u _ {j - 1 } ^ {n} }{2h} } = 0, \\ \overline{u}\; {} _ {j} ^ {n} = { \frac{1}{2} } ( u _ {j + 1 } ^ {n} + u _ {j - 1 } ^ {n} ), \\ \end{array} \right \} $$

$$ \tag{4 } { \frac{u _ {j} ^ {n + 1 } - u _ {j} ^ {n} } \tau } + a { \frac{u _ {j} ^ {n} - u _ {j - 1 } ^ {n} }{h} } = 0. $$

If the law governing the limit transition is

$$ { \frac \tau {h} } = \textrm{ const } , $$

then the finite-difference scheme (3) approximates equation (2), while if that law is

$$ { \frac \tau {h} ^ {2} } = \textrm{ const } , $$

it approximates the equation

$$ \frac{\partial u }{\partial t } + a \frac{\partial u }{\partial x } = \mu \frac{\partial ^ {2} u }{\partial x ^ {2} } ,\ \ \mu = \frac{h ^ {2} }{2 \tau } . $$

The finite-difference scheme (4) absolutely approximates equation (2).

In a similar manner, finite-difference schemes may be subdivided into conditionally-stable schemes and absolutely-stable schemes. Thus, the finite-difference scheme (4) will be stable if the following condition (Courant's condition) is met:

$$ \frac{\tau a }{h } \leq 1, $$

i.e. if it is conditionally stable. On the other hand, the implicit finite-difference scheme

$$ \frac{u _ {j} ^ {n + 1 } - u _ {j} ^ {n} } \tau + a \frac{u _ {j + 1 } ^ {n + 1 } - u _ {j - 1 } ^ {n + 1 } }{2h } = 0 $$

is stable for all relations between $ \tau $ and $ h $, i.e. it is absolutely stable.

Explicit finite-difference schemes are simple in realization, but are either conditionally stable or conditionally approximating. In the case of an absolutely-approximating finite-difference scheme the condition of stability of the explicit scheme usually has the form

$$ \tau \leq \textrm{ const } \cdot h ^ \beta \ \ ( \beta \geq 1), $$

which results in too small a step $ \tau $ and involves unjustifiable labour in computations. All absolutely-stable absolutely-approximating schemes are contained in the class of implicit schemes.

Implicit finite-difference schemes are more complex in realization on passing from one time level to another, but then the step size $ \tau $ may often be chosen arbitrary large, and is thus determined by considerations of accuracy only.

The convergence theorems for finite-difference schemes approximating linear differential equations make it possible to reduce the study of the convergence of such schemes to a study of their stability.

The study of the approximation of a finite-difference scheme corresponding to a hyperbolic equation is rather simple in the case of smooth solutions, has a local character, and in fact amounts to expansion into Taylor series; if the solution is discontinuous, the problem becomes difficult and consists of the verification of integral conservation laws. The study of stability is much more complicated.

The stability of finite-difference schemes which approximate hyperbolic equations with constant coefficients is studied by the Fourier method — i.e. an estimate is made of the norm of the Fourier transform of the step operator of the finite-difference scheme. Since the spectral radius of the matrix of the Fourier transform of the step operator does not exceed the norm of the matrix, a necessary criterion for stability follows: For a finite-difference scheme to be stable it is necessary that the spectral radius of the Fourier transform of the step operator does not to exceed $ 1 + O ( \tau ) $, where $ \tau $ is the step of the scheme along the $ t $- axis. This is also a necessary condition for finite-difference schemes with variable coefficients and, if certain supplementary restrictions are imposed, this is also a sufficient condition. The following methods are used in the study of finite-difference schemes with variable coefficients, as well as for certain non-linear equations: the method of majorizing or of a priori estimates, and a local algebraic method.

The method of a priori estimates is analogous to the corresponding method for differential equations, but in the case of finite differences its realization involves major difficulties, owing to the specific features of finite-difference analysis in which — unlike in the method of a priori estimates in the theory of differential equations — many relationships take a tedious form.

The simplest majorizing estimate is the estimate for finite-difference schemes with positive coefficients.

For instance, consider the following finite-difference scheme for equation (2) with $ a = a( x) $:

$$ \tag{5 } { \frac{u _ {j} ^ {n + 1 } - u _ {j} ^ {n} } \tau } + a _ {j} { \frac{u _ {j} ^ {n} - u _ {j - 1 } ^ {n} }{h} } = 0,\ \ a _ {j} = a ( jh). $$

Then, if

$$ 0 \leq 1 - \kappa _ {j} \leq 1,\ \ \kappa _ {j} = { \frac{\tau a _ {j} }{h} } , $$

the following estimate holds:

$$ \| u ^ {n + 1 } \| \leq \| u ^ {n} \| , $$

where

$$ \| u ^ {n} \| = \max _ { j } | u _ {j} ^ {n} | . $$

It follows that the scheme (5) is uniformly stable in the space $ C $. This estimate can be applied to finite-difference schemes which approximate hyperbolic systems of equations in invariants.

A very important, though limited, class of finite-difference schemes is represented by schemes with positive coefficients and matrices (so-called majorant schemes). If the coefficients of such schemes are symmetric, positive and Lipschitz continuous with respect to $ x $, then they are stable in the space $ L _ {2} $. As a rule, these are finite-difference schemes of first-order approximation in which the derivatives are approximated by one-sided differences. In higher-order approximations, when central differences are taken, one need not obtain positive coefficients. In such cases a priori estimates of a more general type in the spaces $ W _ {2} ^ {p} $ are used.

For instance, let (1) be an acoustic scheme where

$$ w = \left ( \begin{array}{c} u \\ v \end{array} \right ) ,\ \ A = \left ( \begin{array}{cc} 0 & a \\ 1 & 0 \\ \end{array} \right ) ; $$

and let the functions $ a( x, t) $, $ u( x, t) $ and $ v( x, t) $ be periodic with respect to $ x $ with period $ 2 \pi $. The a priori estimate for the finite-difference scheme

$$ \tag{6 } \left . \begin{array}{c} { \frac{u ^ {n + 1 } ( x) - u ^ {n} ( x) } \tau } = a { \frac{v ^ {n + 1 } ( x + h) - v ^ {n + 1 } ( x) }{h} } , \\ { \frac{v ^ {n + 1 } ( x) - v ^ {n} ( x) } \tau } = \ { \frac{u ^ {n + 1 } ( x) - u ^ {n + 1 } ( x - h) }{h} } \\ \end{array} \right \} $$

has the form

$$ \| E ^ {n + 1 } \| ^ {2} \leq \ \frac{1 + b _ {1} \tau }{1 + b _ {2} \tau } \| E ^ {n} \| ^ {2} , $$

where

$$ \| E \| ^ {2} = \ \int\limits _ {- \pi } ^ \pi ( u ^ {2} + av ^ {2} ) dx, $$

$$ b _ {1} = \sup _ {x, n, \tau } \ \left | \frac{a ^ {n + 1 } ( x) - a ^ {n} ( x) }{\tau a ^ {n} ( x) } \right | , $$

$$ b _ {2} = \sup _ {x, n, h } \left ( { \frac{1}{\sqrt {a ^ {n + 1 } ( x + h) }} } \left | \frac{a ^ {n + 1 } ( x + h) - a ^ {n + 1 } ( x) }{h} \right | \right ) . $$

The given estimate proves the stability of the finite-difference scheme (6) and is analogous to the energy inequality for the system of acoustic equations.

The local algebraic method is based on the study of the properties of the local finite-difference operator obtained from the respective finite-difference operator with variable coefficients by "freezing" the coefficients. In this way the analysis of the stability of a finite-difference operator with variable coefficients is replaced by the analysis of an entire family of operators with constant coefficients. The local stability criterion is a generalization of the method of "freezing" the coefficients employed in the theory of differential equations.

The local stability criterion is closely connected with the dissipative stability criterion. A finite-difference scheme is called dissipative of order $ \nu $, where $ \nu $ is an even number, if there exists a $ \delta > 0 $ such that

$$ | \rho | \leq 1 - \delta | \xi | ^ \nu , $$

$$ | \xi | = | kh | \leq \pi , $$

where $ \rho $ is the eigen value of maximum modulus of the transition matrix (the Fourier transform of the step operator) of the finite-difference scheme, and $ k $ is a dual variable. Then if the order of approximation of the scheme is $ 2r + 1 $( $ 2r + 2 $), $ r = 0, 1 \dots $ and if it is dissipative of order $ 2r + 2 $( $ 2r + 4 $), the scheme will be stable in $ L _ {2} $ for hyperbolic systems of first-order differential equations with a Hermitian matrix.

In studying stability of finite-difference schemes for non-linear hyperbolic equations (in particular, for the equations of gas dynamics), the differential approximation method, in which the analysis of the finite-difference scheme is replaced by the analysis of its differential approximation, is employed.

For instance, the differential approximation of the finite-difference scheme (4) for equation (2) is constructed as follows. The expansion in (4) of the functions

$$ u _ {j} ^ {n + 1 } = u ( jh, ( n + 1) \tau ), $$

$$ u _ {j - 1 } ^ {n} = u (( j - 1) h, n \tau ) $$

into Taylor series with respect to the point $ x = jh $, $ t = n \tau $ by the parameters $ \tau $ and $ h $ yields the $ \Gamma $- form of the differential representation of the finite-difference scheme:

$$ \tag{7 } \frac{\partial u }{\partial t } + \sum _ {l = 2 } ^ \infty { \frac{\tau ^ {l - 1 } }{l!} } \frac{\partial ^ {l} u }{\partial t ^ {l} } + a \frac{\partial u }{\partial x } + a \sum _ {l = 2 } ^ \infty { \frac{(- 1) ^ {l - 1 } h ^ {l} }{l!} } \frac{\partial ^ {l} u }{\partial x ^ {l} } = 0. $$

The elimination of the derivatives

$$ \frac{\partial ^ {2} u }{\partial t ^ {2} } ,\ \frac{\partial ^ {3} u }{\partial t ^ {3} } \dots $$

from (7) yields the $ \Pi $- form of the differential representation of the scheme (4):

$$ \tag{8 } \frac{\partial u }{\partial t } + a \frac{\partial u }{\partial x } = \ \sum _ {l = 2 } ^ \infty C _ {l} \frac{\partial ^ {l} u }{\partial x ^ {l} } , $$

where $ C _ {l} $ are certain coefficients which depend on $ \tau $, $ h $, $ a $, and $ C _ {l} = O ( \tau ^ {l-} 1 , h ^ {l-} 1 ) $. Elimination of the terms of order $ O ( \tau ^ {2} , h ^ {2} ) $ from (7) and (8) yields, respectively, the $ \Gamma $- form of the first differential approximation of the scheme (4):

$$ { \frac \tau {2} } \frac{\partial ^ {2} u }{\partial t ^ {2} } + \frac{\partial u }{\partial t } + a \frac{\partial u }{\partial x } - { \frac{ah}{2} } \frac{\partial ^ {2} u }{\partial x ^ {2} } = 0 , $$

and the $ \Pi $- form of the first differential approximation of the scheme (4):

$$ \tag{9 } \frac{\partial u }{\partial t } + a \frac{\partial u }{\partial x } = \ C _ {2} \frac{\partial ^ {2} u }{\partial x ^ {2} } , $$

$$ C _ {2} = { \frac{ah}{2} } \left ( 1 - \frac{a \tau }{h} \right ) . $$

It has been proved in the linear case, for several finite-difference schemes, that if the first differential approximation is correct, then the respective finite-difference scheme is stable. Thus, in the case of scheme (4) above, the fact that (9) is correct means that $ C _ {2} \geq 0 $, i.e. that a necessary and sufficient condition, $ a \tau /h \leq 1 $, for the stability of the scheme is satisfied. The terms with even derivatives in equation (8) ensure the dissipative properties of the finite-difference scheme, while those with odd derivatives are responsible for its dispersive properties.

The dissipation of the finite-difference scheme (4) is the magnitude

$$ d = 1 - | \rho | ^ {2} = \ 1 - \mathop{\rm exp} \left [ \sum _ {l = 1 } ^ \infty (- 1) ^ {l} \frac \tau {h ^ {2l} } C _ {2l} \xi ^ {2l} \right ] , $$

where

$$ \rho = 1 - \frac{a \tau }{h} ( 1 - e ^ {- i \xi } ) $$

is the amplification factor of the scheme; the dispersion of the scheme (4) is the magnitude

$$ \kappa = \kappa a \xi - \mathop{\rm arg} \rho = \ \sum _ {l = 1 } ^ \infty (- 1) ^ {l + 1 } \frac \tau {h ^ {2l + 1 } } C _ {2l + 1 } \xi ^ {2l + 1 } . $$

The dissipative terms in (8) determine the properties of approximate viscosity of the scheme (i.e. a certain mechanism of smoothing in the finite-difference scheme). The form of the dissipative terms is determined both by artificial terms introduced into the initial differential equation and by the structure of the finite-difference scheme itself. The first differential approximation yields the principal term of the approximate viscosity. The method of differential approximation is widely employed in the study of differential schemes for non-linear equations and makes it possible to explain the instability effects of various finite-difference schemes which can be encountered in specific computations and which are not locally revealed by the Fourier method.

The construction of finite-difference schemes in multi-dimensional cases is based on splitting methods (weak approximation) and fractional-step methods, by which the integration of the initial multi-dimensional equation is reduced to the integration of equations of a simpler structure (cf. Fractional steps, method of).

The method of solving hyperbolic equations may be developed by using the finite-element method, which may be considered as a finite-difference method on a special irregular grid.

References

[1]	S.K. Godunov, V.S. Ryaben'kii, "The theory of difference schemes" , North-Holland (1964) (Translated from Russian) MR0163421
[2]	R.D. Richtmeyer, K.W. Morton, "Difference methods for initial value problems" , Wiley (1967) MR220455
[3]	B.L. Rozhdestvenskii, N.N. Yanenko, "Systems of quasilinear equations and their applications to gas dynamics" , Amer. Math. Soc. (1983) (Translated from Russian) MR0694243
[4]	A.A. Samarskii, A.V. Gulin, "Stability of difference schemes" , Moscow (1973) (In Russian) MR2263771 MR1800050 MR1473135 MR1441771 MR0405838 MR0431671 MR0314262 MR0268694 MR0264870 MR0232229 MR0146975 MR0156471 Zbl 1002.65101 Zbl 1189.34107 Zbl 0939.65082 Zbl 0814.65092 Zbl 0368.65032 Zbl 0368.65031 Zbl 0334.65072 Zbl 0303.65079 Zbl 0313.65081
[5]	N.N. Yanenko, Yu.I. Shokin, "The first differential approximation to finite-difference schemes for hyperbolic systems of equations" Siberian Math. J. , 10 : 5 (1969) pp. 868–880 Sibirsk. Mat. Zh. , 10 : 5 (1969) pp. 1173–1187 Zbl 0196.11701
[6]	S.I. Serdyukova, "A necessary and sufficient condition for the stability of a class of difference boundary-value problems" Soviet Math. Dokl. , 14 : 1 (1973) pp. 50–54 Dokl. Akad. Nauk SSSR , 208 : 1 (1973) pp. 52–55

Comments

The Fourier method for analyzing stability of difference schemes is often called normal mode analysis. A useful reference in this connection is [a8].

There is a rapidly growing theory devoted to finding physically meaningful and acceptable approximations of, e.g., shocks (cf. Shock waves, mathematical theory of). In particular, the introduction of total-variation diminishing schemes (TVD-schemes) and flux-limiter schemes have resulted in finite-difference methods that are second order in space and time, and have these TVD or limited flux properties as well, when applied to conservation laws.

References

[a1]	G.E. Forsythe, W.R. Wasow, "Finite difference methods for partial differential equations" , Wiley (1960) MR0130124 Zbl 0099.11103
[a2]	P.R. Garabedian, "Partial differential equations" , Wiley (1964) MR0162045 Zbl 0124.30501
[a3]	I. Gladwell (ed.) R. Wait (ed.) , A survey of numerical methods for partial differential equations , Clarendon Press (1979) MR0569444 Zbl 0417.65047
[a4]	A.R. Mitchell, D.F. Griffiths, "The finite difference method in partial differential equations" , Wiley (1980) MR0562915 Zbl 0417.65048
[a5]	A.A. Samarskii, "Theorie der Differenzverfahren" , Akad. Verlagsgesell. Geest u. Portig K.-D. (1984) (Translated from Russian)
[a6]	G.D. Smith, "Numerical solution of partial differential equations" , Oxford Univ. Press (1977) MR0827497 MR0509636 MR0391467 MR1534099 MR0282538 Zbl 0576.65089 Zbl 0389.65040 Zbl 0123.11806
[a7]	N.N. Yanenko, "The method of fractional steps; the solution of problems of mathematical physics in several variables" , Springer (1971) (Translated from Russian) MR0307493 Zbl 0356.65109 Zbl 0209.47103 Zbl 0183.18201 Zbl 0099.33502
[a8]	R. Vichnevetsky, J.B. Bowles, "Fourier analysis of numerical approximations of hyperbolic equations" , SIAM (1982) MR0675265 Zbl 0495.65041
[a9]	A. Harten, "On a class of high resolution total-variation-stable finite-difference schemes" SIAM J. Numer. Anal. , 21 (1984) pp. 1–23 MR0731210