A method of solving boundary value problems

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 n-dimensional 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, 1031-1036,

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 AT-O, 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').

n-k,

(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 =

n-k,

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 n-dimensional

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 non-degenerate

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=M-EE.

(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 n-dimensional

(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