Evaluation of branching points for nonlinear boundary-value problems based on the GPM technique

Evaluation of Branching Points for Nonlinear Boundary-Value Problems Based on the GPM Technique M. KubiEek Department

of Chemical

Engineering

Institute of Chemical Technology 16628 Prague 6, Czechoslovakia

Transmitted by R. Kalaba

ABSTRACT A novel general direct iteration algorithm for the evaluation of branching points of nonlinear boundary-value problems is suggested. The method proposed takes advantage of the CPM concept. The analysis presented in the paper encompasses also the determination of branching points in nonlinear algebraic equations. The method can be applied to problems arising in a number of physical and mathematical applications. The technique is tailored to problems of diffusion and heat conduction accompanied by a chemical reaction.

1.

INTRODUCTION

Frequently the nonlinear boundary-value problems in ordinary differential equations may possess (for prescribed values of parameters) multiple solutions. To be able to perform a parametric study, we are often faced with the problem of locating all branches of the solution. The points where the solution branches are refered to as the branching points. To establish numerically the solutions near a branching point is a difficult task, since the solutions evaluating

can be very close all branching

each

points,

other. the

Evidently,

computer

time

if we are capable expenditure

of

may be

significantly lowered. Clearly, for an approximate determination of the branching points, the GPM concept [4,5] combined with an interpolation procedure can be used. Since the branching points are important from the physical and computational point of view and, moreover, the literature dealing with this subject is very sparse, a detailed examination of this question is merited. APPLIED MATHEMATICS

AND COMPl TATZON I,341352

0 American Elsevier PublishingCompany, Inc., 1975

(1975)

341

342

M.KUBicEK

The gist of this paper is to propose a new iteration algorithm which takes advantage of both the Newton method and the GPM technique. In the first part of this communication an algorithm for evaluation of branching points in nonlinear algebraic equations will be developed, while in the second part a generalization to nonlinear differential equations which makes use of the GPM method will be presented. The effectiveness of the procedure is demonstrated on three calculated examples. 2.

BRANCHING POINTS FOR A SYSTEM OF NONLINEAR ALGEBRAIC EQUATION%NEWTON METHOD

Consider a system of nonlinear algebraic actual parameter (Y: Fi(x

equations depending on an

i=1,2 ,..., N.

~,...,~~I~)=o,

(I)

Let us study the dependence of the solution x(e) of Eq. (1) on the parameter (Y,with regard to determination of branching points. The necessary condition to establish the branching points follows from the implicit-fimction theorem (e.g. [71): x,,a)=detG(r, F~+~(~~>...,

,..., xN,a)=O.

(2)

Here G is the Jacobian matrix with the elements

aFi(xl,...,x,,a) gi=

axi

i,i= 1,2 ,..., N.

’

It is evident that Eqs. (1) and (2) are capable of determining the coordinates (x?, * ** , x&a*) of the branching point in question. Eqs. (1) and (2) constitute a system of N + 1 nonlinear algebraic equations for N + 1 unknowns; to solve them we can make use of the Newton method: Xk+l=Xk-XII~l(Xk)F(Xk).

(4)

Here we have denoted X = (xl, . . . , r,, CX)~,F = (F,, . . . , FN+ JT, A E elements of the matrix l?r, yii,are aFi

yii =

j=1,2 ,..., N,

ax, ’

aFi Yi,N+l-

(0,1). The

--

aa

’

i=1,2 ,...,N+l.

(5)

Branching Points in Bow&

y-Value

343

Problems

To avoid the divergence of the iteration process, the condition

is to be tested. Each iteration is started with A= 1, and if the above presented inequality is not fulfilled, the value of X is, e.g., halved. The expression for yij, i+ iV+ 1, can be easily developed. From the definition of the determinant we obtain 8detG 7=(-l)i+‘Gii, I

(6)

where Gir is a minor, i.e., a determinant of the N- 1 by N - 1 matrix which arises from G by eliminating the row i and the column i. For the partial derivatives of Eq. (2) we have

YN+

3.

1,m

=!!!&=

2 m

(-‘)i+‘c,i(&-

i,i-1

axi ax, ’

DEVELOPMENT OF THE METHOD FOR A NONLINEAR BOUNDARY-VALUE

m=1,2

,...,

N,

(7)

PROBLEM

Consider a set of ordinary differential equations

dYi

,,=f;(t,yl,...,yn,a),

i=1,2 ,..., n,

(9)

subject to linear two-point boundary conditions 5

aiiyj(u)=ci,

i=1,2 ,..., 7,

(104

2

biiyi(b)=di,

i=1,2 ,..., 12-r.

(lob)

j-1

j-1

M. KUBiCEK

344

A development of the technique for problems described by one or two second-order differential equations instead of (9) is presented in our recent paper [6]. On choosing y*(u)=?Jr,

i=l,2

,...) n-r,

(11)

and after solving the equations (lOa) for yd(a), i = n - r + 1,. . . , n, we can obtain the relations n--T

yi(a)=hd+

x $Vj, j=l

i=n-r+l,...,n.

(12)

Here I+ and htr are the calculated coefficients, described in detail in our recent paper [4]. Differentiation of Eq. (9) with respect to n1 gives rise to auxiliary differential equations

p;i=

c ay n

%

m-l

P,f,

i=l ,...,n,

i=l,...,n-r,

Pa)

m

subject to the iinitial conditions PfiW =

1,

P&4

=o,

i,i=

i=n-r+l,...,n,

Pfj ( 4 = hp

1,2 ,...,n-r,

(14a)

i=1,2 ,..., n-r.

(14b)

Here we have denoted

p&)=%$L

(15)

I

The system of equations which must be satisfied at a branching point is

Fi (~1,.. ., ?n-r,LY)=

5 b~jyj(b>V9a)-di=o>

i=l,2

j=l

,...,n-r

(Isa)

and

F

n-,+1

(1)1,***,77n_r,(Y)=detG(171,....l)n-r,(Y)=01

(16h)

Branching Points in Bout&y-Value

345

Problems

where the elements gm = aFi/ CIqmof the Jacobian matrix G are

gm=

2 bifpe(b,q,a),

,...,

i,m=l,2

n-r,

(17)

j=l

To adopt the Newton method it is necessary to evaluate the second derivatives of the functions Fi. To develop the equations for the second derivatives we differentiate the set of differential equations (9) and (13) with respect to (Y and Q:

(184

i=1,2 ,...,n,

i,k=l,2

,..., n-r,

j
piik=pikj

Pb)

1

&"p,*+-$qjk+'$

Pjk

i=1,2 ,..., n,

7

k=l,2

,,.., n-r.

(18~)

The initial conditions are %(a) =o,

Pijkca)

=‘*

(19)

%k(a)=O*

We have denoted

p+(t)=

$$ 1 k

2Yi qiklt)=

2

*

(20)

M. KUBicEK

343

Clearly, the second and first derivatives are

i,m,k=1,2

,..., n-r,

(214

w-4 and

$= i:bijpi,(b,T,a), m

2

(224

j=l

= i:

biiqj(b,77,a).

Pb)

j=l

We have all variables necessary for making use the Newton method (4) where N=n-r, xi=qi, i=l, 2 ,..., n - r. To recapitulate the algorithm, the following steps are to be performed: (1) Make initial guess of $, . . . ,T&,, a’; set k =O. (2) Integrate Eqs. (9), (13), (18) with initial conditions (ll), (12), (14), (19) from t = a to t= b. The overall number of differential equations to be integrated is n(n-r+l) [2+(n-r)/2]. (3) Calculate the values Fi, i = 1, 2,. . . ,n - r + 1, from Eq. (16). (4) Determine the partial derivatives defined by Eqs. (21), (22), (7) and (8).

(5) Perform one step of the Newton method with an appropriately chosen

A. a k I> E, set increase k by 1 and go to step 2. (6) If l/~k+l-qkkl+Iak+l(E is the prescribed tolerance.) The method can be easily extended to problems with nonlinear boundary conditions containing at t = b the parameter a-i.e., instead of (lo), the conditions qi( yl(b),...,y,(b),a)=O,

i=1,2 ,.*., n-r,

(lob’)

Branching Points in Boundary-Value

Problems

347

are considered. Now Eqs. (Isa), (17), (21a), (2Ib), (22a) and (22b) have to be rewritten considering Eq. (lob’), e.g.,

-=

3%

+ayiPjmkP) ’

i,m,k=1,2

(2la’)

,..., n-r.

1 For mixed boundary conditions cp,(!@)V..3

tln(a),yl(b),...,y”(b)Ya)=O,

i=l,2

,..., n,

(IO’)

it seems convenient to choose al1 initial conditions yi(“)=17i9

i=1,2 ,a**, n;

(II’)

and the system of n nonlinear equations

(16a’) can be handled in a similar way. The method is both straightforward and simple; however, the dimensionality may be a limiting factor for large systems of equations. 4.

APPLICATION

The method developed will be illustrated with three calculated examples arising in chemical engineering.

348

M.KUBiCEK

EAMPLE 1. Heat and mass transfer is described by [3]

dY

in a porous catalyst with a plate shape

(0)

7

y(l)=l.

‘0,

(24)

For y(0) = 17 we can readily develop a set of differential equations (13) and (18) after rewriting Eq. (23) as a system of two first-order equations. The relevant initial-value problems have been integrated by means of Merson’s method (e.g., [I]). Some results are shown in Table 1.

TABLE 1 COURSE

OF ITERATION

Iteration

FOR

EXAMPLE

1 (y = %),p

=

0.6,h = 1)

17

a

F,

F2

1 2 3 4 5 0 1 2 3 4 5

0.8000 0.9480 0.8738 0.8778 0.8743 0.8742 0.0500 0.0967 0.1394 0.1618 0.1654 0.1654

0.1600 0.1752 0.0437 0.0815 0.0860 0.0859 0.0100 0.0156 0.0165 0.0169 0.0170 0.0170

0.14116 0.07120 - 0.05618 - 0.00535 0.00013

- 1.10595 - 0.07091 0.40819 0.05825 -0.00011

- 0.45920

4.02920

-0.09500

1.38898 0.38211 0.04914 0.00097

There it is evident

that

0

F&~~)=p(l,v,~)>

- 0.02516 - 0.00428 -060011

where

aY p=-. a77

We can obtain two different branching points, e.g., CY:= 0.0170, (II; = 0.0859 for y =20, p=O.S. The boundary-value problem (23), (24) possesses a unique solution in the intervals (YE (0, a:) and CY E (cuz, co); on the other hand, in the interval OL~(cu:,ad) three different solutions exist [3].

Branching Points in Boundary-Value Pmblems

349

EXAMPLE 2. Axial heat and mass transfer in tubular reactors is described by the second-order differential equations [2] 1 d2y 4 ---dt

Pe dt

+a(l- y)exp@=O

(254

-J-d28-~-/3(B-B,)+aB(1-y)exp8=0

Pe dt

subject to boundary conditions

dY(l) -=dt

W) =. dt

’

Pey(O)-

dY (0) 7 =0,

de(0)=O.

Pee (0) - 7

(26)

Here Pe, /3, 0,, B are governing parameters, and y and 0 are the conversion and the dimensionless temperature respectively. We denote the missing initial conditions by e(l)=772.

Y(l)=%

(27)

To find the branching points it is necessary to integrate 36 first-order differential equations (initial-value problem) from t = 1 to t =O in each iteration, The course of the iteration process is reported in Table 2. The particular branching points are displayed in Fig. 1, where the dependences qr( LX)and q2( CX)are plotted. EXAMPLE 3. The explosion of solid explosives is described by (e.g. [3]) $++$+aexp(&)=O,

(28) d8(0) -g-- =o,

Nue(l)+

a(l)

7

=o.

(29)

Here y, Nu and p ( = 0,1,2) are physical parameters and 8 is a dimensionless temperature. We choose the missing initial condition at t =0:

e (0) =

71.

(30)

The resulting (Y* and v* for a sequence of values of the parameter Nu are shown in Table 3. The boundary-value problem (28), (29) has two solutions in the interval (YE (0, LX*) and no solution in the remaining interval (YE (a*, @J).

M. KUBicEK

350 TABLE 2 COURSE

OF ITERATION

(B=12,P=2,Bc=0, Iteration

FOR EXAMPLE

2

Pe=2,X=l).

?l

92

a

0

0.0000

0.0000

0.0000

1

2 3 4 5

0.1890 0.3394 0.3818 0.3871 0.3871

0.8519 1.6410 1.8353 1.8647 1.8652

0.1890 0.1299 0.1324 0.1319 0.1319

0 1 2 3 4

0.8000 0.7759 0.7820 0.7858 0.7858

5.0000 4.0535 4.0232 4.0299 4.0296

0.1000 0.1215 0.1115 0.1115 0.1115

0 1 2 3 4 5 6 7

0.9900 0.9793 0.9763 0.9728 0.9704 0.9697 0.9696 0.9696

4.0000 3.0865 3.2581 3.3170 3.3447 3.3524 3.3529 3.3529

0.1000 0.1133 0.1136 0.1149 0.1154 0.1155 0.1155 0.1155

0 1 2 3 4

0.9000 0.9281 0.9199 0.9185 0.9184

4.0000 3.8692 3.8950 3.8992 3.8993

0.1200 0.1226 0.1215 0.1214 0.1214

TABLE 3 RESULTING

BRANCHING

POINTS

FOR EXAMPLE

3 (-f = 20, p = 2).

NU

a*

v*

100 50 20 10 5 2 1 0.5

3.457 3.389 3.196 2.906 2.430 1.564 0.951 0.526

1.820 1.820 1.815 1.801 1.753 1.577 1.400 1.271

Bmnching Points in Bow&y-Value

Pmblems

FIG. 1. The dependences tag

5.

351

and ~)*(a) for Example 2.

CONCLUSIONS

In physical and engineering practice there are a number of situations in which we address ourselves to the evaluation of the branching points, i.e., the determination of boundaries where multiple solutions occur. So far the location of branching points has been established only in an indirect way, i.e., it has been necessary to investigate the neighborhood of the branching point in question by mapping or other similar methods (e.g., the CPM technique). However, the technique proposed has proved economical because of its straightforward character. Of course, this method can also be employed advantageously towards the evaluation of branching points of nonlinear algebraic systems. After an appropriate discretization procedure the method can also be adapted to establish the branching points of elliptic partial differential equations.

352

M. KUBiGEK

REFERENCES 1

L. Fox, Numerical Solution of Ordinary and Partial Differential Equations, Pergamon, New York, 1962. 2 V. HlavaEek, H. Hofmann and M. KubiEek, Modeling of chemical reactorsXXIV. Transient axial heat and mass transfer in tubular reactors. The stability considerations-II, Chem. Eng. Sci. 26 (1971), 16291634. 3 V. HlavaEek, M. Marek and M. KubiEek, Modeling of chemical reactors-X. Multiple solutions of enthalpy and mass balances for a catalytic reaction within a porous catalyst particle, Chem. Eng. Sci. 23 (1968), 1083-1097. 4. M. KubiEek and V. Hlavaeek, General parameter mapping technique-a procedure for solution of non-linear boundary value problems depending on an actual parameter, 1. ht. Math. A&. 12 (1973), 287-293. 5 M. KubiEek and V. HlavPEek, Solution of nonlinear boundary value problems-Va and Vb. A novel method: general parameter mapping (GPM) and predictorcorrector GPM method, Chem. Eng. Sci. 27 (1972), 743-750, 2095-2098. 6 M. KubiEek and V. Hlavacek, Solution of nonlinear boundary value problemsVIII. Evaluation of branching points based on shooting method and GPM technique, Chem. Eng. Sci. 29 (19743, 16951699. 7 T. L. Saaty and J. Bram, Non-linear Mathematics, McGraw-Hill, New York, 1964.