Non-linear boundary value problem, numerical methods

Methods replacing a boundary value problem by a discrete problem (see Linear boundary value problem, numerical methods and Non-linear equation, numerical methods). In many cases, especially in the discussion of boundary value problems for systems of ordinary differential equations, the description of numerical methods usually proceeds without indication of a discretization of the original problem, but such a discretization is nevertheless implicit and is realized relative to a known model.

Suppose, for example, that a two-point boundary value problem is to be discussed for a system of ordinary differential equations

$$ \tag{1 } y ^ \prime = f ( x , y ) \ \ \textrm{ for } 0 < x < X , $$

$$ \tag{2 } B ( y ( 0) ) = 0 ,\ D ( y ( X) ) = 0 , $$

where

$$ y \equiv ( y _ {1} \dots y _ {l} ) ^ {T} ,\ \ f \equiv ( f _ {1} \dots f _ {l} ) ^ {T} , $$

$$ B = ( B _ {1} \dots B _ {l-} r ) ^ {T} ,\ D \equiv ( D _ {1} \dots D _ {r} ) ^ {T} , $$

$$ 0 \leq r \leq l . $$

Suppose that the vector $ G ( y) = ( G _ {1} ( y) \dots G _ {r} ( y) ) $ is such that the system of equations

$$ \tag{3 } B ( y) = 0 ,\ \ G ( y) = g $$

uniquely determines a vector $ y = \alpha ( g) $ for every vector $ g \equiv ( g _ {1} \dots g _ {r} ) ^ {T} $; for example, it is often expedient to take

$$ G _ {1} \equiv y _ {l - r + 1 } \dots G _ {r} \equiv y _ {l} . $$

Then the problem (1)–(2) can be reduced to the system (operator equation)

$$ \tag{4 } \psi ( g) = 0 , $$

where $ \psi ( g) = D ( \omega ( g) ) $ and $ \omega ( g) $ is the value at $ x = X $ of the solution of (1) for the initial condition $ y ( 0) = \alpha ( g) $. To obtain the values of $ \psi ( g) $ one has to solve numerically (see [1], [2]) the corresponding system of differential equations, that is, to use a discretization of the original problem. To solve the system (4) one can use various iteration methods for the solution of non-linear equations. Having found a solution $ g ^ {*} $ of (4) one can determine the vector $ y ^ {*} = \alpha ( g ^ {*} ) $ from the system (3) and, by solving (1) for the initial condition $ y ( 0) = y ^ {*} $, obtain a solution of the problem (1)–(2). In a number of methods, apart from the values of $ \psi $ one also uses the value of its derivative, to be found either by integrating the equation

$$ \tag{5 } \eta ^ \prime = \ \frac{\partial f }{\partial y } \eta $$

(the variational equation for (1)) or by means of some formula for numerical differentiation. This way of reducing a boundary value problem to (4) is called the shooting method, as applied to a problem of more general form

$$ y ^ \prime = f ( x , y , g ) \ \ \textrm{ for } 0 < x < X , $$

$$ y ( 0) = \alpha ( g) ,\ D ( y ( X ( g) , g )) = 0, $$

containing a certain vector parameter $ g $ explicitly, in particular, to an eigen value problem; here the function $ \alpha ( g) $ is assumed to be given. However, if the solution of (5) grows rapidly with $ x $, then the error of the numerical integration leads to large errors in the values of $ \psi ( g) $ and its derivative and ultimately to a large error of the resulting solution.

In the method of linearization (see, for example, [1], [3]) approximations to a solution of the problem (1)–(2) are determined as solutions of a sequence of linear equations

$$ \tag{6 } y _ {n+} 1 ^ \prime = \ A _ {n} ( x) ( y _ {n+} 1 - y _ {n} ) + f ( x , y _ {n} ) \ \ \textrm{ for } 0 < x < X $$

with the linear boundary conditions

$$ \tag{7 } \left . \begin{array}{c} B _ {n} ( y _ {n+} 1 ( 0) - y _ {n} ( 0) ) + B ( y _ {n} ( 0) ) = 0 , \\ D _ {n} ( y _ {n+} 1 ( X) - y _ {n} ( X) ) + D ( y _ {n} ( X) ) = 0 , \\ \end{array} \right \} $$

where $ A _ {n} ( x) $ is a square matrix of order $ l $ consisting of functions and $ B _ {n} $ and $ D _ {n} $ are numerical $ (( l- r ) \times n ) $- and $ ( r \times l ) $- matrices, $ n = 0 , 1 , . . . $. The initial approximation $ y _ {0} ( x) $ is assumed to be known. In the case of the Newton–Kantorovich method (cf. Kantorovich process; Newton method) the equations (6) and (7) take the form

$$ \tag{8 } y _ {n+} 1 ^ \prime = \ \frac{\partial f }{\partial y } ( x , y _ {n} ) ( y _ {n+} 1 - y _ {n} ) + f ( x , y _ {n} ( x) ) \ \ \textrm{ for } 0 < x < X , $$

$$ \tag{9 } \left . \begin{array}{c} \frac{\partial B }{\partial y } ( y _ {n} ( 0)) ( y _ {n+} 1 ( 0) - y _ {n} ( 0)) + B ( y _ {n} ( 0)) = 0, \\ \frac{\partial D }{\partial y } ( y _ {n} ( X)) ( y _ {n+} 1 ( X) - y _ {n} ( X) ) + D ( y _ {n} ( X) ) = 0. \\ \end{array} \right \} $$

If $ f ( x , y ) $, $ B ( y) $ and $ D ( y) $ are twice continuously differentiable in $ y $, if the problems (8)–(9) are well-posed, and if $ y _ {0} ( x) $ is sufficiently close to the required solution $ y ( x) $, then

$$ \| y _ {n} - y \| _ {C} \leq K ^ {-} 1 ( K \| y _ {0} - y \| _ {C} ) ^ {2n} , $$

where $ K > 0 $ is a constant and $ \| y \| _ {C} \equiv \max _ {0 \leq x \leq X } | y ( x) | $( see [1]).

The methods mentioned generalize to the more general case of conditions of the form

$$ A ( y ( X _ {0} ) , y ( X _ {1} ) \dots y ( X _ {m} ) ) = 0 , $$

where

$$ A \equiv ( A _ {1} \dots A _ {l} ) ^ {T} ,\ \ 0 = X _ {0} < X _ {1} < \dots < X _ {m} \equiv X . $$

Many important theoretical and applied problems lead to the need of solving non-linear boundary value problems (and related problems) for equations and systems of equations of elliptic type (see, for example, [4]–[8]). For such a class of problems the basic numerical methods are projection methods (projection-grid, variational-difference, finite element) and difference methods (see [7]–[17]). Their constructions are in many respects analogous to those in the corresponding methods for linear boundary value problems (see Laplace equation, numerical methods; Poisson equation, numerical methods). However, in these methods the resulting discrete (grid) analogues to boundary value problems, being systems of non-linear equations

$$ \tag{10 } L _ {N} ( \overline{u}\; _ {N} ) = \overline{f}\; _ {N} , $$

where

$$ \overline{u}\; _ {N} \equiv ( u _ {1} \dots u _ {N} ) ^ {T} \ \ \textrm{ and } \ \ \overline{f}\; _ {N} \equiv ( f _ {1} \dots f _ {N} ) ^ {T} , $$

often involve considerable additional difficulties both with respect to the analysis of the systems themselves and the proximity in one sense or another of their solution to that of the original problem, and with respect to the input of numerical work for specifying the systems (10) and the search for their solutions. The specific nature of the non-linearity of a boundary value problem forces one to pay special attention to the choice of Banach spaces and their grid analogues in which the problem itself and its finite-dimensional approximations lend themselves to analysis (see, for example, [5]–[12], [15]). But the most detailed study, especially of algorithms, has been made only of numerical methods for classes of non-linear problems related in a certain sense to linear problems (weakly non-linear equations, equations with restricted non-linearity), thus making it possible to use the theory of Hilbert spaces and their finite-dimensional analogues (Euclidean spaces).

Let $ H _ {N} \equiv H $ be a Euclidean space of vector-valued functions $ \overline{u}\; _ {N} \equiv u $ with scalar product $ ( u , v ) \equiv \sum _ {i=} 1 ^ {N} u _ {i} v _ {i} $, and let $ l ( H) $ be the set of self-adjoint positive-definite linear operators mapping $ H $ to $ H $; alongside with $ H $ one also uses the Euclidean space $ H _ {B} $, $ B \in l ( H) $, that differs from $ H $ only in the form of scalar product:

$$ ( u , v ) _ {H _ {B} } \equiv \ ( u , v) _ {B} \equiv \ ( B u , v ) ,\ \ \| u \| _ {B} \equiv \ ( B u , u ) ^ {1/2} . $$

Next, let

$$ S _ {B} ( r) \equiv \ \{ {u } : {\| u \| _ {B} \leq r } \} ,\ \ r > 0 . $$

If the operator $ L _ {N} $ is continuous on $ S _ {B} ( r) $ and such that for any $ u $ with $ \| u \| _ {B} = r $,

$$ \tag{11 } ( L _ {N} ( u) - \overline{f}\; _ {N} , u ) \geq 0 , $$

then the system (10) has at least one solution in $ S _ {B} ( r) $. This assertion (see, for example, [9]) is one of the consequences of the classical existence theorems for solutions of non-linear operator equations based on topological principles (see, for example, [14]) and on the monotonicity of the operators (see, for example, [6], [11], [14]). In particular, if $ L _ {N} $ is continuous everywhere and if

$$ \tag{12 } ( L _ {N} ( u) , u ) \geq \ \delta \| u \| _ {B} ^ {2} ,\ \ \delta > 0 , $$

for any $ u $, then the system (10) always has a solution and any one of them belongs to $ S _ {B} ( r) $ with $ r = \delta ^ {-} 1 \| \overline{f}\; _ {N} \| _ {B} ^ {-} 1 $. Moreover, if for any $ u , v , w $, the inequalities

$$ \tag{13 } ( L _ {N} ( u) - L _ {N} ( v) , u- v ) \geq \ \delta _ {0} \| u - v \| _ {B} ^ {2} ,\ \ \delta _ {0} > 0, $$

$$ \tag{14 } | ( L _ {N} ( u) - L _ {N} ( v) , w ) | \leq \ \delta _ {1} \| u - v \| _ {B} \| w \| _ {B} $$

hold, which indicate the strong monotonicity and Lipschitz continuity of $ L _ {N} $, then (10) and its unique solution can be found by the iteration method

$$ \tag{15 } A u ^ {n+} 1 = \ A u ^ {n} - \gamma ( L _ {N} ( u ^ {n} ) - \overline{f}\; _ {N} ) , $$

where $ A \in l ( H) $ and the iteration parameter $ \gamma > 0 $ is determined by constants like $ \delta _ {0} $ and $ \delta _ {1} $, which are obtained when in (13) and (14) the operator $ B $ is replaced by $ A $( see, for example, [8], [9]). Hence, for difference methods estimates of the error result as a consequence of the fact that the system (10) is well-posed and that an approximation of the original boundary value problem by a difference boundary value problem is available. It is frequently expedient to choose the operator $ B $ so that the space $ H _ {B} $ becomes a discrete (grid) analogue of the corresponding Sobolev space (for example, $ W _ {2} ^ {1} ( \Omega ) $ for second-order equations), and so it becomes possible to use either the imbedding theorems themselves (in the case of projection approximations) or their grid analogues (see, for example, [5]–[9]) for the verification of inequalities like (12)–(14). Projection methods and, in particular, projection-grid (finite-element) methods have the merit that (12)–(14) follow from the corresponding inequalities for differential operators, and the task of estimating the error of the method easily reduces to that of approximating the solutions of the original problem by elements of a selected finite-dimensional subspace (see [6]–[9], [11]). Instead of the validity of (13)–(14) in the whole space it suffices that they hold, for example, in some ball in $ H _ {B} $ surrounding the required solution. Sometimes useful inequalities are obtained from (12)–(14) when the scalar products on the left-hand sides are replaced by scalar products in $ H _ {B} $ and the norm $ \| z \| _ {B} $ on the right-hand sides by $ \| B z \| $. In a number of problems, on the basis of a priori estimates of their solutions, one can replace them by equivalent boundary value problems for which inequalities like (12)–(14) are satisfied (see [9]).

When (13) and (14) hold, then the iterative method (15) with $ A = B $ converges for any initial approximation, and if the constants $ \delta _ {0} $ and $ \delta _ {1} $ are independent of $ N $, then one can obtain a solution of (10) with accuracy $ \epsilon $ in $ O ( | \mathop{\rm log} \epsilon | ) $ iterations (usually $ \epsilon \approx N ^ \alpha $, $ \alpha < 0 $, that is, $ 0 < c _ {0} \leq \epsilon N ^ {- \alpha } \leq c _ {1} $). At the expense of a successful construction of a grid method one sometimes manages to choose the operator $ B $ so that the solution of the system (10) to within $ N ^ \alpha $ by means of an iteration method like (15) is achieved by performing only $ O ( N \mathop{\rm log} ^ {2} N ) $ or even $ O ( N \mathop{\rm log} N ) $ arithmetic operations (see Minimization of the labour of calculation). The use of a sequence of refining grids makes it possible to reduce these estimates to $ O ( N \mathop{\rm log} N ) $ and $ O ( N) $ operations and so to obtain asymptotically-optimal numerical methods (relative to the labour involved) for the solution of certain elliptic boundary value problems. Here it is assumed that the numerical residual $ L ( u ^ {n} ) - \overline{f}\; _ {N} $ itself requires at most $ K N $ operations. Nevertheless, for certain problems (with large $ K $) it is expedient to reduce the frequency of such computations. This can be achieved provided one can find a more accurate initial approximation by the use of some linearization of $ L _ {N} $ and by methods of the type of continuation with respect to a parameter (see Non-linear equation, numerical methods). Most often linearization methods in conjunction with methods of continuation for non-linear operators with respect to a parameter (cf. Continuation method (to a parametrized family, for non-linear operators)) are used in the most difficult situations, for example, when (13) fails to hold. Certain eigen value problems (see, for example, [9], [12], [14]) can also be considered as non-linear boundary value problems of elliptic type; the construction of the relevant numerical methods is in many respects analogous to the ones already mentioned.

Some other stationary boundary value problems can be treated numerically in the same manner as well (cf. [8], [9], [17]). The numerical methods for non-stationary non-linear boundary value problems include problems for equations and systems of parabolic type, hyperbolic type, mixed type, and others, although discrete analogues of elliptic operators similar to the ones considered above have also been used; however, for them it is even more important to have approximation methods for the individual time layers.

References

[1]	N.S. Bakhvalov, "Numerical methods: analysis, algebra, ordinary differential equations" , MIR (1977) (Translated from Russian)
[2]	G. Sewell, "The numerical solution of ordinary and partial differential equations" , Acad. Press (1988)
[3]	R.E. Bellman, R.E. Kalaba, "Quasilinearization and nonlinear boundary-value problems" , Elsevier (1965) (Translated from Russian)
[4]	A.A. Berezovskii, "Lectures on the non-linear boundary value problems of mathematical physics" , 1 , Kiev (1976) (In Russian)
[5]	O.A. Ladyzhenskaya, N.N. Ural'tseva, "Linear and quasilinear elliptic equations" , Acad. Press (1968) (Translated from Russian)
[6]	J.-L. Lions, "Quelques méthodes de résolution des problèmes aux limites nonlineaires" , Dunod (1969)
[7]	R. Glowinski, J.-L. Lions, R. Trémolières, "Numerical analysis of variational inequalities" , North-Holland (1981) (Translated from French)
[8]	R. Glowinski, "Numerical methods for nonlinear variational problems" , Springer (1984)
[9]	E.G. D'yakonov, "Minimization of computational work. Asymptotically-optimal algorithms" , Moscow (1989) (In Russian)
[10]	M.M. Karchevskii, A.D. Lyashko, "Difference schemes for non-linear problems of mathematical physics" , Kazan' (1976) (In Russian)
[11]	R.S. Varga, "Functional analysis and approximation theory in numerical analysis" , Reg. Conf. Ser. Appl. Math. , 3 , SIAM (1971)
[12]	G. Strang, J. Fix, "An analyse of the finite element method" , Prentice-Hall (1973)
[13]	J.T. Oden, "Finite elements of nonlinear continua" , McGraw-Hill (1972)
[14]	M.A. Krasnosel'skii, G.M. Vainikko, P.P. Zabreiko, et al., "Approximate solution of operator equations" , Wolters-Noordhoff (1972) (Translated from Russian)
[15]	F. Brezzi, T. Rappaz, P. Raviart, "Finite dimensional approximation of nonlinear problems III" Numer. Math. , 38 (1981) pp. 1–30
[16]	A.A. Amosov, N.S. Bakhvalov, Yu.I. Osipik, "Iterative processes for the problem of stationary heat exchange in a system of absolutely black bodies" J. Comput. Math. Math. Phys. , 20 : 1 (1980) pp. 110–119 Zh. Vychisl. Mat. Mat. Fiz , 20 : 1 (1980) pp. 104–111
[17]	F. Thomasset, "Implementation of finite element methods for Navier–Stokes equations" , Springer (1981)

References

[a1]	A.M. Ascher, "Numerical solution for boundary value problems for ordinary differential equations" , Prentice-Hall (1988)