Linear boundary value problem, numerical methods

Methods that make it possible to obtain a solution of a linear boundary value problem in the form of a table of approximate values of it at points of a grid, without using preliminary information about the expected form of the solution. In the theory of these methods it is typical to assume that the solution of the original problem exists and has sufficiently many derivatives. Because of the absence of other assumptions numerical methods differ in their universality.

Fundamental to numerical methods for solving a linear boundary value problem is the replacement of the original system of equations by its grid approximation. In the case of integro-differential equations such an approximation is usually constructed by means of difference schemes and quadrature formulas. The following problems arise:

1) How quickly does the exact solution of the grid problem converge to the solution of the original problem under refinement of the grid?

2) How sensitive is the solution of the grid problem to changes in the initial data?

3) How to find, at least approximately, the solution of the grid problem?

In the solution of the first and second problems the following apparatus is used. Suppose that in a closed domain $ D $ one is given the equation

$$ \tag{1 } L u = f , $$

and on its boundary $ \Gamma $, consisting of components $ \Gamma _ {i} $, one is given the boundary conditions

$$ \tag{2 } l _ {i} u = \phi _ {i} \ \ \mathop{\rm on} \Gamma _ {i} ,\ \ i = 1 \dots s . $$

Here the given continuous functions $ f $ and $ \phi _ {i} $ belong to some normed linear spaces $ F $ and $ \Phi _ {i} $, respectively, and the linear operators $ L $ and $ l _ {i} $ transform a certain linear subset $ W $ of the normed linear space $ U $ of functions that are continuous in $ D $ into $ F $ and $ \Phi _ {i} $, respectively.

Suppose that for an approximation of equation (1) one has chosen grids $ D _ {h} $ and $ D _ {h} ^ {0} \subset D $ that depend on a positive parameter $ h $, normed linear spaces $ U _ {h} $ and $ F _ {h} $ of functions $ u _ {h} $ and $ f _ {h} $ defined at the points of the grids $ D _ {h} $ and $ D _ {h} ^ {0} $, respectively, and also a linear operator $ L _ {h} $ that transforms $ U _ {h} $ into $ F _ {h} $. For an approximation of the boundary condition (2) one chooses a grid $ \Gamma _ {ih} \subset \Gamma $, a normed linear space $ \Phi _ {ih} $ of functions $ \phi _ {ih} $ defined at the points of the grid $ \Gamma _ {ih} $, and a linear operator $ l _ {ih} $ that transforms $ U _ {h} $ into $ \Phi _ {ih} $. Furthermore, one chooses linear operators $ P _ {h} $ and $ p _ {ih} $ that transform $ F $ and $ \Phi _ {i} $ into $ F _ {h} $ and $ \Phi _ {ih} $, respectively. As a result one obtains a grid approximation

$$ \tag{3 } L _ {h} u _ {h} = P _ {h} f , $$

$$ \tag{4 } l _ {ih} u _ {h} = p _ {ih} \phi _ {i} \ \mathop{\rm on} \Gamma _ {ih} ,\ i = 1 \dots s , $$

of the original problem (1), (2). Let $ [ u ] _ {h} $, $ \{ f \} _ {h} $, $ \{ \phi _ {i} \} _ {ih} $ denote the traces (restrictions) of the functions $ u $, $ f $, $ \phi _ {i} $ belonging to $ U $, $ F $, $ \Phi _ {i} $ on the grids $ D _ {h} $, $ D _ {h} ^ {0} $, $ \Gamma _ {ih} $, respectively. It makes sense to carry out an investigation of the convergence of the grid approximations $ u _ {h} $ to the solution $ u $ of the original problem as $ h \rightarrow 0 $ only in grid norms compatible with the norms in $ U $, $ F $, $ \Phi _ {i} $, that is, on the condition that for any functions $ u $, $ f $, $ \phi _ {i} $ the limit relations

$$ \| [ u ] _ {h} \| _ {U _ {h} } \rightarrow \| u \| _ {U} ,\ \ \| \{ f \} _ {h} \| _ {F _ {h} } \rightarrow \| f \| _ {F} , $$

$$ \| \{ \phi _ {i} \} _ {ih} \| _ {\Phi _ {ih} } \rightarrow \| \phi _ {i} \| _ {\Phi _ {i} } $$

hold. If $ u \in W $, $ f \in F $, $ \phi _ {i} \in \Phi _ {i} $, then the quantity

$$ z ( h) = \| L _ {h} [ u ] _ {h} - \{ L u \} _ {h} \| _ {F _ {h} } + \| P _ {h} f - \{ f \} _ {h} \| _ {F _ {h} } + $$

$$ + \sum _ { i= } 1 ^ { s } \| l _ {ih} [ u ] _ {h} - \{ l _ {i} u \} _ {ih} \| _ { \Phi _ {ih} } + \| p _ {ih} \phi _ {i} - \{ \phi _ {i} \} _ {ih} \| _ {\Phi _ {ih} } $$

is called the error of approximation of the problem (1), (2) by the problem (3), (4). If $ u $ is the solution of the problem (1), (2), then the quantity

$$ \rho ( h) = \| L _ {h} [ u] _ {h} - P _ {h} f \| _ {F _ {h} } + \sum _ { i= } 1 ^ { s } \| l _ {ih} [ u] _ {h} - p _ {ih} \phi _ {i} \| _ {\Phi _ {ih} } $$

is called the error of approximation on this solution. One says that the problem (3), (4) approximates the problem (1), (2) (or approximates it on the solution $ u $) if $ z ( h) \rightarrow 0 $( or $ \rho ( h) \rightarrow 0 $) as $ h \rightarrow 0 $. The order of smallness of $ z ( h) $( or $ \rho ( h) $) is called the order of approximation (or the order of approximation on the solution $ u $).

Suppose, for example, that the heat equation is approximated. In the quadrant $ D $( $ x \geq 0 , t \geq 0 $) equation (1) is to be solved, where

$$ L u = \ \frac{\partial u }{\partial t } - \frac{\partial ^ {2} u }{\partial x ^ {2} } , $$

under the initial condition $ l _ {1} u = \phi _ {1} ( x) $ on the half-line $ \Gamma _ {1} $( $ t = 0 , x \geq 0 $) and the boundary condition $ l _ {2} u = \phi _ {2} ( t) $ on the half-line $ \Gamma _ {2} $( $ x = 0 , t \geq 0 $), where $ l _ {1} u = u $, $ l _ {2} u = \partial u / \partial x $, and $ \phi _ {1} ^ \prime ( 0) = \phi _ {2} ( 0) $. For $ F $ one chooses the space of continuous functions $ f ( t , x ) $ with the norm

$$ \| f \| _ {F} = \ \sup _ { D } \ | f ( t , x ) | , $$

and for $ \Phi _ {1} $ and $ \Phi _ {2} $ one chooses the spaces of continuous functions $ \phi _ {1} ( x) $ and $ \phi _ {2} ( t) $ with the norms

$$ \| \phi _ {1} \| _ {\Phi _ {1} } = \ \sup _ {\Gamma _ {1} } | \phi _ {1} ( x) | \ \textrm{ and } \ \ \| \phi _ {2} \| _ {\Phi _ {2} } = \ \sup _ {\Gamma _ {2} } \ | \phi _ {2} ( t) | , $$

respectively. One may assume that $ U $ coincides with $ F $, and for $ W $ one chooses the subset of functions $ u \in U $ that have continuous partial derivatives $ \partial u / \partial t $, $ \partial u / \partial x $, $ \partial ^ {2} u / \partial x ^ {2} $ in $ D $. Let $ D _ {h} $ be the set of points $ ( n \tau , m h ) $, $ n = 0 , 1 , . . . $; $ m = - 1 , 0 , 1 \dots $ where the steps of the grid are connected by a functional dependence of the form $ \tau = \tau ( h) $ such that $ \tau \rightarrow 0 $ as $ h \rightarrow 0 $. The set $ D _ {h} ^ {0} $ will consist of the points $ ( ( n + 1 /2 ) \tau , m h ) $, $ n = 0 , 1 , . . . $; $ m = 0 , 1 , . . . $. For $ F _ {h} $ and $ U _ {h} $ one chooses the spaces of grid functions $ f _ {h} ( t , x ) $ and $ u _ {h} ( t , x ) $ with the compatible norms

$$ \| f _ {h} \| _ {F _ {h} } = \ \sup _ {D _ {h} ^ {0} } \ | f _ {h} ( t , x ) | \ \ \textrm{ and } \ \ \| u _ {h} \| _ {U _ {h} } = \ \sup _ {D _ {h} } \ | u _ {h} ( t , x ) | . $$

The operator $ L _ {h} $ can be defined by the relation

$$ L _ {h} u _ {h} ( t , x ) = \frac{1} \tau \left [ u _ {h} \left ( t + \frac \tau {2} , x \right ) - u _ {h} \left ( t - \frac \tau {2} , x \right ) \right ] + $$

$$ - \frac{1}{2 h ^ {2} } \left [ u _ {h} \left ( t + \frac \tau {2} ,\ x + h \right ) - 2 u _ {h} \left ( t + \frac \tau {2} , x \right ) + u _ {h} \left ( t + \frac \tau {2} , x - h \right ) \right ] + $$

$$ - \frac{1}{2 h ^ {2} } \left [ u _ {h} \left ( t - \frac \tau {2} ,\ x + h \right ) - 2 u _ {h} \left ( t - \frac \tau {2} , x \right ) + u _ {h} \left ( t - \frac \tau {2} , x - h \right ) \right ] ; $$

let

$$ P _ {h} f = \{ f \} _ {h} \ \ \left ( \textrm{ or } P _ {h} f ( t , x ) = \ \int\limits _ { x- } h ^ { x+ } h f ( t , \xi ) d \xi \right ) . $$

For an approximation of the initial condition one chooses the grid $ \Gamma _ {1h} $ consisting of the points $ ( 0 , m h ) $, $ m = 0 , 1 \dots $ the space $ \Phi _ {1h} $ of grid functions $ \phi _ {1h} ( x) $ with the compatible norm

$$ \| \phi _ {1h} \| _ {\Phi _ {1h} } = \ \sup _ {\Gamma _ {1h} } | \phi _ {1h} ( x) | $$

and

$$ l _ {1h} u _ {h} ( x) = u _ {h} ( 0 , x ) ,\ p _ {1h} \phi _ {1} = \ \{ \phi _ {1} \} _ {1h} . $$

For an approximation of the boundary condition one chooses the grid $ \Gamma _ {2h} $ consisting of the points $ ( n \tau , 0 ) $, $ n = 1 , 2 , . . . $; the space $ \Phi _ {2h} $ of grid functions $ \phi _ {2h} ( t) $ with the compatible norm

$$ \| \phi _ {2h} \| _ {\Phi _ {2h} } = \ \sup _ {\Gamma _ {2h} } \ | \phi _ {2h} ( t) | ; $$

the operator $ l _ {2h} $ defined by

$$ l _ {2h} u _ {h} ( t) = \frac{1}{2} [ u _ {h} ( t , h ) - u _ {h} ( t , - h ) ] ; $$

and the operator $ p _ {2h} \phi _ {2} = \{ \phi _ {2} \} _ {2h} $. Then

$$ z ( h) = \| L _ {h} [ u ] _ {h} - \{ L u \} _ {h} \| _ {F _ {h} } + \| l _ {2h} [ u ] _ {h} - \{ l _ {2} u \} _ {2h} \| _ {\Phi _ {2h} } $$

and $ z ( h) \rightarrow 0 $ as $ h \rightarrow 0 $. If for $ W $ one takes the more restricted set of functions $ u \in U $ having, apart from those mentioned above, bounded partial derivatives $ \partial ^ {2} u / \partial t ^ {2} $ and $ \partial ^ {4} u / \partial x ^ {4} $ in $ D $, then for $ \tau = h $ one obtains an approximation of order two.

The grid problem

$$ \tag{3'} L _ {h} u _ {h} = f _ {h} , $$

$$ \tag{4'} l _ {ih} u _ {h} = \phi _ {ih} \ \mathop{\rm on} \Gamma _ {ih} $$

is said to be well-posed if for sufficiently small $ h $ the following conditions are satisfied: The problem is solvable for any given functions $ f _ {h} \in F _ {h} $, $ \phi _ {ih} \in \Phi _ {ih} $ and there is a function $ \psi ( \lambda ) $, independent of $ h $, that tends to zero together with $ \lambda $ and is such that

$$ \tag{5 } \| u _ {h} \| _ {U _ {h} } \leq \psi \left ( \| L _ {h} u _ {h} \| _ {F _ {h} } + \sum _ { i= } 1 ^ { s } \| l _ {ih} u _ {h} \| _ {\Phi _ {ih} } \right ) $$

for any function $ u _ {h} \in U _ {h} $. The solution of a well-posed grid problem depends continuously on the given functions, uniformly with respect to $ h $. The fact that a grid problem is well-posed is a necessary condition for small sensitivity of its solution to rounding in the course of computation. If $ u $ is the solution of (1), (2) and if the problem (3'}), (4'}) is well-posed, then

$$ \tag{6 } \| u _ {h} - [ u ] _ {h} \| _ {U _ {h} } \leq \psi ( \rho ( h ) ) . $$

Moreover, if the problem (3), (4) approximates the problem (1), (2) on the solution $ u $, then

$$ \| u _ {h} - [ u ] _ {h} \| _ {U _ {h} } \rightarrow 0 $$

as $ h \rightarrow 0 $. If the condition for compatibility of the norms, the condition $ z ( h) \rightarrow 0 $ and condition (5) are satisfied, then by putting $ u _ {h} = [ u ] _ {h} $ and by a limit transition in (5) as $ h \rightarrow 0 $, one obtains the inequality

$$ \tag{7 } \| u \| _ {U} \leq \psi \left ( \| L u \| _ {F} + \sum _ { i= } 1 ^ { s } \| l _ {i} u \| _ {\Phi _ {i} } \right ) , $$

which is true for any function $ u \in W $. Thus, to investigate whether boundary value problems of the form (1), (2) are well-posed one can proceed as follows: first obtain an estimate (5), and from it estimate (7). By means of estimates of the type (5) one can often prove the existence of a solution $ u $ of the problem (1), (2) as the limit of the grid approximations $ u _ {h} $ as $ h \rightarrow 0 $.

For a complete solution of the problem of finding an approximate solution of the problem (1), (2) it is necessary to construct an exact or approximate method for finding a solution of the approximation (3), (4) that has the property of stability to rounding. To investigate stability it is useful to use the concept of closure of a computational algorithm. In the given example it is advisable to solve the equation by the double-sweep method with respect to the variable $ x $.

The estimate (6) of the error of the grid method has the disadvantage that in the expression of the function $ \psi ( \lambda ) $ derivatives of the exact solution $ u $ usually occur. In some cases one can estimate these derivatives a priori (that is, before giving a solution of the problem), but such estimates usually turn out to be coarse. There are somewhat more accurate estimates in which the derivatives are replaced by difference quotients of the approximate solution $ u _ {h} $. A practical estimation of the error of grid methods is mainly carried out by means of repeated solutions of the problem (3), (4) with various $ h $ and subsequent isolation of the main part of the error of the form $ h ^ \gamma [ z ] _ {h} $, where $ \gamma $ is the known order of smallness of the error. Thus, if one has the asymptotic relation

$$ \| [ u ] _ {h} - u _ {h} - h ^ \gamma [ z ] _ {h} \| _ {U _ {h} } = \ o ( h ^ \gamma ) , $$

then

$$ [ z ] _ {h} \sim \frac{u _ {h} - u _ {H} }{H ^ \gamma - h ^ \gamma } . $$

One can sometimes obtain equations for $ z $ that contain derivatives of the exact solution $ u $. Then one can solve them numerically on a coarser grid after solving the original problem, obtain the main error term and add it to the approximate solution $ u _ {h} $, thereby refining it.

In some cases, using a special choice of the coordinate functions in variation or projection methods for solving the problem (1), (2) one obtains equations of the form (3), (4) that ensure convergence not only to the classical but also to the generalized solution. This method of constructing approximations, sometimes called the finite-element method, allows great freedom in the choice of the grid. The possibility of appropriately arranging the nodes makes it possible to attain the required accuracy with fewer nodes of the grid.

References

[1]	I. [I. Babushka] Babuška, M. Práger, E. Vitásek, "Numerical processes in differential equations" , SNTL (1966)
[2]	N.S. Bakhvalov, "Summary of the course "Foundations of computational mathematics" " , 4 , Moscow (1968) (In Russian)
[3]	N.S. Bakhvalov, "Numerical methods: analysis, algebra, ordinary differential equations" , MIR (1977) (Translated from Russian)
[4]	R.S. Varga, "Functional analysis and approximation theory in numerical analysis" , Reg. Conf. Ser. Appl. Math. , 3 , SIAM (1971)
[5]	M.K. Gavurin, "Lectures on computing methods" , Moscow (1971) (In Russian)
[6]	S.K. Godunov, V.S. Ryaben'kii, "The theory of difference schemes" , North-Holland (1964) (Translated from Russian)
[7]	E.G. D'yakonov, "Difference methods for solving boundary value problems" , 1–2 , Moscow (1971–1972) (In Russian)
[8]	L.V. Kantorovich, G.P. Akilov, "Functional analysis" , Pergamon (1982) (Translated from Russian)
[9]	G.I. Marchuk, "Methods of numerical mathematics" , Springer (1975) (Translated from Russian)
[10]	S.G. Mikhlin, "Approximate methods for solution of differential and integral equations" , American Elsevier (1967) (Translated from Russian)
[11]	R.D. Richtmeyer, K.W. Morton, "Difference methods for initial-value problems" , Wiley (1967)
[12]	V.S. Ryaben'kii, A.F. Filippov, "Ueber die Stabilität von Differenzengleichungen" , Deutsch. Verlag Wissenschaft. (1960) (Translated from Russian)
[13]	A.A. Samarskii, "Theorie der Differenzverfahren" , Akad. Verlagsgesell. Geest u. Portig K.-D. (1984) (Translated from Russian)
[14]	A.A. Samarskii, E.S. Nikolaev, "Numerical methods for grid equations" , 1–2 , Birkhäuser (1989) (Translated from Russian)

References

[a1]	U.M. Ascher, R.M.M. Mattheij, R.D. Russell, "Numerical solution for boundary value problems for ordinary differential equations" , Prentice-Hall (1988)
[a2]	H.B. Keller, "Numerical solution of two point boundary value problems" , SIAM (1976)