Difference scheme, variational

A difference scheme based on a variational problem corresponding to a boundary value problem for a differential equation. The basic idea behind the construction of a variational difference scheme is to obtain a system of linear algebraic equations with the same structure as a system of difference equations (cf. Difference equation) by choosing special coordinate functions in the Ritz method; usually the unknown parameters are approximate values to the exact solution, and, possibly, to certain of its derivatives at the nodes of a grid. It is possible to use piecewise-linear, semi-linear and other functions as such coordinate functions.

Difference schemes can also be obtained by a special choice of the coordinate functions in the Galerkin method. The method in which difference schemes are obtained via Galerkin's method is called the variational difference method (or projection difference method). The variational difference method is sometimes called the finite-element method, although the latter phrase is used also in a more general sense.

Let the boundary value problem

$$ \tag{1 } - ( p ( x) u ^ \prime ( x)) ^ \prime = f ( x),\ \ 0 < x < 1, $$

$$ \tag{2 } u ( 0) = u ( 1) = 0, $$

be given, where $ f $ is a continuous function, $ p $ is continuously differentiable and $ p ( x) \geq p _ {0} > 0 $.

Multiplying (1) by an arbitrary function $ \phi $ satisfying conditions (2) and integrating with respect to $ x $,

$$ - \int\limits _ { 0 } ^ { 1 } ( pu ^ \prime ) ^ \prime \phi dx = \ \int\limits _ { 0 } ^ { 1 } f \phi dx, $$

leads to the identity

$$ \tag{3 } L ( u , \phi ) \equiv \ \int\limits _ { 0 } ^ { 1 } pu ^ \prime \phi ^ \prime dx = \ \int\limits _ { 0 } ^ { 1 } f \phi dx, $$

which is satisfied by the solution to the problem (1), (2). The converse is also true: A function $ u $ satisfying the boundary conditions (2) and the identity (3) for arbitrary functions $ \phi $, $ \phi ( 0) = \phi ( 1) = 0 $, is a solution to the problem (1), (2). The identity (3) is used to construct an approximate solution by the Galerkin method. The interval $ [ 0, 1] $ is divided into $ N $ parts by the points $ x _ {i} = ih $, $ i = 1 \dots N - 1 $, $ h = N ^ {-} 1 $. The set $ \{ x _ {i} \} $ is called the grid, the points $ x _ {i} $ are called the nodes of the grid and $ h $ is called the step of the grid. The coordinate functions in the Galerkin method are taken to be the functions

$$ \phi _ {i} ( x) = \ \psi ( h ^ {-} 1 x - i),\ \ i = 1 \dots N - 1, $$

where

$$ \psi ( t) = \ \left \{ \begin{array}{ll} 1 - | t | & \textrm{ if } | t | \leq 1, \\ 0 & \textrm{ if } | t | > 1. \\ \end{array} \right .$$

Clearly $ \phi _ {i} \equiv 0 $ outside the interval $ [ x _ {i - 1 } , x _ {i + 1 } ] $. This property of coordinate functions is customarily called the property of locality or the property of local support. Let an approximate solution to the problem be sought for in the form

$$ \tag{4 } v = \ \sum _ {i = 1 } ^ { {N } - 1 } v _ {i} \phi _ {i} ( x) , $$

where the required parameters $ v _ {i} $ are the values of the approximate solution at the nodes of the grid:

$$ v _ {i} = \ v ( x _ {i} ),\ \ i = 1 \dots N - 1. $$

Let $ K $ be the set of functions of the form (4). The functions in $ K $ are linear on the intervals $ [ x _ {i} , x _ {i + 1 } ] $, continuous on $ [ 0, 1] $ and vanish for $ x = 0 $ and $ x = 1 $. In the Galerkin method a system is obtained by substituting in (3) the function $ v $ for $ u $ and the functions $ \phi _ {i} $ for $ \phi $:

$$ \tag{5 } L ( v, \phi _ {i} ) = \ \sum _ {j = 1 } ^ { {N } - 1 } v _ {j} \int\limits _ { 0 } ^ { 1 } p ( x) \phi _ {j} ^ \prime ( x) \phi _ {i} ^ \prime ( x) dx = \ \int\limits _ { 0 } ^ { 1 } f ( x) \phi _ {i} ( x) dx, $$

$$ i = 1 \dots N - 1. $$

Here

$$ \phi _ {i} ^ \prime ( x) = \ \left \{ \begin{array}{rl} 1/h & \textrm{ if } x _ {i - 1 } < x < x _ {i} , \\ - 1/h & \textrm{ if } x _ {i} < x < x _ {i + 1 } , \\ \end{array} \right .$$

and

$$ L ( \phi _ {j} , \phi _ {i} ) = \ \int\limits _ { 0 } ^ { 1 } p \phi _ {j} ^ \prime \phi _ {i} ^ \prime dx = 0 $$

only for $ j = 1, i, i + 1 $, so that in each equation there are at most three unknowns.

The system (5) can be written in the form

$$ { \frac{1}{h ^ {2} } } [- \alpha _ {i - 1/2 } v _ {i - 1 } + ( \alpha _ {i - 1/2 } + \alpha _ {i + 1/2 } ) v _ {i} - \alpha _ {i + 1/2 } v _ {i + 1 } ] = f _ {i} , $$

$$ i = 1 \dots N - 1, $$

where

$$ \alpha _ {i - {1 / 2 } } = \ { \frac{1}{h} } \int\limits _ { x _ {i - 1 } } ^ { {x _ i } } p ( x) dx, $$

$$ f _ {i} = { \frac{1}{h} } \int\limits _ {x _ {i - 1 } } ^ { {x } _ {i + 1 } } f ( x) \phi _ {i} ( x) dx,\ v _ {0} = v _ {N} = 0. $$

This system is similar in structure to an ordinary difference system. Systems of equations obtained in this way are also called variational difference schemes. In contrast to the usual difference schemes, the coefficients $ \alpha _ {i - 1/2 } $ and $ f _ {i} $ are not values of the functions $ p $ and $ f $ at fixed points but averages of these. This property enables one to use variational difference schemes for equations with "bad" (for example, discontinuous) coefficients.

Let $ L _ {h} = \{ L ( \phi _ {j} , \phi _ {i} ) \} $ be the matrix of the system (5). Since $ L ( \phi _ {j} , \phi _ {i} ) = L ( \phi _ {i} , \phi _ {j} ) $, the matrix $ L _ {h} $ is symmetric. The equation

$$ ( L _ {h} w _ {h} , w _ {h} ) = \ \int\limits _ { 0 } ^ { 1 } p ( w ^ \prime ) ^ {2} dx $$

holds, where $ w _ {h} = \{ w _ {i} \} _ {i = 1 } ^ {N - 1 } $ is an arbitrary vector in the Euclidean space $ E ^ {N - 1 } $ of dimension $ N - 1 $, $ ( \cdot , \cdot ) $ is the scalar product in $ E ^ {N - 1 } $ and

$$ w = \ \sum _ {i = 1 } ^ { {N } - 1 } w _ {i} \phi _ {i} ( x). $$

From the inequality

$$ \max _ { x } u ^ {2} ( x) = \ \max \left ( \int\limits _ { 0 } ^ { x } u ^ \prime ( t) dt \right ) ^ {2} \leq \ \int\limits _ { 0 } ^ { 1 } ( u ^ \prime ) ^ {2} dx, $$

which is valid for an arbitrary differentiable function $ u $ satisfying $ u ( 0) = 0 $, one deduces the estimate

$$ \int\limits _ { 0 } ^ { 1 } p ( x) ( w ^ \prime ( x)) ^ {2} dx \geq \ p _ {0} \max _ { x } x ( x) ^ {2} \geq \ p _ {0} h \sum _ {i = 1 } ^ { {N } - 1 } w _ {i} ^ {2} , $$

from which one derives the inequality

$$ \tag{6 } ( L _ {h} w _ {h} , w _ {h} ) \geq \ p _ {0} h ( w _ {h} , w _ {h} ). $$

The matrix $ L _ {h} $ is positive definite; the system (5) has a unique solution.

For small values of $ h $ the system (5) consists of a large number of equations. The accuracy of the solution of a system of algebraic equations and the amount of work necessary to find it depends to a large extent on the size of the so-called condition number $ P = \lambda _ \max / \lambda _ \min $ of the matrix of the system. Here $ \lambda _ \max $ and $ \lambda _ \min $ are the largest and the smallest eigen value of $ L _ {h} $. The inequality (6) implies that $ \lambda _ \min \geq p _ {0} h $. The estimate

$$ \lambda _ \max \leq \ { \frac{1}{4h } } \max _ { x } p ( x) $$

is also valid.

The condition number $ P = O ( h ^ {-} 2 ) $, which is the same order in $ h $ as the known estimates for the matrices of ordinary difference schemes.

The convergence of the approximate solution to the exact solution is proved by the usual scheme for the Galerkin method. For an arbitrary function $ \phi $ in $ K $ equations (3) and (5) imply that

$$ \tag{7 } L ( u - v, \phi ) = 0, $$

whence

$$ \tag{8 } L ( u - v, u - v) = \ L ( u - v, u - w), $$

where $ w $ is an arbitrary function in $ K $. The right-hand side of (8) is estimated with the help of the inequality

$$ [ L ( \phi , \psi )] ^ {2} \leq \ L ( \phi , \phi ) L ( \psi , \psi ). $$

Thus,

$$ L ( u - v, u - v) = \ \int\limits _ { 0 } ^ { 1 } p ( x) ( u ^ \prime - v ^ \prime ) ^ {2} dx \leq $$

$$ \leq \ \inf _ {w \in K } \ \int\limits _ { 0 } ^ { 1 } p ( x) ( u ^ \prime - w ^ \prime ) ^ {2} dx. $$

The notation

$$ | u | = \ \left ( \int\limits _ { 0 } ^ { 1 } pu ^ {\prime 2 } dx \right ) ^ {1/2} $$

(the number $ | u | $ is called the energy norm of the function $ u $) enables one to rewrite the latter inequality in the form

$$ | u - v | \leq \ \inf _ {w \in K } \ | u - w | . $$

Estimating the error of the variational difference scheme leads to estimating the best approximation to the exact solution by functions in the class $ K $. If one takes $ w $ to be the piecewise-linear function

$$ w ( x) = \ \sum _ {i = 1 } ^ { {N } - 1 } u ( x _ {i} ) \phi _ {i} ( x), $$

which coincides with the function $ u $ at the nodes of the grid, then the estimate

$$ | u - v | \leq \ Ch ^ {2} \max p ( x) \int\limits _ { 0 } ^ { 1 } ( u ^ {\prime\prime} ) ^ {2} dx $$

is valid for some constant $ C $.

In the example above certain characteristic features of the variational difference method are obvious: the coordinate functions are local, ensuring that the structure of the variational difference scheme is close to that of difference schemes and that one can apply the technique of projection methods to investigate the convergence of the variational difference scheme.

The choice of local coordinate functions having the required approximation properties is basic for the construction of the variational difference scheme. The approximation problem can be posed in various function spaces. For problems in mathematical physics the Sobolev spaces $ W _ {p} ^ {l} ( \Omega ) $ are important, that is, linear sets of functions with finite norm

$$ \| u \| _ {p,l} = \ \sum _ {| \alpha | \leq l } \left ( \int\limits _ \Omega | D ^ \alpha u | ^ {p} \ dx _ {1} \dots dx _ {n} \right ) ^ {1/p} , $$

where $ \Omega $ is a region in $ E ^ {n} $, $ p \geq 1 $, $ l $ is a non-negative integer, $ \alpha = ( \alpha _ {1} \dots \alpha _ {n} ) $ is a vector with integer coordinates, $ | \alpha | = \alpha _ {1} + {} \dots + \alpha _ {n} $, and

$$ D ^ \alpha u = \ \frac{\partial ^ {| \alpha | } u }{\partial x _ {1} ^ {\alpha _ {1} } \dots \partial x _ {n} ^ {\alpha _ {n} } } . $$

Many classes of local coordinate functions can be constructed according to the following scheme. Let functions $ \phi ^ {1} ( \xi ) \dots \phi ^ {r} ( \xi ) $ belonging to $ W _ {p} ^ {k} ( E _ {n} ) $ and vanishing outside the $ n $- dimensional cube $ | \xi _ {j} | < R $, $ j = 1 \dots n $, be given. Let $ h = ( h _ {1} \dots h _ {n} ) $ be a fixed vector with positive coordinates, let $ i = ( i _ {1} \dots i _ {h} ) $ be an arbitrary integer vector and let

$$ h ^ {-} 1 ( x) \equiv \ ( h _ {1} ^ {-} 1 x _ {1} \dots h _ {n} ^ {-} 1 x _ {n} ). $$

Let $ I $ denote the set of vectors $ i $ such that the $ n $- dimensional parallelopipedon $ | h _ {j} ^ {-} 1 x _ {j} - i _ {j} | < R $, $ j = 1 \dots n $, intersects $ \Omega $. For the given region $ \Omega $ as coordinate functions one chooses

$$ \phi _ {i} ^ \mu ( x) = \ \phi ^ \mu ( h ^ {-} 1 x - i),\ \ i \in I,\ \mu = 1 \dots r, $$

that is, functions obtained from the original functions by a change of scale for the variables and a shift by the vector $ i $. Such coordinate functions are called regular. Let the class $ K $ be the set of functions of the form

$$ w ( x) = \ \sum _ {i \in I } \sum _ {\mu = 1 } ^ { r } w _ {i} ^ \mu \phi _ {i} ^ \mu ( x). $$

If any polynomial $ P _ {l - 1 } $ in $ \xi $ of degree $ l $ can be represented as a linear combination of $ \phi ^ {1} ( \xi ) \dots \phi ^ {r} ( \xi ) $, then for an arbitrary function $ u \in W _ {p} ^ {l} ( \Omega ) $ it is possible to find a function $ w \in K $ such that the approximation inequality

$$ \tag{9 } \| u - w \| _ {p,k} \leq \ C \overline{h}\; {} ^ {l - k } \| u \| _ {p,l} ,\ \ 0 \leq k < l, $$

is valid, where $ \overline{h}\; = \max _ {1 \leq j \leq n } h _ {j} $ and $ C $ does not depend on $ h $ and $ u $.

The construction of a variational difference scheme for boundary value problems for elliptic equations is based on equivalent problems of finding functions satisfying integral identities. Many of these problems consist of finding functions $ u \in W _ {2} ^ {m} ( \Omega ) $ satisfying, for an arbitrary function $ \phi \in W _ {2} ^ {m} ( \Omega ) $, an integral identity

$$ \tag{10 } L ( u, \phi ) = \ \sum _ {| \alpha | ,| \beta | \leq m } \ \int\limits _ \Omega a _ {\alpha \beta } D ^ \alpha u D ^ \beta \phi d \Omega + $$

$$ + \sum _ {| \rho | ,| \tau | < m } \int\limits _ { S } b _ {\rho \tau } D ^ \rho u D ^ \tau \phi dS = \int\limits _ \Omega f \phi d \Omega , $$

where $ S $ is the boundary of $ \Omega $ and $ a _ {\alpha \beta } $, $ b _ {\rho \tau } $ and $ f $ are given functions. It is assumed that

$$ a _ {\alpha \beta } = \ a _ {\beta \alpha } ,\ \ b _ {\rho \tau } = \ b _ {\tau \rho } $$

and

$$ L ( u , u ) \geq \gamma \ \| u \| _ {2} ^ {2,m} ,\ \ \gamma = \textrm{ const } > 0. $$

The application of Galerkin's method to (10) using the coordinate functions $ \phi _ {i} ^ \mu $ gives a variational difference scheme for the problem (10). Let the solution $ u $ to (10) belong to $ W _ {2} ^ {l} ( \Omega ) $, $ l > m $, and let the functions $ \phi _ {i} ^ \mu $ satisfy conditions under which inequality (9) is valid. To estimate the error of the variational difference scheme one uses the standard technique for the Galerkin method:

$$ \gamma \| u - v \| _ {2,m} ^ {2} \leq \ L ( u - v, u - v) \leq $$

$$ \leq \ \inf _ {w \in K } L ( u - w, u - w) \leq C h ^ {2 ( l - m) } \| u \| _ {2l} ^ {2} , $$

where $ v $ is the approximate solution.

Problems of the form (10) for which an arbitrary function from $ W _ {2} ^ {m} ( \Omega ) $ can be taken as the function $ \phi $ are called problems with natural boundary conditions. There exists another class of boundary value problems in which on the boundary $ S $ one poses boundary conditions of the form

$$ \tag{11 } \left . \left ( \sum _ {| \alpha | < m } b _ \alpha ^ {j} D ^ \alpha u \right ) \right | _ {S} = 0,\ \ 0 \leq j \leq v \leq m - 1. $$

In this case the boundary conditions (11) should also be satisfied by the functions $ \phi \in W _ {2} ^ {m} ( \Omega ) $ in the identity (10). For an approximate solution of such problems by the Galerkin method it is necessary that the coordinate functions satisfy the conditions (11). The coordinate functions $ \phi _ {i} ^ \mu ( x) $, introduced above by the very method of construction, are, generally speaking, not suitable for representing an approximate solution satisfying conditions (11).

A method for constructing a variational difference scheme for problems with boundary conditions of the form (11) uses the penalty method. Suppose, for example, that it is required to solve the Dirichlet problem for the Poisson equation. This problem is equivalent to determining a function $ u $ with $ u \mid _ {S} = 0 $, which, for an arbitrary function $ \phi $ with $ \phi \mid _ {S} = 0 $, satisfies the integral equation

$$ \int\limits _ \Omega \left ( \frac{\partial u }{\partial x } \frac{\partial \phi }{\partial x } + \frac{\partial u }{\partial y } \frac{\partial \phi }{\partial y } \right ) d \Omega = \ \int\limits _ \Omega f \phi d \Omega . $$

In the penalty method one introduces a function $ v $ which, for an arbitrary function $ \phi $ and $ \epsilon > 0 $, satisfies the integral equation

$$ \int\limits _ \Omega \left ( \frac{\partial v }{\partial x } \frac{\partial \phi }{\partial x } + \frac{\partial v }{\partial y } \frac{\partial \phi }{\partial y } \right ) d \Omega + { \frac{1} \epsilon } \int\limits _ { S } v \phi dS = \ \int\limits _ \Omega f \phi d \Omega , $$

The function $ v $ is the solution to a problem with a natural boundary condition. It is proved that for small values of $ \epsilon $ the solutions $ u $ and $ v $ are close. For an approximate solution to the latter problem it is possible to apply a variational difference scheme using regular coordinate functions.

The general method for constructing coordinate functions is as follows.

For an arbitrary positive number $ h $ let there be given a set of points $ z _ {i} $ in $ \Omega $, $ i = 1 \dots N $, called the nodes of the grid, such that the distance of every point of the domain from some node is at most $ h $. For each node $ z _ {i} $ let a collection of functions $ \phi _ {1} ^ {i} \dots \phi _ {v ( i) } ^ {i} $ in $ W _ {2} ^ {m} ( \Omega ) $ be defined satisfying the given boundary conditions (11), where $ v ( i) \leq M $ and $ M $ does not depend on $ i $ and $ h $. For each $ i $ and all $ j $ let the integral

$$ \int\limits _ \Omega \phi _ {j} ^ {i} \phi _ {m} ^ {k} \ d \Omega $$

be non-zero only for a number of indices $ k $ bounded by some number independent of $ i $ and $ h $( the locality condition for coordinate functions). Let $ K $ be the class of functions of the form

$$ w = \ \sum _ {i = 1 } ^ { N } \sum _ {j = 1 } ^ { {v } ( i) } \lambda _ {j} ^ {i} \phi _ {j} ^ {i} , $$

where $ \lambda _ {j} ^ {i} $ are numerical parameters.

If the solution $ u $ to the boundary value problem can be approximated by functions in $ K $ with accuracy characterized by the inequality

$$ \inf _ {w \in K } \ \| u - w \| _ {2,m} \leq \ C h ^ {l - m } \| u \| _ {2,l} ,\ \ m < l, $$

then for a solution obtained via the variational difference scheme the error estimate

$$ \| u - v \| _ {2,m} \leq \ C h ^ {l - m } \| u \| _ {2,l} $$

is valid.

Irregular grids are sometimes applied to give a more complete calculation of the properties of the problem. For example, for a more accurate reproduction of a function in a neighbourhood of a corner of the boundary it is possible to arrange the nodes on a radial-ring grid.

For numerical computations the matrix of a variational difference scheme should not be too bad conditioned.

For problems of the form (10) the optimal conditioning is considered to be expressed by the relation $ P = O ( N ^ {2m/n} ) $, where $ P $ is the condition number of the matrix of the variational difference scheme, $ N $ is the number of nodes of the grid and $ n $ is the dimension of the space containing the region $ \Omega $. For many practical problems such a conditioning really does hold.

The use of a variational difference scheme combines the merits of the grid method and projection methods. The structure of a variational difference scheme enables one to use economical solution methods of various kinds. The solvability of variational difference schemes is easily established: the matrix of the variational difference scheme is positive definite if the difference operator is positive definite. The question of convergence reduces to the question of the approximation of the exact solution by the coordinate functions of the variational difference scheme and, consequently, the speed of convergence is determined by the differential properties of the exact solution. The variational difference scheme can be applied to problems under very weak restrictions.

Research into variational difference schemes is being carried out in the following basic directions:

1) creating coordinate functions satisfying the boundary conditions, investigating their approximation properties;

2) obtaining estimates for the accuracy in various norms;

3) constructing a variational difference scheme for problems which have some singularities (lines of discontinuity in the coefficients, corners on the boundary, etc.);

4) elaborating methods for solving variational difference schemes and methods for optimizing the methods of solution;

5) solving non-linear equations;

6) applying variational difference schemes to non-stationary equations.

References

[1]	L.A. Oganesyan, V.Ya. Rivkind, L.A. Rukhovets, "Variational-difference methods for solving elliptic equations" , 1–2 , Vilnius (1973–1974) (In Russian)
[2]	G. Fix, "An analyse of the finite element method" , Prentice-Hall (1973)
[3]	J.-P. Aubin, "Approximation of elliptic boundary-value problems" , Wiley (1972) MR0478662 Zbl 0248.65063
[4]	R.G. Varga, "Functional analysis and approximation theory in numerical analysis" , Reg. Conf. Ser. Appl. Math. , 3 , SIAM (1971) (Translated from Russian) MR0310504 Zbl 0226.65064
[5]	S.G. Mikhlin, "Variational-grid approximation" Zap. Nauchn. Sem. Leningrad. Otdel. Mat. Inst. , 48 (1974) pp. 32–188 (In Russian) MR0570097
[6]	G.I. Marchuk, "Methods of numerical mathematics" , Springer (1982) (Translated from Russian) MR0661258 Zbl 0485.65003
[7]	E.G. D'yakonov, "Difference methods for solving boundary value problems" , 1 , Moscow (1971) (In Russian) Zbl 0434.65078 Zbl 0308.65051

References

[a1]	P.G. Ciarlet, "Introduction to the numerical analysis of the finite element method" , North-Holland (1977)
[a2]	P.G. Ciarlet, "The finite element method for elliptic problems" , North-Holland (1978) MR0520174 Zbl 0383.65058
[a3]	A.R. Mitchell, R. Wait, "The finite element method in partial differential equations" , Wiley (1973) MR0615080 MR0483547 MR0488853 MR0383792 Zbl 0344.35001 Zbl 0228.65069