- Email: [email protected]
Numerical solutions of two-point boundary value problems
0362-546X;79/0301-0395 12.00/O
Nonlinear Analysis. Theory. Merhods & Applicarions, Vol. 3, No. 3, pp. 395-418 Q Pngamon Press Ltd. 1979. Printed in Great Britain.
NUMERICAL
SOLUTIONS OF TWO-POINT VALUE PROBLEMS*
BOUNDARY
M. N. O&JZTBRELI and F. PASCAL Department of Mathematics, University of Alberta, Edmonton, Alberta, Canada, T6G 2Gl and K. V. LEUNG* Department of Computer Science, Concordia University, Montreal, Quebec, Canada, H3G lM8 (Received22
May 197X; revised28 Augusf 197X)
Key words: Nonlinear ordinary differential equations, two-point boundary value problem with mixed conditions, fourth order approximations, iterative method.
1. INTRODUCTION
paper is concerned with the numerical solution of a two-point boundary value problem with mixed boundary conditions for a nonlinear second order ordinary differential equation of the form XEI Y” = f’k Y, Y’), (1.1)
THIS
ay(a) - y’(a) = 6,
(1.2)
PY@) + Y’(b) = Y
where I = [a, b], a, /3, y and 6 given real numbers such that a B 0,
LX+/?>0
B 2 0,
(1.3)
and f(x, y, z) is a real valued function defined in a region 9 c R3, W=IQ9 (l-4) 9 = HY,41 IYJG PI,
I4 G Pz>
in which it satisfies a Lipschitz condition of the form If(&Y,,Y;)
-f(&Yz,Y;)I
G qy,
- Y,l + M/Y;
Y;l
(1.5)
where pt, p2, M and K are constants. The special case in which the function f does not depend explicitly on y’ has been investigated in a number of monographs and papers partially listed in the bibliography. Recently, in [l] and [2], fourth-order tridiagonal finite difference methods for the general problem have been proposed. In [3] the methods of collocation and finite differences have been compared and 0(h4) convergence discussed. In a very recent paper, [4], the existence of solutions for the general nonlinear boundary value problem has been investigated in a Banach space. In the present work we propose an iterative finite difference method which gives globally * This work was partially supported by the National Research Council of Canada by Grant NRC A-4342 to K.V.L. and Grant NRC A-4345 to M.N.O. *&Deceased on I June 1978.
396
M. N. OGUZT~~KELI. F. PASCAL and K. V. LEUNG
0(h4)-convergence over the mesh points of a uniform mesh of width h. We calculated the derivative y’ at the mesh points by a new scheme which differs from those in [l-3 1. We also extended some of the results of [5]. In the next two sections we introduce some preliminary notation and relationships on which our iterative method is based. The iterative procedure is described in Section 4 and is used in Section 5 in connection with the existence and uniqueness of the solutions of the approximating equation which appear in a natural way in our analysis. This section is concluded with our first theorem about approximate solutions and their iterative constructions. In Section 6 we deal with the error analysis in connection with the exact solution of equations (I. l), (1.2) and the approximating iterative solutions considered in the previous section. The results of this section are summarized by our second theorem. In Section 7 a brief explanation of our computer programming is given and the implementation is described in Flowchart I. In the last section we present a short summary of our simulations related with three boundary value problems previously considered in [2] by Stepleman. In Appendix A we establish one of the elementary results used in the second section. 2. PRELIMINARIES
Let N be a positive
odd integer.
Put
h = (b - a)/N,
xn = a + nh
Y’k’= Y’k’(XJ
(2.1)
f, = f(X,, Yn, YA) n = 0, 1, . . . , N).
(k = 0, 1,2;
By virtue of (l.l), (1.2) we have the relationships Yb = UY, -
Y:: = f,?
Using
the following
well known
formula
4
(2.2)
y)N= -py,+s.
of numerical
differentiation
(cf. [6])
+ 2Y” - Yn+i = - g [Y;_ 1 + 10~; + y;+ J + R,,(h)
-y,_i
(2.3)
where R,(h) = O(h6) and the formulae
(A6) and (A8) of Appendix
(1 + ah) Y, - YI = - ;
(n = 1,2,. . , N -
A, (1.1) by the following
[;f,
_Y”_ 1 + 2Yn - Y,+1 = - &_I
-Y,_, with
n = 1,2,. . . , N -
1)
+ 3f,
- $fi]
+ lofn +
system of equations
+ h6 + R,(h)
f,+Il + R,(h)
+ 3fjve1 +
+ (1 + Bh)Y, = - ;[-+r,_,
(2.4)
;&I + b
(2.5)
+ R,(h)
1, where R,(h) = 0(h5),
R,(h)
= 0(h5).
(2.6)
Two-point
boundary
397
value problems
We now put l+ah
- 1
0
2 -1
1
- 1
0
0.. .
0
0
0
0
0.. .
0
0
0
0
.
0
0
0
0
2 -1..
\
A=
I$
0
B=
1
0
2 -1
0
0..
0
0
0..
o-1
0
0
0.. .
0
2 - 1 O-l
1 + bh
(N+
1) x (N+l)
.o
0
0
o\
,(h)
10.
.o
0
0
0
(4
101.
.o
0
0
0
3 -*
110
-1
0
0.
R(h) =
Pp.,1 .l
10
.O
(4
10
1 10 1
.0-t
35
!
(N-1 1) x (N+ 1)
,(h)
Y’=
(N+l)
Then the system (2.5) can be written
x 1
in the matrix
h= AY = - $3F(I: It has been shown in [5] that the inverse A-’
form
Y’) + Hh + of the matrix
R(h).
(2.8)
A is of the form
(2.9) x (N+l)
M. N. O~LT~RELI,
398
F. PASCAL and K. V. LEUNG
with
if 0 < i < j <
(l + ‘ah)[l + (N - j) phi
1
’
h[u + p + Nhccfi]
Uij
=
’
N
’
’ (2.10)
ifO
(l+jah)C1+(N-i)Bhl
’
h[u + p + Nhc$]
We now multiply
’ ’
’
’
.
(2.8) by A- ‘, obtaining
Y = - $‘BF(Y,
Y’) +
hA-‘Ii
+ A-'R(h)
or
Y = - f+(Y,
Y’) + hD[l]
+ s
(2.11)
where
c = Ppql(N+l) x(N+l)
= A-’ B
(2.12)
with C
py
=
2 aprbrq
(2.13)
*=O
and 0
D = [dp.ql(N+l) x (Ntl,
(2.14) 0 0
0
0
d N-2,N-3
. ..o
0 d N-l.Nm2 d N,N-2
where if p = q = 0 if 4 = P -
1,
if p = N,
q = N - 2
for p = 1,2,
, iv -
1
(2.15)
Two-point
boundary
399
value problems
and
(2.16)
(2.17)
with r n = jto
(n
anjRj(h)
=
l,...,N).
0,
(2.18)
Equation (2.11) expresses Y in terms of F(Y, Y’). In order to derive an iterative procedure to determine Y, it is necessary to couple (2.5) with another set of equations giving Y’ in terms of F( Y, Y’). This will be accomplished in the next section. Before closing this section we note that all cP q’s are nonnegative, and Y, = - ;
yN = - ;
: cOnfn + W4) n=O
5 cNn_fn+ 0(h4) n=O
and N-l
max[IaioI+ laiNI]= ‘fN),
C
which can be verified by direct calculations 3. DETERMINATION
As the number
N is an odd number
s sxN
Iaij(
(2.20)
0fN2)
=
j=l
i
using (2.7), (2.9), (2.13) and (2.17). OF
Y’ IN
by hypothesis,
TERMS
OF
F(Y,
using Simpson’s
Y’)
quadrature
formula
we find
x*i
Y;i - Yb =
f(x, Y, Y') dx
x0
= ;c.f, YX-Y;i-l
+ 4fi + 2f2 + 4& + . . . + 2f2i_2
=
+ 4f2,_i
+
f2iI+ R”;,(h)
(3.1)
f[x,y,y’)dx
x2,-1
+4f2i + 2f2i+l + 4f2i+2+ ...
+ 2fN_2 + 4fN-i
+
f,]
- R;,_,(3.2)
400
M. N. O~IZT~RELI, F. PASCALand K. V. LEUNG
where
RTi_,(h) = 0(h4).
R;,(h) = O(V), Further,
by virtue of the boundary
conditions
(3.3)
(1.2) and (2.1 l), we have the following
relationships
Yb = “y, - 6
=-p
-;
; n=O
[
CN”fn+h[u,063+u,,y)
+y-/I+,.
(3.5)
1
Hence we have the system of equations yb = u - g [
cOnfn+ h(a,,6 +
i
u,,y)] - 6 + ar,
n=O
y;i = i IIf0 + 4fr + 2fz + 4f3 + . . . + 2fzi-2
+ 4fzi_1
+ fzi]
+
R:,(h)
(3.6) y;i-l
=
-~[&,-1
-
+
4fJ_i
+
+4fii+z
2fzi_l
[ 2 $“f”
R;i_l(W + B ;
+ . . . + 2.L2
- &a,,~ + U,,Y) + 1
n=O
Y-
CNnfn- Na,,~ + U,,Y) +
0
0
0 -1
1) x 0
0 -4
(Iv+
-2
1)
0 -4
Y-
0
-2 0
1
0
0
0 -1
1
4
2
4
1
0
0
0
0
0 -1
1
424241
-4
0 -4
0
4
Br,
=
0
1
+ f,]
Br,
1
C4= [4&v+
1 + f,+
+ 4.L
-2
0 -2
0 -4
0
-4 0
-2 0
-4
0
-2. 0
-4 0
-2
0. -2.
0 -4
0
0.
0. -2.
0
0.
.oooooooooo .-4
-2
-4
-2
-4
-2
-4
-2 0
.ooooooo .-4
-2
-4
-2
-4
-2
-4
.oooo
-2
-2 0
.oooooo .-4
-4
-4
-2
-4
-2 0
-2 0
0 -4
0 -4
0
-1
-1 0
-4 0
0
0 -1
0
0,
Two-point
1
4
2
4
2
4
2
0000000 1
4
2
4
2
2
4
0000000 1
4
2
4
2
4
2
0000000
boundary
value problems
0
0
0
4
2
4.
.
2
4
2
4
1
0
0
0.
.
0
0
0
0
0 -1
4
2
4.
.
2
4
2
4
2
4
1
0
0
0.
.
0
0
0
0
0
0
0 -1
4
2
4.
.
2
4
2
4242410
0
0
0.
.
0
0
0
0000000
-4
-2
0 -4
0
0 -1
0 -4
0 -1
If p = even = 2i, we have
AP4 = AZi 4 = ’ for p = 2,4,6, . . . , N -
I
if q = 0 or q = p = 2i
4,
if q = 1,3,5, . . . , p - 1
2,
if q = 2,4,6, . . . , p - 2
0,
if q > p = 2i
1 and q = 0, 1,2,. . . , N. If p = odd = 2i -
;1P.4 = A2i-l,q =
(3.8)
1, we have
-1,
ifq=p=2i-l,orifq=N
-4,
ifq=p+l,p+3
,...,
N-l
-2,
ifq=p+2,p+4
,...,
N-2
[ for p = 2i -
1,
0,
(3.9)
ifq
1 = 1,3,5,. . . , N - 2, and q = 0, 1,2,. . . , N. Let also
[Z”‘]
= [z;;]
=
0
0
0
0
o...o
0
0
0
0
1
1
1
1
l...l
1
1
1
1
0
0
0
0
o...o
0
0
0
0
1
1
1
1
l...l
1
1
1
1
;
(3.10)
0
0
0
0
o...o
0
0
0
0
1
1
1
1
l...l
1
1
1
1
0
0
0
0
o...o
0
0
0
0
1
1
1
1
l...l
1
1
1
1
with 1,
p = 0, 1,2,. . . , N,
if p = odd number = 2i -
Z”’ = 1 0, if p = even number = 2i q=0,1,2 ,..., N
1
IN+
1)
x (N+)
M. N. O~IZT~RELI, F. PASCAL and K. V. LEUNG
402
[P2’] = [I:;]
1
1
1
1
1 ..1
1
1
1
1’
0
0
0
0
0 ..o
0
0
0
0
1
1
.l
1
1 ..1
1
1
1
1
0
0
0
0
0 ..o
0
0
0
0 (3.11)
= 1
1
1
1
0 ..o
0
0
0
0
1
1 ..l
1
1
1
1
0
0 ..o
0
0
0
0
1
1
1
1
1
0
0
0
0
1
1
1
0
0
0
1
(N+
1) x (N+ 1)
where 1’2’ = P.4
1,
if p = even number = 2i
0,
if p = odd number = 2i -
1
and
-CN,oO 0
'N,
0
0
1
0
.
0
0
0
0
.
0
0
0
‘N, 2 ’
0
0
0 (3.12)
[C(N)] = [CT;] = 0
0
0.
‘N,N-2
0
0
0
0
0.
0
CN,N-1
0
0
0
0.
0
0
‘N. ?v_ (.Y+ 1) x (N+
1)
where
0, if p # 4
C(N)= P.4
ifp-q
(‘N, p’
for p = 0, 1,2,. . . , N and q = 0, 1,2,. . . , N, and
[c(O)] =
Co,0
0
0
.
0
0
0
0
Co1
O
.
0
0
0
0
0
0
0
0
Co2.
[cr;] =
(3.13, 0
0
0.
‘O,N-2
0
0
0
0
0.
0
C
0
0
0
0.
0
0
O,N-1
CO.1 (N+l)
x (N+li
Two-point
boundary
403
value problems
where p
P.4
=
1
0, if
P f
c o,p
ifP =
4
4
for p = 0, 1,2,. . . , N and q = 0, 1,2,. . . , N. k,
0
0
0 ...O
0
0
0
0
k,
0
0 ...O
0
0
0
0
0
k,
0 ...O
0
0
0
0
0
0
k,...O
0
0
0
(3.14)
with
0, if p # k P.
4
=
4
k, = -iWa,,~ + a,,~) + k, = ah(a,,d
+ a,,~)
- 6,
Y,
ifp=q=oddnumber=2i-1
if p = q = even number
(3.15) = 2i.
and, finally, let aoro -R?(h) R;(h) T =
kl(NN+ 1)x 1=
- Pr, + ar,
: R:(h)+ two
(3.16)
- Br, with
aor0 tn =
-R’;,._
I
,(h) - /?r,
Rzi(h) + ccrO
I - Dr,
for n = 0 for n = 2i for n = 2i for n = N.
1 (3.17)
M. N. O~LIZT~RELI, F. PASCAL and K. V. LEUNG
404
Then, (3.6) can be written
y’
in the following
h[A]
=
ph”z(1)
+
3 =
matrix
M
form:
W) - $
12
I”’ M”’
[F] + [K] [l]
(3.18)
[Ll t-F1+ [Kl [II + T
where
h
[A]
/3h’1(1) c(N) _
+
12
3 Formula (3.18) determines The iterative procedure section.
~(2) ~‘0’
the iterative
(3.19)
.
Y’ in terms of F = F( Y, Y’), as required. based on the formulas (2.11) and (3.18) will be described
4. ITERATIVE
To describe concisely following matrices:
c12
in the next
PROCEDURE
procedure
associated
with (2.11) and (3.18) we introduce
the
z0
[Z]=
(4.1)
zN
Z N+l
Z
ZN+I
(2N+2)xl
with ifO
L‘n
Z” =
i l‘:,_N_l
[F] =
ifN+ldn<2N+l
iN
(4.2)
(4.3)
'v+1
F 2N+l
(2N+2) X I
with F
=
n
f,,
if 0 d n d N
1,
ifN+l
(4.4)
Two-point
PI
(2Nt2) X (2N+2)
boundary
$1
405
value problems
x (N+l)
(N+l)
=
PI (N+l)
x(N+l)
h[Dl
(N+
1) x (Nf 1)
CKI(N+l)
and
x (N+ 1)
1
(4.5)
(4.6) Then, (2.11) and (3.18) can be combined
together
and rewritten
as
[Z] = [Ql t-W)] + [El Equation
(4.7) can be used as an iterative
formula
to obtain
(4.7) Y,, Y,, . . . , I’,, Yd, Y;, . . . , YA from
Compute L Computea I=l
iii
z=zo
z=newz F = F(Z)
Flowchart
I. Iterative
procedure
to find y,,, y,,
, yn.
M. N. OEULT~RELI, F. PASCAL and K. V. LEUNG
406
an initial guessed value of [Z], say [Z(O)] = [ YdO’,Y;o’, . . . ) YAO),Y;(o), Y;(O),. . . ) Y;‘“‘]‘. Then, the successive
iterations
(4.8)
are defined by the formula
[Z(‘f I)] = [sz][Z”‘] for r = 0, 1,2, . . . . The computer Flowchart I.
implementation
+ [E]
(4.9)
of this iterative
5. APPROXIMATE
procedure
is shown
in the
SOLUTIONS
In this section we deal with the approximate values F and r’ of the solution and its first derivative of the boundary value problem (l.l), (1.2) at the grid points xk. For this purpose we shall investigate the approximation equation [Z] = [Q][F(Z)],
[Z] =
(5.1)
[I ;,
which is obtained from (4.7) omitting the remainder term [El. Let W be a subset of RzN+ ’ with elements W = (W,, WI,. . , WN, W,, l,. . . , W,,, with the norm
1) equipped
(5.2) where K and M are the Lipschitz constants which appear in the inequality (1.5). It can easily be verified that the space W is complete with norm (5.1). Clearly the operator [!A][F( W)] maps W into itself. Let
w=[f,],w*$] be any two elements
(5.3)
of W. Since -y = -(/?/12)CF(I: Y’) + hD[l] -P = LF(Y, Y’) + K[l]
by (2.11), (3.18) and (4.7) with neglected [Q][F(W)
-
F(W”)]
[El, we find =
p;12’cl [F(K
7’) - F(F;*, Y*‘,].
(5.4)
Thus, by virtue of (5.2), we have II[n][F(
W) -
F( W*)] 11d (h2/12)K max (C) om;:N {If,(I:
L’) - f,(r*,
r*‘)l>
. . + M max (L) oyna:N {lf,(X Y’) - f,(p*,
. .
F*‘)I>
(5.5)
Two-point
boundary
407
value problems
where
max(Cl = max n
max (L) = m;x
jfo lCnjl f
(5.6)
ILnjl.
j=O
Further,
the Lipschitz
condition (f”(T V
(1.5) yields r*‘)l < K(j,
- f,(p*,
- y,*j + M\y: - Y,*‘J
from which we find w*ll.
(5.7)
IIPlCmv - eW*)l I/G PIIw - w* II
(5.8)
p = (h2/12) K max (C) + M max (L).
(5.9)
T*‘)( d 11 w -
omn%; Jf”(K P’) - f,(p*, . . Hence, by virtue of (5.5) and (5.7), we have
where
According to the principle of contraction mappings, (5.1) has a unique for any w W* E $f. We shall return to this question after evaluating Determination ofmax (C) andmax (L) From (2.7), (2.9), (2.12) and (2.13), remembering arranging the terms that
j$o
l2
j$o
l2
'n, j =
solution in YT if 0 < p < 1 max (C) and max (L).
that all c~,~‘s are nonnegative,
2an, - 6[an,+ ‘n, Nl’ j
0
we find by re-
(5.10)
j=O
If n = 0, (5.10) becomes 'O,j =
i aO,j j=O
6[a0,0
+
'O,Nl.
(5.11)
Since
1 + (N - j)Ph ‘Ofj= h[a + p + ~$(b - u)]
(OdjdN)
by (2.10) we obtain (N + 1)[2 + Nj3h-j iY ao*j = j=O
h[a
+
p +
ap(b
-
(5.12)
a)3
and 2 + BPh
ao.0 + a0.N = h[a + fl + c$(b-
(5.13)
M. N. OEUZT~RELI, F. PASCAL and K. V. LEUNG
408
Accordingly,
we have 6(b - a)[2 + B(b - a)]
;
co, j = ~___ h2[a
j=O
Further,
+
p +
c$(b
-
(5.14)
a)] .
if 0 < n < N, we have i
a
=
(n + l)[l + (N - nWh][2 + n&h]
n, J
2h[a + p + crfi(b - a)]
j=O
(5.15)
.
and (N - n)(l + nah)[2 + (N - n - l)Bh] 2h[a + p + a&
j=n+l
together
(5.16)
- a)]
with an.0
2 + nah + (N - n)/?h +
'n,N =
~ h[LY +
p +
(5.17)
a)]
a/3(b -
which, for 0 -C n < N, yield the equality 6(b - a)
2 + ah 2n - x ( “3
h2[a + B + c$(b - a)]
j=O
Finally,
[l + b(b - a)] +
for n = N we find a
1 +jah N~j= h[a + p + c&b - a)]
?
(N + 1)(2 + Nah) 'N,j =
j=O
(5.19)
2h[a + B + t$(b - a)] 2 + Nah
aN,O +
'N,N
=
h[a + fi + ap(b - a)]
yielding the equality 6(b - a) IF 'N,j = j=O
h2[a
We now write (5.18) in the following N
6(b - a)
’ C”~j=&$+/?+a/?(b-a)] j=O
+
b +
Q(b
_
a)J[2
+
'db -
form: [2 + P(b - a)] + [2ah + aB(b - a)h]n
-
g+p,+
lb - a)
4W N- 4h
1.
(5.21)
Xc” j has a maximum.
It can be
n2
1 Since the coefficient
(5.20)
~11.
of n2 on the left hand side of (5.21) is negative
Two-point
boundary
409
value problems
easily shown that this maximum occurs for N a[2 + fl(b - a)] =_--p 2 a + p + aB(b - a)
“0
(5.22)
Note that 0 < no < N. Therefore, substituting n by n, on the right hand side of (5.21), and putting 8=
2 + j?(b - a) a + fl + a@ - a)
we find after some simple calculations that max (C) d
c?(b - a) e 6(b hZ@ [ l + 4
(5.23)
1
(5.24)
Similarly, from (3.7) we find max 2 jAP,ql= F jIzI,q( = 5 IAN_l,q( = 3(N - 3) + 6 < 3N. P
(5.25)
q=o
q=o
q=o
Further, we see from (3.10) and (3.12) that max i P
[Z(l)C(N)]p,q= 5 [PC(Q+
q=o
l,q = j$o CN,j = FCa + T+;Jb
_ a)] c2 + a@ -
41.
q=o
(5.26) We also see from (3.11) and (3.13) that max 5 [I(2)C’o’]p,q = 5 [I(%Y]2i,q P
q=O
= i j=O
q=o
6(b - a) [2 + P(b - a)]. ‘0.j = h2[a + p + @(b - a)]
(5.27) Note that 5 [z(r)C(Nqp,q =0
ifp=even
q=o
;
(5.28) =0
[zw(o)]p,q
ifp=odd
q=o
Combining the relationships (5.25H5.28) and (3.19), we find max(L)=(b-a)
l+ [
v[2 + v(b - a)] 2[” + B + aB(b - a)]
1
(5.29)
where v = max{a,fl}.
(5.30)
This concludes the evaluations of max (C) and max (L). Consider now the number p defined by (5.9). Substituting the expressions (5.24) and (5.29) for
410
M.N.
OCXJZTBRELI, F.P~sc~~andK.V.LEUNC
max(C) and max(L) into (5.9) we obtain
P=(b-a$ i
(
Thus, we can state the following 1. If the function
THEOREM
a2(b - a) 4 8
l+
+M 1
theorem
l+V [
f(x, y, y’) satisfies the Lipschitz
conditions
(1.5) in W = I 0 9, and if
2 tl + B + c&b - a) where 8 and v are defined by (5.23) and (5.30), then the system (5.1) has a unique Cw,, w,,.“,
and this solution
wN’ wN+I’.“’
can be constructed
w,,+,)l(~,
WN+i+l)E~‘i
by the iterative
to s.
Moreover,
= O, l,,“‘,N)
(5.32) U in %, (5.33)
procedure r=0,1,2,...
Ur+ l) = [fl][F(U)“‘)], where U(O) is any vector belonging expressed by the inequality
(5.31)
:
2 + v(b - u)
~. = [w=
II’
2 + v(b - a) 2a+p+ap(b-a)
(5.34)
the committed
error in the r-th iteration
is
(5.35) withr
= 0, 1,2,.
..
6.ERRORANALYSI.S
In this section
we prove the following
theorem:
THEOREM 2. Assume that the boundary value problem (l.l), (1.2) has a unique solution of class C6[u, b], the conditions of Theorem 1 are satisfied and that the function f(x, y, y’) is sufficiently
smooth in 9 as required by the Simpson’s rule with remainder term of order 0(h4). Then the unique solution of (5.1) approximates the solution of the problem (l.l), (1.2) as well as its derivative at the grid points {x,},” with an error of order 0(h4). Further, if 2 is the solution of (l.l), (1.2) and U”’ is the r-th iterate associated with the solution U of (5.1), then j/z Proof: Let Z defined
solution
VI/
d &O(V)
+ p’l/U’”
by (4.1) be the exact solution
-
U’O’JI}.
of the problem
(l.l),
(6.1) (1.2) and let U be the
of (5.1). Put E= Z -
Subtracting
U (the error vector).
(6.2)
(5.1) from (4.7) we find s = [Ql[HZ)
- F(U)]
+ [El
(6.3)
and
II4 6 IIPlmz) - ~WIII + IIC~lll~
(6.4)
Two-point
It follows from the inequalities
boundary
411
value problems
(5.8) and (6.4) that (6.5)
Further,
by virtue of (2.17), (3.16), (4.6) and (5.2), we have
II[El II = K OQn
(6.6)
It follows from (2.18) that
max
[Y,\ =
max
(an, + ‘nN) max {lR,@)l, IR,(h)l}
Using(5.17)
jR,(h)(. (6.7)
max
+
> l
04ndN
OQflSN
we find 2 + v(b - a)
(6.8)
max (ano + unN)= h[” + p + c$(b - a)] OQnQN where v is defined by (5.30). Further,
we have
N-l
max
C
G,
anj =
+
4G,G, (6.9)
8h*[a + j? + a/W - a)]
OSn
where Go =
h[a + P + a/?(!~- a)] (6.10)
G, = a@ - a)(2+ l/N) + 3j?h + a&b - a)*(1 + 2/~)
G, = 2W - 1) - 2a(b - a) + j?(b - a)(1 + ~/N)(N - 2) - c$(b - a)*(~ + l/~). It follows from (6.9) and (6.10) that N-l
max
C
unj = 0
Further,
according
$
(6.11)
.
0
OQnQNj=l
to (2.4) and (2.6) we have
R,(h) + O(@),
R,(h) = 0(h5),
R,(h) = 0(h6) (n= 1,2,
,N -
1).
Hence
max Ir,l = O(P)O(l/h) + O(P)O(l/h*)
(6.12)
= 0(h4).
OSnQN
To evaluate
max
It,( consider
(3.17). Accordingly
OQnSN
(6.13) where
R,*'s are the remainder
terms of Simpson’s 1
1 Pxl
quadrature
= W4).
formula
and by assumption (6.14)
412
M. N. OC~UZTBRELI, F. PASCALand K. V. LEUNG
Hence, on account
of(1.12)-(6.14),
we have max
It,1 = 0(h4).
(6.15)
O
Thus, combining
(6.6), (6.12) and (6.14), we find
IIL-ElI/= ah41 and, therefore,
(6.16)
we have (6.17)
ljsli = 0(h4) by virtue of the relationships (6.5) and (6.19). Finally, let us observe that I/z -
VII
which, together
< llz -
UII + IIU -
[II [El II+
V-11 < &
with (6.16), imply the inequality
7. COMPUTER
P’llv
(6.1). This completes
u(l)- q1
the proof of Theorem
2.
PROGRAMMING
The iterative method described in Section 4 has been implemented in the Concordia University Montreal, Quebec, in the PASCAL language for their CDC 6400 Computer System. The most important part of this program is the construction of the matrix OMEGA. This matrix consists of four submatrices which are shown in the diagram in Fig. I, along with their respective sizes. OMEGA is defined once in the program and multiplied by the vector Z several times.
OMEGA (2N + 2) x (2N + 2)
FIG. 1.
In order to save storage locations as well as to minimize execution time we eliminated some trivial and redundant operations by simply taking advantage of the structure of the matrices OMEGAD and OMEGAK. In this way the size of the matrix OMEGA is reduced considerably. Using the notation introduces in the preceeding sections and employing the approximate iteration formula
L-4 = [~l[w)1 we can easily calculate
the i-th element
(7.1)
Zi of [Z] :
~2 N
Zi = -
& ,c J-0
CijFj + Di
(i = 0, 1,. .) N)
(7.2)
413
Two-point boundary value problems
where Di = h i
(i = 0, 1,.
dij
(7.3)
. . , N).
j=O
Similarly we find N
zi = 1
(i = N + 1, N + 2,.
LijFj + Ki
.
,2N + 1)
(7.4)
j=O
where N
Ki =
1 kij
(i = N + 1, N + 2,. . . ,2N + 1).
(7.5)
j=O
Note that the expressions Di and Ki are independent of Fj and can be calculated once and for all. Using (7.2) and (7.4), we construct matrix OMEGA as shown in Fig. 2.
OMEGA OMEGAL
[Ll(N+l)x(N+l)
OMEGAK
L
I
(ZN+Z)x(N+2)
FIG. 2. Further,
let us recall
that
the matrix [L] is defined by (3.19). Since we have
r0 ‘NO
0 z(Uc(N)
=
. ..o
0
0
'Nl 0
'N2 'NN 0 ...o
c 'NN NO CN1 'N2 . ... ... .. ... .. ... .. . 0
0
0
'NO
CN1
CN2... 'NN
(7.6)
0
by virtue of (3.10) and (3.12), and Coo CO1 0
ZWC’O’=
0
Coo CO1 0
0
Co2 0
'ON 0
C02"'CON 0 0
... .. .. ... ... .. ... .. Coo 0
Co1 0
Co2 0
'ON 0
(7.7)
M. N. OC~UZT~RELI,F. PASCAL and K. V. LEUNG
414
by virtue of(3.11) and (3.13), the elements
Lij= for i, j = 0, 1,2, . . . , N. We now introduce new indices OMEGA in Fig. 2. is constructed
Lij of the matrix [L] are given by the equalities if i even
(h/3)Aij - (ah2/12)c,i
(7.8)
if i odd
i(h/3)Aij - (/9?/12)c,,
I = i + 1, J = j + 1. Then, as follows:
in our program,
each block of
OMEGAC:
10 OMEGAD
-(h2/12) NPl = N + 1, NPl 1, NPI (I, J) = - H212* C(1, J)
-H212 = DO10 I = DOlOJ = OMEGAC
: This matrix is a column
vector defined
1
by
NPl
OMEGAD
(I) = h c
(I = 1, N + 1)
d(1, J)
J=l
where a,, 16 + a,,‘v+1y d(I,
J) =
a, ,a + a, ,,,+I”? 1, QN+l,ld
1, NPl (I) = 0,O 1,NPl (I) = OMEGAD (I) = OMEGAD
40
20
DO40 I = DO201 = DO20J = OMEGAL
1, NPl 1,NPl 1,NPl (I, J) = EL(1, J)
where EL(I, J)‘s are computed OMEGAK
vector defined by
NPl
OMEGAK(I)
=
c
K,,
I = 1, NPl
J=l
OMEGAK (1) = ay, - 6 OMEGAK(N + 1) = - By, + y
A- ’ .) In our program
(I) + D(1, J) (I)*h
by (7.8).
: This matrix is a column
1
if1 = N + 1, J = N -
aij of the matrix
DO401 = OMEGAD DO30J = OMEGAD OMEGAD CONTINUE
OMEGAL:
ifJ = I -
+ ’ ‘N+l.N+ly
(Note the shift in the indices of the elements
30
if1 - J = 1
1. we have
Two-point
boundary
415
value problems
and 6 + aN+l,N+l,N+ly)
-fih(%+,,, K(1, J) =
E&,,b 0
The coding
+ a, N+lY) - 6
+ yifI
= J = even
if I = J = odd
if1 #J
of this part is as follows:
DO60 I = 1, NPl OMEGAK (I) = 0 DO 50 J = 1, NPl 50 OMEGAK (I) = OMEGAK (I) + RK(1, J) 60 CONTINUE RK(1, J) = K(1, J), as above. where Computer implementation of Zr + i = QZr, the main program, is as follows: (1) First OMEGA constructed CALL OMEGA (2) y; and y)N+1 are found ZYP(1) = ALFA*ZY(l) - DELTA ZYP(NP1) = -BETA*ZY(NPl) + GAMMA (3) Fi, i = 1,. . .) N + 1 arecomputed DO301 = 1,NPl 30 FZ(I) = FUN(X(I), ZY(I), ZYP(I)) (i.e. fi = f(x, Yi, YI)) (4) multiply OMEGA by F part by part (4.1) multiply OMEGAC by F using OMEGAD as initial value of Zi (i.e., yi): DO 501 = 1,NPl ZNEWY(1) = OMEGAD(1) DO40J = 1,NPl 40 ZNEWY(1) = ZNEWY(1) + OMEGAC(1, J) * FZ(J) 50 CONTINUE (4.2) multiply OMEGAL by F using OMEGAK as initial value of Zi (i.e. yi), i = 2,. . . , N (Y; and YX+ 1 are already known): DO701=2,N ZNEWYP(1) = OMEGAK(1) DO 605 = 1,NPl 60 ZNEWYP(1) = ZNEWP(1) + OMEGA(1, J) * FZ(J) 70 CONTINUE (5) Find lIZnew - Zoldl/2. (6) If /IZ,_ - Zold/l 2 d s STOP otherwise copy Zn_ into Zold and GO TO STEP 2. The above implementation has been diagrammatically shown in Flowchart program can be obtained from the first author.
I. The details of the
8. EXAMPLES
In 121 Stepleman
considered
the following
boundary
value problems
for comparison
purposes:
416
M. N. O~UZT~RELI, F. PASCAL and K. V. LEUNG
(A) y” = t e -“[y’
+ ~‘~1
Y’(O)- Y(O)= 0, y” = -$[e-2Y
y’(1) + y(1) = 2e
+ y’“]
(B)
y’(1) + y(1) = $ + In2 i Y’(O) - Y(O) = 1, (C) Y” = (Y + +qV(l + x) i ‘(0) - $y’(O) = i,
y’(l), + +y(l) = +e.
The solutions of (A) and (C) are y(x) = ex, while (B) has the solution the small difference between the problems (B) here and in [2].) Let E denote the order of the desired accuracy of the computations. Taking N = 51 and the initial guesses
we made a number of computations with different results are summarized in Table 1, below.
E
using our computer
program.
Our
TABLE 1,
E Number of iterations
values of
y(x) = In (1 + x). (Note
10-l
10-Z
10-S
10-h
10-S
10-e
10-7
10-s
10-g
A
6
9
12
14
17
20
23
25
28
B
2
4
6
7
8
9
10
11
11
C
10
16
22
28
34
40
46
53
59
Taking N = 11,21,51,101 we repeated our calculations with different E and different initial guesses. We observed that for a certain E and initial guess for each of the problems (A), (B), (C), the number of iterations are the same. TABLE 2. Number Initial guesses
- 50 -2.0 - 1.0 0.0 0.3 0.5 0.69 1.0 1.5 2.0 2.5 2.7 3.0 4.0 5.0 10.0
B
A CPUEXIT 4
of iterations
2
C
CPUEXIT _
2
2
; 2 6” 5 4
CPUEXIT CPUEXIT 12 12 12
:
i:
2 3
10 7 7 10 10
z t : CPUEXIT CPUEXIT
2 6 2 2
;:
CPUE6XIT
Remark. CPUEXIT 2 means that a very large number the computer stops computation.
1 2
2
occurs
CPU;XIT
during
2
the iteration
process,
so
417
Two-point boundary value problems Keeping N and E fixed we made calculations our results for N = 51 and s = 10-l.
with different initial
APPENDIX Consider
the following
well-known
numerical
We wish to determine the multipliers will be zero in the expression
guesses. Table 2 summarizes
A
differentiation
formulas:
a, b and c in such a way that the coefficients
hC@Y, -
of y,, y, and y,
642)
61+ ah’& + bh2fl+ ch2f2.
Hence a, b and c satisfy the linear system 19a + b - 5c = 6 - 14~ + b + 4c = - 4 lla The unique
solution
- b - c = 3.
of this system is a=$,b=L
4’c
With the above values of a, b, and c, the expression h(ay, - 6) +
$I’&
=
-1
24’
(A2) becomes
+ $h2fl - 2$h2f2 = -y,
y,.
Hence y&l
ah)
y, = - ;
[$f,
+ 3f, -
h6 +
646)
(48)
(A9)
418
M. N. OGUZTBRELI, F. PASCAL and K. V. LEUNG
A&nowledgement---It is the authors’ use of its computer facilities. Obituary-The of the present
first two authors work.
very pleasant
inform
duty to express their sincere thanks
with great sorrow
the sudden
to Concordla
death of Professor
University
for the
K. V. Leung in the final stage
REFERENCES 1 CHAWLA M. M., A fourth-order trigonal finite difference method for general nonlinear two-point boundary-value problems with mixed boundary conditions, J. Inst. Math. Applic. 21, 83-93 (1978). 2 S~EPLEMANR. S., Tridiagonal fourth order approximations to general two points nonlinear boundary value problems with mixed boundary conditions, Math. Comp. 30, 92-103 (1976). 3 RUSSELL R. D., A comparison of collocation and finite differences for two-point boundary value problems, SIAMJ. Numer. Analysis 14, 19-39 (1977). 4 CHANDRA J., LAKSHMIKAN~HAMV. & MITCHELL A. R., Existence of solutions of boundary value problems for nonlinear second order systems in a Banach space, Nonlinear Analysis T.M.A. 2, 157-168 (1978). 5. USMANI R. A., A note on the numerical integration of nonlinear equations with mixed boundary conditions. Utilitas math. 9, 181-192 (1976). Berlin-Heidelberg-New York 6. COLLATZ L., The Numerical Treatment of Differential Equations. Springer-Verlag, (1966). 7. AZIZ A. K. (ed.), Numerical Solutions of Boundary-Value Problems for Ordinary Dtyferential Equations. Academic Press, New York (1975). 8. BAILEY P. B., SHAMPINE L. F. & WALTMAN P. E., Nonlinear Two Point Boundary-Value Problems. Academic Press, New York and London (1968). 9. BELLMANR. & KALABA R., Quasilinearization andNonlinear Boundary- ValueProblems. American Elsevier, New York (1965). 10. BELLMANR., On the iterative solution of two-point boundary value problems, Boll. (in. mat. ital. 16,145-149 (1961). 11. BERNFELD S. & LAKSHMIKAN~HAMV., An Introduction to Nonlinear Boundary Value Problems. Academic Press, New York (1974). 12. DANIEL J. W. & MARTIN A. J., Numerov’s method with deferred corrections for two-point boundary-value problems, SIAM J. Numer. Analysis 14, 1033-1050 (1977). 13. FOX L., The Numerical Solution of Two-Point Boundary Value Problems in Ordinary Differential Equations. Oxford Univ. Press, London and New York (1957). boundary-value problems for nonlinear second order differential systems, SIAM J. 14. HARTMAN P., On two-point math. Analysis 5, 172-177 (1974). 15. HENRICI P., Discrete Variable Methods in Ordinary Differential Equations. Wiley, New York (1962). Two-Point Boundary Value Problems: Shooting Methods. American Elsevier, 16. KELLER H. B., Numerical Methodsfor New York (1972). 17. KELLER H. B. & WHITE A. B. JR., Difference methods for boundary value problems in ordinary differential equations, SIAM J. Numer. Analysis 12, 791-802 (1975). for nonlinear operator equations, Numer. Math. 10, 316-323 (1967). 18. PEREYRA V., Iterated deferred corrections
Nonlinear Analysis. Theory. Merhods & Applicarions, Vol. 3, No. 3, pp. 395-418 Q Pngamon Press Ltd. 1979. Printed in Great Britain.
NUMERICAL
SOLUTIONS OF TWO-POINT VALUE PROBLEMS*
BOUNDARY
M. N. O&JZTBRELI and F. PASCAL Department of Mathematics, University of Alberta, Edmonton, Alberta, Canada, T6G 2Gl and K. V. LEUNG* Department of Computer Science, Concordia University, Montreal, Quebec, Canada, H3G lM8 (Received22
May 197X; revised28 Augusf 197X)
Key words: Nonlinear ordinary differential equations, two-point boundary value problem with mixed conditions, fourth order approximations, iterative method.
1. INTRODUCTION
paper is concerned with the numerical solution of a two-point boundary value problem with mixed boundary conditions for a nonlinear second order ordinary differential equation of the form XEI Y” = f’k Y, Y’), (1.1)
THIS
ay(a) - y’(a) = 6,
(1.2)
PY@) + Y’(b) = Y
where I = [a, b], a, /3, y and 6 given real numbers such that a B 0,
LX+/?>0
B 2 0,
(1.3)
and f(x, y, z) is a real valued function defined in a region 9 c R3, W=IQ9 (l-4) 9 = HY,41 IYJG PI,
I4 G Pz>
in which it satisfies a Lipschitz condition of the form If(&Y,,Y;)
-f(&Yz,Y;)I
G qy,
- Y,l + M/Y;
Y;l
(1.5)
where pt, p2, M and K are constants. The special case in which the function f does not depend explicitly on y’ has been investigated in a number of monographs and papers partially listed in the bibliography. Recently, in [l] and [2], fourth-order tridiagonal finite difference methods for the general problem have been proposed. In [3] the methods of collocation and finite differences have been compared and 0(h4) convergence discussed. In a very recent paper, [4], the existence of solutions for the general nonlinear boundary value problem has been investigated in a Banach space. In the present work we propose an iterative finite difference method which gives globally * This work was partially supported by the National Research Council of Canada by Grant NRC A-4342 to K.V.L. and Grant NRC A-4345 to M.N.O. *&Deceased on I June 1978.
396
M. N. OGUZT~~KELI. F. PASCAL and K. V. LEUNG
0(h4)-convergence over the mesh points of a uniform mesh of width h. We calculated the derivative y’ at the mesh points by a new scheme which differs from those in [l-3 1. We also extended some of the results of [5]. In the next two sections we introduce some preliminary notation and relationships on which our iterative method is based. The iterative procedure is described in Section 4 and is used in Section 5 in connection with the existence and uniqueness of the solutions of the approximating equation which appear in a natural way in our analysis. This section is concluded with our first theorem about approximate solutions and their iterative constructions. In Section 6 we deal with the error analysis in connection with the exact solution of equations (I. l), (1.2) and the approximating iterative solutions considered in the previous section. The results of this section are summarized by our second theorem. In Section 7 a brief explanation of our computer programming is given and the implementation is described in Flowchart I. In the last section we present a short summary of our simulations related with three boundary value problems previously considered in [2] by Stepleman. In Appendix A we establish one of the elementary results used in the second section. 2. PRELIMINARIES
Let N be a positive
odd integer.
Put
h = (b - a)/N,
xn = a + nh
Y’k’= Y’k’(XJ
(2.1)
f, = f(X,, Yn, YA) n = 0, 1, . . . , N).
(k = 0, 1,2;
By virtue of (l.l), (1.2) we have the relationships Yb = UY, -
Y:: = f,?
Using
the following
well known
formula
4
(2.2)
y)N= -py,+s.
of numerical
differentiation
(cf. [6])
+ 2Y” - Yn+i = - g [Y;_ 1 + 10~; + y;+ J + R,,(h)
-y,_i
(2.3)
where R,(h) = O(h6) and the formulae
(A6) and (A8) of Appendix
(1 + ah) Y, - YI = - ;
(n = 1,2,. . , N -
A, (1.1) by the following
[;f,
_Y”_ 1 + 2Yn - Y,+1 = - &_I
-Y,_, with
n = 1,2,. . . , N -
1)
+ 3f,
- $fi]
+ lofn +
system of equations
+ h6 + R,(h)
f,+Il + R,(h)
+ 3fjve1 +
+ (1 + Bh)Y, = - ;[-+r,_,
(2.4)
;&I + b
(2.5)
+ R,(h)
1, where R,(h) = 0(h5),
R,(h)
= 0(h5).
(2.6)
Two-point
boundary
397
value problems
We now put l+ah
- 1
0
2 -1
1
- 1
0
0.. .
0
0
0
0
0.. .
0
0
0
0
.
0
0
0
0
2 -1..
\
A=
I$
0
B=
1
0
2 -1
0
0..
0
0
0..
o-1
0
0
0.. .
0
2 - 1 O-l
1 + bh
(N+
1) x (N+l)
.o
0
0
o\
,(h)
10.
.o
0
0
0
(4
101.
.o
0
0
0
3 -*
110
-1
0
0.
R(h) =
Pp.,1 .l
10
.O
(4
10
1 10 1
.0-t
35
!
(N-1 1) x (N+ 1)
,(h)
Y’=
(N+l)
Then the system (2.5) can be written
x 1
in the matrix
h= AY = - $3F(I: It has been shown in [5] that the inverse A-’
form
Y’) + Hh + of the matrix
R(h).
(2.8)
A is of the form
(2.9) x (N+l)
M. N. O~LT~RELI,
398
F. PASCAL and K. V. LEUNG
with
if 0 < i < j <
(l + ‘ah)[l + (N - j) phi
1
’
h[u + p + Nhccfi]
Uij
=
’
N
’
’ (2.10)
ifO
(l+jah)C1+(N-i)Bhl
’
h[u + p + Nhc$]
We now multiply
’ ’
’
’
.
(2.8) by A- ‘, obtaining
Y = - $‘BF(Y,
Y’) +
hA-‘Ii
+ A-'R(h)
or
Y = - f+(Y,
Y’) + hD[l]
+ s
(2.11)
where
c = Ppql(N+l) x(N+l)
= A-’ B
(2.12)
with C
py
=
2 aprbrq
(2.13)
*=O
and 0
D = [dp.ql(N+l) x (Ntl,
(2.14) 0 0
0
0
d N-2,N-3
. ..o
0 d N-l.Nm2 d N,N-2
where if p = q = 0 if 4 = P -
1,
if p = N,
q = N - 2
for p = 1,2,
, iv -
1
(2.15)
Two-point
boundary
399
value problems
and
(2.16)
(2.17)
with r n = jto
(n
anjRj(h)
=
l,...,N).
0,
(2.18)
Equation (2.11) expresses Y in terms of F(Y, Y’). In order to derive an iterative procedure to determine Y, it is necessary to couple (2.5) with another set of equations giving Y’ in terms of F( Y, Y’). This will be accomplished in the next section. Before closing this section we note that all cP q’s are nonnegative, and Y, = - ;
yN = - ;
: cOnfn + W4) n=O
5 cNn_fn+ 0(h4) n=O
and N-l
max[IaioI+ laiNI]= ‘fN),
C
which can be verified by direct calculations 3. DETERMINATION
As the number
N is an odd number
s sxN
Iaij(
(2.20)
0fN2)
=
j=l
i
using (2.7), (2.9), (2.13) and (2.17). OF
Y’ IN
by hypothesis,
TERMS
OF
F(Y,
using Simpson’s
Y’)
quadrature
formula
we find
x*i
Y;i - Yb =
f(x, Y, Y') dx
x0
= ;c.f, YX-Y;i-l
+ 4fi + 2f2 + 4& + . . . + 2f2i_2
=
+ 4f2,_i
+
f2iI+ R”;,(h)
(3.1)
f[x,y,y’)dx
x2,-1
+4f2i + 2f2i+l + 4f2i+2+ ...
+ 2fN_2 + 4fN-i
+
f,]
- R;,_,(3.2)
400
M. N. O~IZT~RELI, F. PASCALand K. V. LEUNG
where
RTi_,(h) = 0(h4).
R;,(h) = O(V), Further,
by virtue of the boundary
conditions
(3.3)
(1.2) and (2.1 l), we have the following
relationships
Yb = “y, - 6
=-p
-;
; n=O
[
CN”fn+h[u,063+u,,y)
+y-/I+,.
(3.5)
1
Hence we have the system of equations yb = u - g [
cOnfn+ h(a,,6 +
i
u,,y)] - 6 + ar,
n=O
y;i = i IIf0 + 4fr + 2fz + 4f3 + . . . + 2fzi-2
+ 4fzi_1
+ fzi]
+
R:,(h)
(3.6) y;i-l
=
-~[&,-1
-
+
4fJ_i
+
+4fii+z
2fzi_l
[ 2 $“f”
R;i_l(W + B ;
+ . . . + 2.L2
- &a,,~ + U,,Y) + 1
n=O
Y-
CNnfn- Na,,~ + U,,Y) +
0
0
0 -1
1) x 0
0 -4
(Iv+
-2
1)
0 -4
Y-
0
-2 0
1
0
0
0 -1
1
4
2
4
1
0
0
0
0
0 -1
1
424241
-4
0 -4
0
4
Br,
=
0
1
+ f,]
Br,
1
C4= [4&v+
1 + f,+
+ 4.L
-2
0 -2
0 -4
0
-4 0
-2 0
-4
0
-2. 0
-4 0
-2
0. -2.
0 -4
0
0.
0. -2.
0
0.
.oooooooooo .-4
-2
-4
-2
-4
-2
-4
-2 0
.ooooooo .-4
-2
-4
-2
-4
-2
-4
.oooo
-2
-2 0
.oooooo .-4
-4
-4
-2
-4
-2 0
-2 0
0 -4
0 -4
0
-1
-1 0
-4 0
0
0 -1
0
0,
Two-point
1
4
2
4
2
4
2
0000000 1
4
2
4
2
2
4
0000000 1
4
2
4
2
4
2
0000000
boundary
value problems
0
0
0
4
2
4.
.
2
4
2
4
1
0
0
0.
.
0
0
0
0
0 -1
4
2
4.
.
2
4
2
4
2
4
1
0
0
0.
.
0
0
0
0
0
0
0 -1
4
2
4.
.
2
4
2
4242410
0
0
0.
.
0
0
0
0000000
-4
-2
0 -4
0
0 -1
0 -4
0 -1
If p = even = 2i, we have
AP4 = AZi 4 = ’ for p = 2,4,6, . . . , N -
I
if q = 0 or q = p = 2i
4,
if q = 1,3,5, . . . , p - 1
2,
if q = 2,4,6, . . . , p - 2
0,
if q > p = 2i
1 and q = 0, 1,2,. . . , N. If p = odd = 2i -
;1P.4 = A2i-l,q =
(3.8)
1, we have
-1,
ifq=p=2i-l,orifq=N
-4,
ifq=p+l,p+3
,...,
N-l
-2,
ifq=p+2,p+4
,...,
N-2
[ for p = 2i -
1,
0,
(3.9)
ifq
1 = 1,3,5,. . . , N - 2, and q = 0, 1,2,. . . , N. Let also
[Z”‘]
= [z;;]
=
0
0
0
0
o...o
0
0
0
0
1
1
1
1
l...l
1
1
1
1
0
0
0
0
o...o
0
0
0
0
1
1
1
1
l...l
1
1
1
1
;
(3.10)
0
0
0
0
o...o
0
0
0
0
1
1
1
1
l...l
1
1
1
1
0
0
0
0
o...o
0
0
0
0
1
1
1
1
l...l
1
1
1
1
with 1,
p = 0, 1,2,. . . , N,
if p = odd number = 2i -
Z”’ = 1 0, if p = even number = 2i q=0,1,2 ,..., N
1
IN+
1)
x (N+)
M. N. O~IZT~RELI, F. PASCAL and K. V. LEUNG
402
[P2’] = [I:;]
1
1
1
1
1 ..1
1
1
1
1’
0
0
0
0
0 ..o
0
0
0
0
1
1
.l
1
1 ..1
1
1
1
1
0
0
0
0
0 ..o
0
0
0
0 (3.11)
= 1
1
1
1
0 ..o
0
0
0
0
1
1 ..l
1
1
1
1
0
0 ..o
0
0
0
0
1
1
1
1
1
0
0
0
0
1
1
1
0
0
0
1
(N+
1) x (N+ 1)
where 1’2’ = P.4
1,
if p = even number = 2i
0,
if p = odd number = 2i -
1
and
-CN,oO 0
'N,
0
0
1
0
.
0
0
0
0
.
0
0
0
‘N, 2 ’
0
0
0 (3.12)
[C(N)] = [CT;] = 0
0
0.
‘N,N-2
0
0
0
0
0.
0
CN,N-1
0
0
0
0.
0
0
‘N. ?v_ (.Y+ 1) x (N+
1)
where
0, if p # 4
C(N)= P.4
ifp-q
(‘N, p’
for p = 0, 1,2,. . . , N and q = 0, 1,2,. . . , N, and
[c(O)] =
Co,0
0
0
.
0
0
0
0
Co1
O
.
0
0
0
0
0
0
0
0
Co2.
[cr;] =
(3.13, 0
0
0.
‘O,N-2
0
0
0
0
0.
0
C
0
0
0
0.
0
0
O,N-1
CO.1 (N+l)
x (N+li
Two-point
boundary
403
value problems
where p
P.4
=
1
0, if
P f
c o,p
ifP =
4
4
for p = 0, 1,2,. . . , N and q = 0, 1,2,. . . , N. k,
0
0
0 ...O
0
0
0
0
k,
0
0 ...O
0
0
0
0
0
k,
0 ...O
0
0
0
0
0
0
k,...O
0
0
0
(3.14)
with
0, if p # k P.
4
=
4
k, = -iWa,,~ + a,,~) + k, = ah(a,,d
+ a,,~)
- 6,
Y,
ifp=q=oddnumber=2i-1
if p = q = even number
(3.15) = 2i.
and, finally, let aoro -R?(h) R;(h) T =
kl(NN+ 1)x 1=
- Pr, + ar,
: R:(h)+ two
(3.16)
- Br, with
aor0 tn =
-R’;,._
I
,(h) - /?r,
Rzi(h) + ccrO
I - Dr,
for n = 0 for n = 2i for n = 2i for n = N.
1 (3.17)
M. N. O~LIZT~RELI, F. PASCAL and K. V. LEUNG
404
Then, (3.6) can be written
y’
in the following
h[A]
=
ph”z(1)
+
3 =
matrix
M
form:
W) - $
12
I”’ M”’
[F] + [K] [l]
(3.18)
[Ll t-F1+ [Kl [II + T
where
h
[A]
/3h’1(1) c(N) _
+
12
3 Formula (3.18) determines The iterative procedure section.
~(2) ~‘0’
the iterative
(3.19)
.
Y’ in terms of F = F( Y, Y’), as required. based on the formulas (2.11) and (3.18) will be described
4. ITERATIVE
To describe concisely following matrices:
c12
in the next
PROCEDURE
procedure
associated
with (2.11) and (3.18) we introduce
the
z0
[Z]=
(4.1)
zN
Z N+l
Z
ZN+I
(2N+2)xl
with ifO
L‘n
Z” =
i l‘:,_N_l
[F] =
ifN+ldn<2N+l
iN
(4.2)
(4.3)
'v+1
F 2N+l
(2N+2) X I
with F
=
n
f,,
if 0 d n d N
1,
ifN+l
(4.4)
Two-point
PI
(2Nt2) X (2N+2)
boundary
$1
405
value problems
x (N+l)
(N+l)
=
PI (N+l)
x(N+l)
h[Dl
(N+
1) x (Nf 1)
CKI(N+l)
and
x (N+ 1)
1
(4.5)
(4.6) Then, (2.11) and (3.18) can be combined
together
and rewritten
as
[Z] = [Ql t-W)] + [El Equation
(4.7) can be used as an iterative
formula
to obtain
(4.7) Y,, Y,, . . . , I’,, Yd, Y;, . . . , YA from
Compute L Computea I=l
iii
z=zo
z=newz F = F(Z)
Flowchart
I. Iterative
procedure
to find y,,, y,,
, yn.
M. N. OEULT~RELI, F. PASCAL and K. V. LEUNG
406
an initial guessed value of [Z], say [Z(O)] = [ YdO’,Y;o’, . . . ) YAO),Y;(o), Y;(O),. . . ) Y;‘“‘]‘. Then, the successive
iterations
(4.8)
are defined by the formula
[Z(‘f I)] = [sz][Z”‘] for r = 0, 1,2, . . . . The computer Flowchart I.
implementation
+ [E]
(4.9)
of this iterative
5. APPROXIMATE
procedure
is shown
in the
SOLUTIONS
In this section we deal with the approximate values F and r’ of the solution and its first derivative of the boundary value problem (l.l), (1.2) at the grid points xk. For this purpose we shall investigate the approximation equation [Z] = [Q][F(Z)],
[Z] =
(5.1)
[I ;,
which is obtained from (4.7) omitting the remainder term [El. Let W be a subset of RzN+ ’ with elements W = (W,, WI,. . , WN, W,, l,. . . , W,,, with the norm
1) equipped
(5.2) where K and M are the Lipschitz constants which appear in the inequality (1.5). It can easily be verified that the space W is complete with norm (5.1). Clearly the operator [!A][F( W)] maps W into itself. Let
w=[f,],w*$] be any two elements
(5.3)
of W. Since -y = -(/?/12)CF(I: Y’) + hD[l] -P = LF(Y, Y’) + K[l]
by (2.11), (3.18) and (4.7) with neglected [Q][F(W)
-
F(W”)]
[El, we find =
p;12’cl [F(K
7’) - F(F;*, Y*‘,].
(5.4)
Thus, by virtue of (5.2), we have II[n][F(
W) -
F( W*)] 11d (h2/12)K max (C) om;:N {If,(I:
L’) - f,(r*,
r*‘)l>
. . + M max (L) oyna:N {lf,(X Y’) - f,(p*,
. .
F*‘)I>
(5.5)
Two-point
boundary
407
value problems
where
max(Cl = max n
max (L) = m;x
jfo lCnjl f
(5.6)
ILnjl.
j=O
Further,
the Lipschitz
condition (f”(T V
(1.5) yields r*‘)l < K(j,
- f,(p*,
- y,*j + M\y: - Y,*‘J
from which we find w*ll.
(5.7)
IIPlCmv - eW*)l I/G PIIw - w* II
(5.8)
p = (h2/12) K max (C) + M max (L).
(5.9)
T*‘)( d 11 w -
omn%; Jf”(K P’) - f,(p*, . . Hence, by virtue of (5.5) and (5.7), we have
where
According to the principle of contraction mappings, (5.1) has a unique for any w W* E $f. We shall return to this question after evaluating Determination ofmax (C) andmax (L) From (2.7), (2.9), (2.12) and (2.13), remembering arranging the terms that
j$o
l2
j$o
l2
'n, j =
solution in YT if 0 < p < 1 max (C) and max (L).
that all c~,~‘s are nonnegative,
2an, - 6[an,+ ‘n, Nl’ j
0
we find by re-
(5.10)
j=O
If n = 0, (5.10) becomes 'O,j =
i aO,j j=O
6[a0,0
+
'O,Nl.
(5.11)
Since
1 + (N - j)Ph ‘Ofj= h[a + p + ~$(b - u)]
(OdjdN)
by (2.10) we obtain (N + 1)[2 + Nj3h-j iY ao*j = j=O
h[a
+
p +
ap(b
-
(5.12)
a)3
and 2 + BPh
ao.0 + a0.N = h[a + fl + c$(b-
(5.13)
M. N. OEUZT~RELI, F. PASCAL and K. V. LEUNG
408
Accordingly,
we have 6(b - a)[2 + B(b - a)]
;
co, j = ~___ h2[a
j=O
Further,
+
p +
c$(b
-
(5.14)
a)] .
if 0 < n < N, we have i
a
=
(n + l)[l + (N - nWh][2 + n&h]
n, J
2h[a + p + crfi(b - a)]
j=O
(5.15)
.
and (N - n)(l + nah)[2 + (N - n - l)Bh] 2h[a + p + a&
j=n+l
together
(5.16)
- a)]
with an.0
2 + nah + (N - n)/?h +
'n,N =
~ h[LY +
p +
(5.17)
a)]
a/3(b -
which, for 0 -C n < N, yield the equality 6(b - a)
2 + ah 2n - x ( “3
h2[a + B + c$(b - a)]
j=O
Finally,
[l + b(b - a)] +
for n = N we find a
1 +jah N~j= h[a + p + c&b - a)]
?
(N + 1)(2 + Nah) 'N,j =
j=O
(5.19)
2h[a + B + t$(b - a)] 2 + Nah
aN,O +
'N,N
=
h[a + fi + ap(b - a)]
yielding the equality 6(b - a) IF 'N,j = j=O
h2[a
We now write (5.18) in the following N
6(b - a)
’ C”~j=&$+/?+a/?(b-a)] j=O
+
b +
Q(b
_
a)J[2
+
'db -
form: [2 + P(b - a)] + [2ah + aB(b - a)h]n
-
g+p,+
lb - a)
4W N- 4h
1.
(5.21)
Xc” j has a maximum.
It can be
n2
1 Since the coefficient
(5.20)
~11.
of n2 on the left hand side of (5.21) is negative
Two-point
boundary
409
value problems
easily shown that this maximum occurs for N a[2 + fl(b - a)] =_--p 2 a + p + aB(b - a)
“0
(5.22)
Note that 0 < no < N. Therefore, substituting n by n, on the right hand side of (5.21), and putting 8=
2 + j?(b - a) a + fl + a@ - a)
we find after some simple calculations that max (C) d
c?(b - a) e 6(b hZ@ [ l + 4
(5.23)
1
(5.24)
Similarly, from (3.7) we find max 2 jAP,ql= F jIzI,q( = 5 IAN_l,q( = 3(N - 3) + 6 < 3N. P
(5.25)
q=o
q=o
q=o
Further, we see from (3.10) and (3.12) that max i P
[Z(l)C(N)]p,q= 5 [PC(Q+
q=o
l,q = j$o CN,j = FCa + T+;Jb
_ a)] c2 + a@ -
41.
q=o
(5.26) We also see from (3.11) and (3.13) that max 5 [I(2)C’o’]p,q = 5 [I(%Y]2i,q P
q=O
= i j=O
q=o
6(b - a) [2 + P(b - a)]. ‘0.j = h2[a + p + @(b - a)]
(5.27) Note that 5 [z(r)C(Nqp,q =0
ifp=even
q=o
;
(5.28) =0
[zw(o)]p,q
ifp=odd
q=o
Combining the relationships (5.25H5.28) and (3.19), we find max(L)=(b-a)
l+ [
v[2 + v(b - a)] 2[” + B + aB(b - a)]
1
(5.29)
where v = max{a,fl}.
(5.30)
This concludes the evaluations of max (C) and max (L). Consider now the number p defined by (5.9). Substituting the expressions (5.24) and (5.29) for
410
M.N.
OCXJZTBRELI, F.P~sc~~andK.V.LEUNC
max(C) and max(L) into (5.9) we obtain
P=(b-a$ i
(
Thus, we can state the following 1. If the function
THEOREM
a2(b - a) 4 8
l+
+M 1
theorem
l+V [
f(x, y, y’) satisfies the Lipschitz
conditions
(1.5) in W = I 0 9, and if
2 tl + B + c&b - a) where 8 and v are defined by (5.23) and (5.30), then the system (5.1) has a unique Cw,, w,,.“,
and this solution
wN’ wN+I’.“’
can be constructed
w,,+,)l(~,
WN+i+l)E~‘i
by the iterative
to s.
Moreover,
= O, l,,“‘,N)
(5.32) U in %, (5.33)
procedure r=0,1,2,...
Ur+ l) = [fl][F(U)“‘)], where U(O) is any vector belonging expressed by the inequality
(5.31)
:
2 + v(b - u)
~. = [w=
II’
2 + v(b - a) 2a+p+ap(b-a)
(5.34)
the committed
error in the r-th iteration
is
(5.35) withr
= 0, 1,2,.
..
6.ERRORANALYSI.S
In this section
we prove the following
theorem:
THEOREM 2. Assume that the boundary value problem (l.l), (1.2) has a unique solution of class C6[u, b], the conditions of Theorem 1 are satisfied and that the function f(x, y, y’) is sufficiently
smooth in 9 as required by the Simpson’s rule with remainder term of order 0(h4). Then the unique solution of (5.1) approximates the solution of the problem (l.l), (1.2) as well as its derivative at the grid points {x,},” with an error of order 0(h4). Further, if 2 is the solution of (l.l), (1.2) and U”’ is the r-th iterate associated with the solution U of (5.1), then j/z Proof: Let Z defined
solution
VI/
d &O(V)
+ p’l/U’”
by (4.1) be the exact solution
-
U’O’JI}.
of the problem
(l.l),
(6.1) (1.2) and let U be the
of (5.1). Put E= Z -
Subtracting
U (the error vector).
(6.2)
(5.1) from (4.7) we find s = [Ql[HZ)
- F(U)]
+ [El
(6.3)
and
II4 6 IIPlmz) - ~WIII + IIC~lll~
(6.4)
Two-point
It follows from the inequalities
boundary
411
value problems
(5.8) and (6.4) that (6.5)
Further,
by virtue of (2.17), (3.16), (4.6) and (5.2), we have
II[El II = K OQn
(6.6)
It follows from (2.18) that
max
[Y,\ =
max
(an, + ‘nN) max {lR,@)l, IR,(h)l}
Using(5.17)
jR,(h)(. (6.7)
max
+
> l
04ndN
OQflSN
we find 2 + v(b - a)
(6.8)
max (ano + unN)= h[” + p + c$(b - a)] OQnQN where v is defined by (5.30). Further,
we have
N-l
max
C
G,
anj =
+
4G,G, (6.9)
8h*[a + j? + a/W - a)]
OSn
where Go =
h[a + P + a/?(!~- a)] (6.10)
G, = a@ - a)(2+ l/N) + 3j?h + a&b - a)*(1 + 2/~)
G, = 2W - 1) - 2a(b - a) + j?(b - a)(1 + ~/N)(N - 2) - c$(b - a)*(~ + l/~). It follows from (6.9) and (6.10) that N-l
max
C
unj = 0
Further,
according
$
(6.11)
.
0
OQnQNj=l
to (2.4) and (2.6) we have
R,(h) + O(@),
R,(h) = 0(h5),
R,(h) = 0(h6) (n= 1,2,
,N -
1).
Hence
max Ir,l = O(P)O(l/h) + O(P)O(l/h*)
(6.12)
= 0(h4).
OSnQN
To evaluate
max
It,( consider
(3.17). Accordingly
OQnSN
(6.13) where
R,*'s are the remainder
terms of Simpson’s 1
1 Pxl
quadrature
= W4).
formula
and by assumption (6.14)
412
M. N. OC~UZTBRELI, F. PASCALand K. V. LEUNG
Hence, on account
of(1.12)-(6.14),
we have max
It,1 = 0(h4).
(6.15)
O
Thus, combining
(6.6), (6.12) and (6.14), we find
IIL-ElI/= ah41 and, therefore,
(6.16)
we have (6.17)
ljsli = 0(h4) by virtue of the relationships (6.5) and (6.19). Finally, let us observe that I/z -
VII
which, together
< llz -
UII + IIU -
[II [El II+
V-11 < &
with (6.16), imply the inequality
7. COMPUTER
P’llv
(6.1). This completes
u(l)- q1
the proof of Theorem
2.
PROGRAMMING
The iterative method described in Section 4 has been implemented in the Concordia University Montreal, Quebec, in the PASCAL language for their CDC 6400 Computer System. The most important part of this program is the construction of the matrix OMEGA. This matrix consists of four submatrices which are shown in the diagram in Fig. I, along with their respective sizes. OMEGA is defined once in the program and multiplied by the vector Z several times.
OMEGA (2N + 2) x (2N + 2)
FIG. 1.
In order to save storage locations as well as to minimize execution time we eliminated some trivial and redundant operations by simply taking advantage of the structure of the matrices OMEGAD and OMEGAK. In this way the size of the matrix OMEGA is reduced considerably. Using the notation introduces in the preceeding sections and employing the approximate iteration formula
L-4 = [~l[w)1 we can easily calculate
the i-th element
(7.1)
Zi of [Z] :
~2 N
Zi = -
& ,c J-0
CijFj + Di
(i = 0, 1,. .) N)
(7.2)
413
Two-point boundary value problems
where Di = h i
(i = 0, 1,.
dij
(7.3)
. . , N).
j=O
Similarly we find N
zi = 1
(i = N + 1, N + 2,.
LijFj + Ki
.
,2N + 1)
(7.4)
j=O
where N
Ki =
1 kij
(i = N + 1, N + 2,. . . ,2N + 1).
(7.5)
j=O
Note that the expressions Di and Ki are independent of Fj and can be calculated once and for all. Using (7.2) and (7.4), we construct matrix OMEGA as shown in Fig. 2.
OMEGA OMEGAL
[Ll(N+l)x(N+l)
OMEGAK
L
I
(ZN+Z)x(N+2)
FIG. 2. Further,
let us recall
that
the matrix [L] is defined by (3.19). Since we have
r0 ‘NO
0 z(Uc(N)
=
. ..o
0
0
'Nl 0
'N2 'NN 0 ...o
c 'NN NO CN1 'N2 . ... ... .. ... .. ... .. . 0
0
0
'NO
CN1
CN2... 'NN
(7.6)
0
by virtue of (3.10) and (3.12), and Coo CO1 0
ZWC’O’=
0
Coo CO1 0
0
Co2 0
'ON 0
C02"'CON 0 0
... .. .. ... ... .. ... .. Coo 0
Co1 0
Co2 0
'ON 0
(7.7)
M. N. OC~UZT~RELI,F. PASCAL and K. V. LEUNG
414
by virtue of(3.11) and (3.13), the elements
Lij= for i, j = 0, 1,2, . . . , N. We now introduce new indices OMEGA in Fig. 2. is constructed
Lij of the matrix [L] are given by the equalities if i even
(h/3)Aij - (ah2/12)c,i
(7.8)
if i odd
i(h/3)Aij - (/9?/12)c,,
I = i + 1, J = j + 1. Then, as follows:
in our program,
each block of
OMEGAC:
10 OMEGAD
-(h2/12) NPl = N + 1, NPl 1, NPI (I, J) = - H212* C(1, J)
-H212 = DO10 I = DOlOJ = OMEGAC
: This matrix is a column
vector defined
1
by
NPl
OMEGAD
(I) = h c
(I = 1, N + 1)
d(1, J)
J=l
where a,, 16 + a,,‘v+1y d(I,
J) =
a, ,a + a, ,,,+I”? 1, QN+l,ld
1, NPl (I) = 0,O 1,NPl (I) = OMEGAD (I) = OMEGAD
40
20
DO40 I = DO201 = DO20J = OMEGAL
1, NPl 1,NPl 1,NPl (I, J) = EL(1, J)
where EL(I, J)‘s are computed OMEGAK
vector defined by
NPl
OMEGAK(I)
=
c
K,,
I = 1, NPl
J=l
OMEGAK (1) = ay, - 6 OMEGAK(N + 1) = - By, + y
A- ’ .) In our program
(I) + D(1, J) (I)*h
by (7.8).
: This matrix is a column
1
if1 = N + 1, J = N -
aij of the matrix
DO401 = OMEGAD DO30J = OMEGAD OMEGAD CONTINUE
OMEGAL:
ifJ = I -
+ ’ ‘N+l.N+ly
(Note the shift in the indices of the elements
30
if1 - J = 1
1. we have
Two-point
boundary
415
value problems
and 6 + aN+l,N+l,N+ly)
-fih(%+,,, K(1, J) =
E&,,b 0
The coding
+ a, N+lY) - 6
+ yifI
= J = even
if I = J = odd
if1 #J
of this part is as follows:
DO60 I = 1, NPl OMEGAK (I) = 0 DO 50 J = 1, NPl 50 OMEGAK (I) = OMEGAK (I) + RK(1, J) 60 CONTINUE RK(1, J) = K(1, J), as above. where Computer implementation of Zr + i = QZr, the main program, is as follows: (1) First OMEGA constructed CALL OMEGA (2) y; and y)N+1 are found ZYP(1) = ALFA*ZY(l) - DELTA ZYP(NP1) = -BETA*ZY(NPl) + GAMMA (3) Fi, i = 1,. . .) N + 1 arecomputed DO301 = 1,NPl 30 FZ(I) = FUN(X(I), ZY(I), ZYP(I)) (i.e. fi = f(x, Yi, YI)) (4) multiply OMEGA by F part by part (4.1) multiply OMEGAC by F using OMEGAD as initial value of Zi (i.e., yi): DO 501 = 1,NPl ZNEWY(1) = OMEGAD(1) DO40J = 1,NPl 40 ZNEWY(1) = ZNEWY(1) + OMEGAC(1, J) * FZ(J) 50 CONTINUE (4.2) multiply OMEGAL by F using OMEGAK as initial value of Zi (i.e. yi), i = 2,. . . , N (Y; and YX+ 1 are already known): DO701=2,N ZNEWYP(1) = OMEGAK(1) DO 605 = 1,NPl 60 ZNEWYP(1) = ZNEWP(1) + OMEGA(1, J) * FZ(J) 70 CONTINUE (5) Find lIZnew - Zoldl/2. (6) If /IZ,_ - Zold/l 2 d s STOP otherwise copy Zn_ into Zold and GO TO STEP 2. The above implementation has been diagrammatically shown in Flowchart program can be obtained from the first author.
I. The details of the
8. EXAMPLES
In 121 Stepleman
considered
the following
boundary
value problems
for comparison
purposes:
416
M. N. O~UZT~RELI, F. PASCAL and K. V. LEUNG
(A) y” = t e -“[y’
+ ~‘~1
Y’(O)- Y(O)= 0, y” = -$[e-2Y
y’(1) + y(1) = 2e
+ y’“]
(B)
y’(1) + y(1) = $ + In2 i Y’(O) - Y(O) = 1, (C) Y” = (Y + +qV(l + x) i ‘(0) - $y’(O) = i,
y’(l), + +y(l) = +e.
The solutions of (A) and (C) are y(x) = ex, while (B) has the solution the small difference between the problems (B) here and in [2].) Let E denote the order of the desired accuracy of the computations. Taking N = 51 and the initial guesses
we made a number of computations with different results are summarized in Table 1, below.
E
using our computer
program.
Our
TABLE 1,
E Number of iterations
values of
y(x) = In (1 + x). (Note
10-l
10-Z
10-S
10-h
10-S
10-e
10-7
10-s
10-g
A
6
9
12
14
17
20
23
25
28
B
2
4
6
7
8
9
10
11
11
C
10
16
22
28
34
40
46
53
59
Taking N = 11,21,51,101 we repeated our calculations with different E and different initial guesses. We observed that for a certain E and initial guess for each of the problems (A), (B), (C), the number of iterations are the same. TABLE 2. Number Initial guesses
- 50 -2.0 - 1.0 0.0 0.3 0.5 0.69 1.0 1.5 2.0 2.5 2.7 3.0 4.0 5.0 10.0
B
A CPUEXIT 4
of iterations
2
C
CPUEXIT _
2
2
; 2 6” 5 4
CPUEXIT CPUEXIT 12 12 12
:
i:
2 3
10 7 7 10 10
z t : CPUEXIT CPUEXIT
2 6 2 2
;:
CPUE6XIT
Remark. CPUEXIT 2 means that a very large number the computer stops computation.
1 2
2
occurs
CPU;XIT
during
2
the iteration
process,
so
417
Two-point boundary value problems Keeping N and E fixed we made calculations our results for N = 51 and s = 10-l.
with different initial
APPENDIX Consider
the following
well-known
numerical
We wish to determine the multipliers will be zero in the expression
guesses. Table 2 summarizes
A
differentiation
formulas:
a, b and c in such a way that the coefficients
hC@Y, -
of y,, y, and y,
642)
61+ ah’& + bh2fl+ ch2f2.
Hence a, b and c satisfy the linear system 19a + b - 5c = 6 - 14~ + b + 4c = - 4 lla The unique
solution
- b - c = 3.
of this system is a=$,b=L
4’c
With the above values of a, b, and c, the expression h(ay, - 6) +
$I’&
=
-1
24’
(A2) becomes
+ $h2fl - 2$h2f2 = -y,
y,.
Hence y&l
ah)
y, = - ;
[$f,
+ 3f, -
h6 +
646)
(48)
(A9)
418
M. N. OGUZTBRELI, F. PASCAL and K. V. LEUNG
A&nowledgement---It is the authors’ use of its computer facilities. Obituary-The of the present
first two authors work.
very pleasant
inform
duty to express their sincere thanks
with great sorrow
the sudden
to Concordla
death of Professor
University
for the
K. V. Leung in the final stage
REFERENCES 1 CHAWLA M. M., A fourth-order trigonal finite difference method for general nonlinear two-point boundary-value problems with mixed boundary conditions, J. Inst. Math. Applic. 21, 83-93 (1978). 2 S~EPLEMANR. S., Tridiagonal fourth order approximations to general two points nonlinear boundary value problems with mixed boundary conditions, Math. Comp. 30, 92-103 (1976). 3 RUSSELL R. D., A comparison of collocation and finite differences for two-point boundary value problems, SIAMJ. Numer. Analysis 14, 19-39 (1977). 4 CHANDRA J., LAKSHMIKAN~HAMV. & MITCHELL A. R., Existence of solutions of boundary value problems for nonlinear second order systems in a Banach space, Nonlinear Analysis T.M.A. 2, 157-168 (1978). 5. USMANI R. A., A note on the numerical integration of nonlinear equations with mixed boundary conditions. Utilitas math. 9, 181-192 (1976). Berlin-Heidelberg-New York 6. COLLATZ L., The Numerical Treatment of Differential Equations. Springer-Verlag, (1966). 7. AZIZ A. K. (ed.), Numerical Solutions of Boundary-Value Problems for Ordinary Dtyferential Equations. Academic Press, New York (1975). 8. BAILEY P. B., SHAMPINE L. F. & WALTMAN P. E., Nonlinear Two Point Boundary-Value Problems. Academic Press, New York and London (1968). 9. BELLMANR. & KALABA R., Quasilinearization andNonlinear Boundary- ValueProblems. American Elsevier, New York (1965). 10. BELLMANR., On the iterative solution of two-point boundary value problems, Boll. (in. mat. ital. 16,145-149 (1961). 11. BERNFELD S. & LAKSHMIKAN~HAMV., An Introduction to Nonlinear Boundary Value Problems. Academic Press, New York (1974). 12. DANIEL J. W. & MARTIN A. J., Numerov’s method with deferred corrections for two-point boundary-value problems, SIAM J. Numer. Analysis 14, 1033-1050 (1977). 13. FOX L., The Numerical Solution of Two-Point Boundary Value Problems in Ordinary Differential Equations. Oxford Univ. Press, London and New York (1957). boundary-value problems for nonlinear second order differential systems, SIAM J. 14. HARTMAN P., On two-point math. Analysis 5, 172-177 (1974). 15. HENRICI P., Discrete Variable Methods in Ordinary Differential Equations. Wiley, New York (1962). Two-Point Boundary Value Problems: Shooting Methods. American Elsevier, 16. KELLER H. B., Numerical Methodsfor New York (1972). 17. KELLER H. B. & WHITE A. B. JR., Difference methods for boundary value problems in ordinary differential equations, SIAM J. Numer. Analysis 12, 791-802 (1975). for nonlinear operator equations, Numer. Math. 10, 316-323 (1967). 18. PEREYRA V., Iterated deferred corrections