A METHOD
OF
SOLVING
BOUNDARY
VALUE
PROBLEMS*
S. A. PIYAVSKII Kuibyshev 24 November
(Received
1969)
THE majority of methods of solving boundary value problems for ordinary differential equations used in engineering practice require one to choose a “good” initial approximation. The method suggested in the present paper requires the assignment of no initial approximation. Ideas of extension with a parameter are used in its development. In Section 1 the method is demonstrated by an example of two differential equations, in Section 2 it is considered for an ndimensional set of equations and in Section 3 its application to linear boundary value problems is investigated in detail. 1. We consider the following boundary value problem: to find functions y”(t), p2(t) on IO,T’l, satisfying the set of two ordinary differential equations ti’ = fit& Y’, Y’),
ti2 = f2(4 Yi, Y2)
(1.1)
and the boundary conditions Y’(O)
(the functions
=
4
y*(P)
=
B
(1.2)
fi(t, yl, y”), f2(t, Y’, y”) are continuous and twice differentiable).
In Fig. 1 the required solution is shown in the phase plane by a curve with the properties: (a) is satisfies (1.1); (b) its initial and final points belong to the initial and final sets (1.2); (c) the time of motion from the initial to the final point T* (for convenience we shall consider the time to be the independent variable). Our method of finding this solution of the problem consists of including the required curve in a family C of curves for which the conditions (a) and (b), but *Zh. uychisl. Mat. mat. Fiz., 10,
4, 10311036,
290
1970.
291
S. A. Piyavskii
Y4
In. Fin.
FIG.
not (c), are satisfied;
1.
in other words the curves of the family differ in their
time of motion from the initial to the final point, Let us consider two close curves (y”‘(t), y”z(t)) and (y’(t), y”(t)), family C whose times of motion are T and T + AZ’
respectively.
a connection between their initial and final values.
Introducing
Ay’
=
y’_
Ayz
y”‘,
=:
y” _
in the
We establish
~“2,
we can write the set of equations in the variations Ad’ = filll(t, y”‘, y”“).Ay’+ fiy~(t, @‘, Ari2 =
fzul(4
.?)AY~,
(1.3)
y”“)Ay’ +
V, f21/2(6
y”‘,
y”2)A~“,
from which AY’(T)
=
Gif(T, A, y"02)Ayoi$ Gz(T,
Ay2(T)
=
Gzi
Gzi(T,
Here Gis(t),
i, s =
i, 2,
are
~02)Ayo2,
A, (7’,
A,
y”08)Ayo’ +
A,
y”02)Ayo2.
(1.4)
the elements of the matrix G of the fundamental
set of solutions of (1.3) for G(0) = E.
(I.9
The solutions of (1.3) depend on the functions g’(f), and p(t), in the coeificients of (1.3), and these sre obtained by integrating (1.1) with initial conditions
292
A method of solving
In tegration
boundary value problems
(1.7)
\
Intwg \
$I \
\ v
FIG.
2.
y”‘(o) = A (since the curve belongs to the family 0, and p(o) = y%. This fact is reflected in the writing of the arguments for the functions G,, in the formula (1.4).
and
The final values of the increments of the phase coordinates Ay’ (T are connected with the initial values by the relations Ayz(T + AT) Ay’(T
+ AZ’) =
Ayil =
Ay’(T)
+ f’(T,
j7si, B)AT,
AT) =
Ayiz =
Ay2(T)
+ f”(T,
gi’,
+ AT)
(1.6) Ay2(T+
@AT
(here y”t*= y”‘(T)., B = F(T)). If both curves belong to the family C, then ,Ay,’ =
0,
BAyi =
(L2a)
0,
and the system (1.6), (1.2a) gives the increments of the initial and final values of the phase coordinates on the passage between two neighbouring curves of the family C which differ in the time of motion along them (by a quantity AT). On substituting (1.4) in (1.6), using (1.2a) and proceeding to the limit as ATO, we obtain $
=
G22
(T,
G2
4,
CT,
,o~,~
4
yo2) s+
=
fiU’

PV,
yil,B),
vi’,B).
293
S. A. Piyauskii
The system (1.7) is a system of two ordinary differential equations with two unknowns y,‘, y,?. Its initial conditions are easily written for 7’ = 0, which represents
the trajectory
Over this trajectory
of the family
the initial
C with zero length
and final values
(the point 0 in Fig.
of the phase
coordinates
1).
are the
same and therefore yl’(O)
If we integrate the solution
(1.7) with initial of the original
= A,
conditions
boundary
(1.8)
yo2(0) = B.
value
(1.8) up to the instant problem (l.l),
(1.2).
(A, y$(T*))
of final
and the pair (yl’(T”), B) solution of the Cauchy
gives the complete set of initial conditions conditions for finding y’(t), 8”(t) by the simple
problem for equations
(1.1).
Note 1. The writing in exceptional
cases.
T*, we obtain In fact the pair
and integration
of (1.7) are analytically
In the numerical
integration
possible
of (1.7) to evaluate
only
the
to integrate system (1.1) with functions Giz(l’, A, yoz). i = 1, 2, it is necessary initial conditions y’(O) = A, ~~(0) = YO* and (1.3) and (1.5) up to the instant
T. This process abscissa
is shown graphically
in Fig. 2: to “make a step”
axis we must “go along the segment”
that “the passage
along these
of them is a solution
segments”
of the original
10,
~1 of
is not always
boundary
value
along the
the ordinate superfluous
axis. work.
Note Each
but with time T
problem,
different from T*, so that where it is required to construct not one solution, their whole class with variation of time, there are generally no superfluous calculations
but
in the method.
Note 2. The question
of the existence
of the solution
of (1.7) in [0, T’]
is rather complicated but is necessary for the development. Equations (1.7) may even not be written in the normal Cauchy form if Gz, vanishes in IO,T*l its solution may not exist. having 2. problem:
This
corresponds
to the original
boundary
value
problem
no solution. We now turn to the solution the vector
function
of the following
y(t) = (.?(t), . . . , p(t))
common boundary must satisfy
d = f(C Y), k conditions
at the instant
(2.1)
t = 0,
gi(YO)= 0, and n  k conditions
value
the system
at the instant
i = 1, . . . , k, t = T”,
Yo = Y(O),
(2.2)
A method
294
of solving
boundary
j = I,...)
Ql(Yi) = 0,
value problems
Yl = Y(T').
nk,
(2.3)
Here the vector function f(t, y) and the functions gi (yo), qj(yt) are continuous and twice differentiable. In addition if we conditionally consider that yo = Yl = 2,
(2.4)
the conditions (2.2) and (2.3) lead to the system gf (2) = 0,
CljCz)= 0
(2.5)
(here are n equations with n unknowns which are the components of the vector z), which has a unique solution z*. As in the previous section we shall consider the set of equations in the variations A$ = fg (t, y”)AY,
fll(4 !a =
vi,,
i, j = 1, . . . , n,
(4 !a),
(2.6)
connecting the increments of the phase coordinates on passage from the curve ‘j;(t) of the family C to its neighbour of the same family, along which the time of motion is greater by AT. Using the matrix G of the fundamental solutions of (2.6) with the condition G(0) = (E
E
(2.7)
is a unit matrix), we obtain AY(T) =
G(T, Yo)Ayo
and as for (1.5) AYE = AY (T +
4
= G(T, fTo)Ayo+ f(T, y"i)AT.
(2.8)
Adding to (2.8) the relations gi,Ayo = 0,
QjvAY1 =O,
i = 1,. . . , k,
j = 1,. . . , n 
giu
k,
=
(gtg4..
Qjv = (Qjv4 *t Q/vn)
and proceeding in (2.8), (2.9) to the limit as AT+O, a set of equations which we write as follows. We introduce the column vector x = (xl,...,
.,gfyn)t
Xzn),
(2.9)
as for (1.7) we can obtain
(2.10)
295
S. A. Piyavskii
the 2n x 2n matrix
k(
gll
Il._, 0
M=
(2.11)
G(T,~)
E Ic Tl
lI
where .!Tv= qlr =
and the matrix of xs instead
1,. . . , k,
i =
(give)7
j =
(4jy8)7
I,....,
s =
s =
nk,
II,
vector (2.12)
f(T,x!)* II
?L
where f(T, X) is obtained
1, . ..)
from the matrix G (t, ya) by the substitution
G(T, x) is obtained
of yOs, s = 1, . . . , n, and the column N = (0, . . , 0,
1,...,
1,. . . , n,
t by 7’ and ys by ~l”+~, s =
from f(t, y) by replacing
n. Then the set of 2n equations
M(T, q;
Note.
from (2.8)
obtained
and (2.9) assumes
= N(T, x).
(2.13)
(2.13) is of order 2n, but its n first
The system
the form
,gi(x) = 0,
integrals
i=l
are substituted for yes) and , ...I k are known (in (2.2) x8, s = 1, . . . , n), (in (2.3) x”+‘, s = 1,. .., .n), qj (x) = 0, j = 1,. . . , n  k are substituted for yl”) so that the order is easily The set of initial
to n.
reduced
conditions
for (2.13) is of the form
x(O)= I I
By integrating the complete n components
(2.13) with initial solution of
Iz* II
condition
of the original
x(2’*)
Z’
(2.14) up to the instant
boundary
give all the initial
j’(Y)
Both the notes
value
values
P(O) = X’(F), and the last n components
.(2.14)
problem,
because
of the phase
s = I,...,
n,
T*, we obtain the first
coordinates (2.15)
all the final values = c+*(P),
at the end of the previous
s = 1,. . . , n.
section
are applicable
(2.16) here: (2.13) is
296
A method of solving
integrated
along the abscissa
In the second
note there
problems
axis and (2.1) and (2.6) along the ordinate
is the question
may not admit of the description
of the degeneracy
of (2.13) in the Cauchy
dx = dT
3.
boundar.v value
M‘(T,
axis.
of the matrix hl which form
x)N(T, x).
We now apply the above method to solve
(2.13a)
the linear
boundary
value
problem (3.1)
ti = A (t)y + B(t), for t = 0 g(yo)
=
aye +
b =
(3.2)
0,
and for t = T” q(yi)
Here y(t) is an ndimensional
vector
of these
a = (aiS), C=
We assume
(Cjs),
d =
matrices
(3.3)
0.
function, k, s
are continuous i =
1, 2.
. . , k,
d =
s =
1, 2,
. . ., n.
(c,),
1I...?
4
functions,
b = (bi),
j=l,
2 ,***I
n 
k,
that the n Y.n matrix D=
is nondegenerate
solution
The system
a
IIII c
and then the system g(z) = 0,
has a unique
CYI +
B(r) = (bk(t)),
A(t) = (%(t)), and the elements
=
q(z)
=
0
z*.
(2.6) assumes
the form $1 = A(l)Ay
(3.4)
and is independent of the function .?t). Therefore the elements of the matrix M are functions of T alone and there is no necessity to integrate (2.1) at each step of the integration of (2.13). Equation (2.13) now takes the form
297
S. A. Piyavskii
M(T);T= A(T)xf
B(T),
(3.5)
where
and it is sufficient to T* to obtain
can here be reduced
of the boundary
(3.5) and nz equations value
by half (see note to Section
We now transform independently
2n equations
to integrate
the solution
problem.
(3.4) from 0
The order of (3.5)
2).
(3.5) in order to write it in the normal Cauchy
of the degeneracy
of the matrix M(T).
We introduce
form the new
variables U r= M(T)%.
(3.6)
For them du = dT
We also introduce
dM(T) ,,x+A(T)x+B(T).
(3.7)
F=MEE.
(3.8)
the matrix
It is easy to verify directly
that
dM(T) = A(T)F.
(3.9)
dT
Then (3.7) assumes
the form du dr = A(T) (F + E)x + B(T) = A(T)u
We have written ndimensional
(3.5) in the normal Cauchy vector
form.
+ B(T).
(3.10)
Moreover if we introduce
v = (vi, . . . , u”), us = ZZ*+~,s = 1, . . . , n,
it follows
the from
(2.10) that us = const,
which is the same as the original boundary
value
problem
system
has been reduced
dv/dT= (3.1).
Au/B,
(3.11)
Thus the solution
to the following
procedure:
of the linear
298
A method of solving
(11 z* is calculated
boundary value problems
and the vector
~(O)=~(O)~(O),
x(o) =
z' 2' % Ii //
M(O)=
determined; (2) the vector uO, UO’= z+“(O),
S = 1, * *. , n;
is determined; (3) the original set of n equations is integrated up to the instant !I’* with the initial condition v(O) = v. and the n* equations in the variations (3.4) with initial conditions (2.7) (to calculate the elements of the matrix M); (4) if the matrix ~(~*) is nondegenera~, the solution of the original boundary value problem can be obtained at the instant T* by the formulae WY
= M‘(T’)u(T*),
u”+~(T*)
= va(T’),
u*(T*) = us(O), s=
1,,..,n.
The author wishes to thank V. V. Korzhenkov for his useful advice and K. L. Polyanskii for collaborating in the analysis of the linear problem. Tr~sla~
by H. F. Cleaves