On Non-linear Boundary Value Problems and Parametrisation at Multiple Nodes

We show how a suitable interval division and parametrisation technique can help to essentially improve the convergence conditions of the successive approximations for solutions of systems of non-linear ordinary differential equations under non-local boundary conditions. The application of the technique is shown on an example of a problem with non-linear integral boundary conditions involving values of the unknown function and its derivative.


Introduction
Recently, boundary value problems with non-local conditions for non-linear differential equations have attracted much attention (see, e. g., the editorial note [1] and the rest of the issue for extensive references).Problems with non-local boundary conditions are usually treated by using equivalent reformulation as a suitable fixed point or coincidence equation, for which purpose, as a rule, one uses Green's operator of a linearised problem.The process of approximation of the solution based directly on this kind of representations, however, may be quite complicated.
A reasonably efficient way to deal with this kind of problems is provided by methods of numerical-analytic type (see, e. g., [3]).Since convergence conditions often involve terms proportional to the length of the time interval, the conditions needed for the applicability of this type of methods can be significantly weakened if one constructs the scheme using a suitable interval division.It turns out that, by introducing a single intermediate point, one can weaken the convergence conditions by half at the cost of one more variable in the parameter list (see [5][6][7]).In [7], where mainly the periodic problem is considered, we also note that it is possible to consider multiple interval divisions.A scheme of this kind, which is applicable to the case of general boundary conditions, is constructed in the present note.
The approach that we are going to discuss is based on a suitable parametrisation, so that the values of approximations to a solution are monitored at multiple time instants.In this way, it can be regarded as an efficient alternative to the multiple shooting [2,10] and may be well applicable also in the cases where shooting procedures fail.The latter may happen either because of the complicated character of the boundary conditions (according to our knowledge, the currently available shooting schemes are designed for the cases where the boundary conditions are local two-point) or, more importantly, due to the failure to satisfy the basic assumptions needed to apply the method.Indeed, one may note that shooting methods require the existence of sufficiently many derivatives of the non-linearity (in particular, because Newton-like methods are commonly used to solve the corresponding numerical equations, see, e.g., [10, p. 516] or [11, p. 375]).Furthermore, in order to carry out shooting, one has to be sure that the initial value problem for the differential equation in question has always a unique solution defined on the entire given time interval.The smoothness of the non-linearity alone is insufficient: consider, e. g., u = u 2 on [a, b] with u(a) = 1/(λ − a), where a < λ < b; then the solution u(t) = 1/(λ − t) is undefined at t = λ.Last, but not least, the existence of a solution is usually assumed a priori when applying shooting methods.In contrast to this, the approach that we suggest here, in many cases, allows one to prove the solvability of the problem in a rigorous way (see, e. g., [7,9]).
Here, we study the non-linear boundary value problem where f : [a, b] × R n → R n , d ∈ R n is a given vector, and φ is a vector functional on the space of absolutely continuous functions (generally speaking, non-linear).Following the idea used in numerical methods for approximate solution of initial value problems for ordinary differential equations, let us fix a natural N and choose N + 1 grid points where 3), however, may allow one to pose better conditions on the non-linearity in the corresponding region.The aim of this note is to present an approach to problems of type (1.1)-(1.2) which is similar in principle to [7] and is also based on reductions to certain simpler problems with unknown parameters.The auxiliary two-point problems are constructed here with multiple interval divisions, which leads one to convergence conditions significantly weaker than in the case of a single intermediate point.Here, in contrast to the case of linear two-point conditions discussed in [6,7], the exact fulfilment of the boundary condition for approximations is not guaranteed any more (of course, the boundary condition is satisfied exactly in the limit).The advantage is, however, that many different types of boundary conditions can be thus handled in a unified way, the specific properties of the problem being transferred to the determining equations.It seems that, in the case of general boundary value problems, interval division for approximations constructed analytically is employed here for the first time.

Notation
We fix an n ∈ N and a bounded closed set D ⊂ R n .For vectors x = col(x 1 , . . ., x n ) ∈ R n the obvious notation |x| = col(|x 1 | , . . . ,|x n |) is used and the inequalities between vectors are understood componentwise.The same convention is adopted implicitly for operations like "max" and "min".
1 n and 0 n are, respectively, the unit and zero matrices of dimension n. r(K) is the maximal, in modulus, eigenvalue of a matrix K.
Similarly, the componentwise -neighbourhood of a set Ω ⊂ R n is defined as For given two bounded connected sets D 0 ⊂ R n and D 1 ⊂ R n , introduce the set for any {τ 1 , and max t∈[τ,τ+l] α 1 (t, τ, l) = l/2.
3 Parametrisation and auxiliary problems

Parameter sets
Let us fix certain closed bounded sets and focus on the absolutely continuous solutions u of problem (1.1)-(1.2) whose values at nodes (1.3) lie in the corresponding sets (3.1), i. e., the solutions u such that Given sets (3.1), we introduce the sets and, for any non-negative vector , put Recall that, according to (2.3), (2.4), D k−1,k is the set of all possible straight line segments joining points of D k−1 with points of D k , whereas Ω k ( ) is the componentwise -neighbourhood of D k−1,k .

Freezing
The idea that we are going to use suggests to replace the original non-local problem (1.1)-(1.2) by a suitable family of model boundary value problems with simpler boundary conditions (see, e. g., [8,9].Let us do this in the following way.Consider the vectors where N is the number of nodes from (1.3).These vectors will be regarded as unknown parameters whose values are to be determined.Let us "freeze" the values of u at the nodes (1.3) by formally putting and consider the restrictions of equation (1.1) to each of the subintervals of the division: Then, in a natural way, we have For any fixed k = 1, 2, . . ., N, relations (3.7), (3.8) can be regarded formally as an overdetermined boundary value problem with two-point boundary conditions containing unknown parameters z (k−1) and z (k) .This leads one to a kind of reduction principle where, instead of the original equation (1.1), one considers the parametrised problems (3.7), (3.8) and tries to determine the appropriate value of z (0) , z (1) , . . ., z (N) .
Due to the form of the boundary condition (3.8), it is natural to apply to (3.7), (3.8) techniques similar to those used in [7] for two-point problems.This is done in Section 4.2 below, where the successive approximations x are constructed.Note that the differential equation (3.7) is considered on an interval of length h k (see (1.3)).

Interval division and successive approximations 4.1 Assumptions
Let us fix the sets D k ⊂ R n , k = 0, 1, . . ., N, from (3.1).We make the following assumptions.
We suppose that f is Lipschitzian, in the space variable, on the sets Ω k ( (k) ), k = 1, 2, . . ., N. Namely, Assumption 4.2.There exist non-negative matrices Finally, we assume in the sequel that the matrices Assumptions 4.1 and 4.2, together with condition (4.3), are used to prove the applicability of the techniques described below.They mean essentially that the non-linearities in the equation are Lipschitzian on sufficiently large domains ( (k) satisfies inequality (4.1)) with sufficiently small constants (condition (4.3)).It should be noted, however, that (4.1) and (4.3) are both satisfied if the number N of nodes in (1.3) is large enough.Thus, the basic and, in fact, the only restrictive assumption in this note is that f is Lipschitzian on a bounded set.

Convergence of successive approximations
It turns out that the sequences {x . ., N given by (4.4) and (4.5) are helpful for the investigation of solutions of the given problem (1.1)-(1.2).Theorem 5.1.Let Assumptions 4.1 and 4.2 hold and, moreover, the corresponding matrices K 1 , K 2 , . . . ,K N satisfy condition (4.3).Then, for any (z (0) , z (1) , . . ., 2. The limit function (5.1) satisfies the conditions (5.2) ) is the unique absolutely continuous solution of the integral equation Proof.The proof is carried out similarly to that of [8,Theorem 3].Let us fix arbitrary vectors z (i) ∈ D i , i = 0, 1, . . ., N, and a number k ∈ {1, 2, . . ., N}.We first show that, under the conditions assumed, for any m ≥ 0. Indeed, the validity of (5.6) for m = 0 is an immediate consequence of (4.4).Let us put where In view of (4.1), this means that x (k) (5.6) holds for m = 1.Using this and arguing by induction with the help of Lemma 2.2, we easily establish that for k = 1, 2, . . ., N and m ≥ 2. Therefore, (5.6) holds for any m ≥ 0.
Theorem 5.1 implies, in particular, that one can introduce the function Then it follows immediately from (5.3) that the following statement holds.

Limit functions and determining equations
It is natural to expect that the limit functions x . ., N, of iterations (4.5) may help one to state general criteria of solvability of problem (1.1), (1.2).Such criteria can be formulated in terms of the respective functions ∆ (k) : D k−1 × D k → R n , k = 0, 1, . . ., N, given by equalities (5.14) that provide such a conclusion.Indeed, Theorem 5.1 ensures that, under the conditions assumed, the functions x we obtain a function u ∞ (•, z (0) , z (1) , . . ., z (N) ) : [a, b] → R n , which is well defined for all the values z (k) ∈ D k , k = 0, 1, . . ., N. This function is obviously continuous because, at the points t = t k , we have x for all k = 1, 2, . . ., N − 1. Equalities (6.2) follow immediately from the fact that the function x ) is a solution of equation ( 5.

Approximate determining equations
Although Theorem 6.1 provides a complete theoretical answer to the question on the construction of a solution of problem (1.1)-(1.2),its application faces complications since it is difficult to find the limit function (5.1) and, as a consequence, the functions (6.3) are usually unknown explicitly.The complication can be overcome if we replace the unknown limit 1) , z (k) ) of form (4.5) for a fixed m and put We see that (7.1) is an approximate version of the unknown function (6.1).Its values can be found explicitly for all t ∈ [a, b] and z (k) ∈ D k , k = 0, 1, 2, . . ., N.
Note that, unlike system (6.3), the mth approximate determining system (7.2) contains only terms involving the functions x 1) , z (j) ), j = 1, 2, . . ., N, and, thus, known explicitly.Let ( z(0) , z(1) , . . ., z(N) ) be a solution of the approximate determining system (7.2) for a certain value of m.Then the function is natural to be regarded as the mth approximation to a solution of the boundary value problem (1.1)-(1.2).In particular, it follows from (5.17) that The existence of a solution can be analysed based on the approximate determining equations (7.2) similarly to [3,9], this topic is not considered here.In relation to estimate (7.3) one may note that, according to Theorem 6.1, the solution necessarily has form (6.4) with certain values of z (k) , k = 0, 1, . . ., N. Thus, we have z (k) ≈ z(k) , k = 0, 1, . . ., N, and, therefore,

Example
Let us demonstrate the approach described above on a model example.Consider the system of differential equations with the integral boundary conditions Clearly, problem (8.1), (8.2) is a particular case of (1.1)-(1.2) with a = 0, b = 1.9, d ≈ col (0.21, 0.428), x → φ(x) := col 1.9 0 Let us choose the grid (1.3) with N = 3 and the nodes Then, obviously, According to (3.5), the scheme will depend on four two-dimensional vector parameters z (k) , 0 ≤ k ≤ 3; their meaning is explained by Table 8.2. Variable Value it approximates x(0) x( Table 8.2:The meaning of the parameters in the example. The number of the solutions of the algebraic determining system (7.2) coincides with the number of the solutions of the given problem.Different solutions have to be detected by changing appropriately the initial domains D k , 0 ≤ k ≤ 3. Let us carry out several steps of iteration with two different choices of the initial domains and the radii of neighbourhoods.

First solution
Let us choose the sets D k , 0 ≤ k ≤ 3, as follows: This choice can be justified by the fact that the zeroth approximate determining system (i.e., (7.2) with m = 0) has roots lying in these sets (see the first column in Table 8.3).Furthermore, sets (8.6) contain the corresponding parts of the graph of the zeroth approximation.The graphs of the components of the latter function, which, according to (4.4), have the form of broken lines, are shown on Figure 8.1.
Using (3.3), we find that the corresponding sets D k−1,k , 0 ≤ k ≤ 3, have the form  In order to construct suitable sets on which Assumptions 4.1 and 4.2 will be verified, we need to choose vectors (1) , 0 ≤ k ≤ 3. Let us put, for example, (1) := col (0.2, 0.3), (2) := col (0.1, 0.4).( Then, according to formula (3.4), the corresponding sets Ω k ( (1) ), 0 ≤ k ≤ 3, have the form A direct computation shows that the Lipschitz condition (4.2) for the right-hand side terms of (8.1) holds in Ω 1 ( (1) ), Ω 2 ( (2) ), Ω 3 ( (3) ), respectively, with the matrices Relations (8.10) show that matrices (8.9) satisfy conditions (4.3) with the step sizes (8.5).Furthermore, in view of (8.5), (8.7), and (8.8), we have 1) , 3) , which means that vectors (8.7) satisfy conditions (4.1) of Assumption 4.1.Thus, we see that that all the conditions of Theorem 5.1 are fulfilled, and the sequences of functions (4.5) for this example are convergent.Using Maple 14 for constructing the iterations and solving the approximate determining equations (7.2) for m = 0, 1, 2, 9, we obtain the numerical results shown in Table 8.3.We may note at this point that, at nodes (8.4), the pair of functions (8.3), which, as has been indicated, is a solution of problem (8.1), (8.2), has the values listed in Table 8.1.Comparing Tables 8.3 and 8.1, we find enough evidence to claim that the results of computation with the present choice of initial domains correspond to solution (8.3).This is further confirmed when we put the components of this function and the first approximation (m = 1) on the same plot (see Figure 8.2).The graphs of higher approximations (we have carried out computations up to m = 9) practically coincide with one another and there is no way to distinguish them in the given resolution.Considering the difference between the approximation and solution (8.3), e. g., for m = 1, we see that the maximal error is about 6 • 10 −4 (see Figure 8.3).All this, together with Tables 8.3 and 8.1, demonstrates a rather high quality of approximation.
Carrying out computations, we see that the approximate determining systems (7.2), along with the solution found in Section 8.1 in sets (8.8) (see Table 8.3), has another solution in sets (8.12).The corresponding approximate values of parameters at several steps of iteration (m = 0, 1, 2, 7, 9) are presented in Table 8. 4. In particular, we see that, as in Section 8.1, the piecewise linear zeroth approximation provides a useful hint as to where the solution should be looked for (see the first column of Table 8.4).The graphs of three approximations to this solution (m = 0, 1, 9) are shown on Figure 8.4.The residual obtained as a result of substitution of the ninth approximation into the given differential system (8.1) is of order of 10 −8 .