Chemical
Engineering
Science,
1974, Vol. 29, pp. 1695-1699.
Pergamon
Press.
Printed
in Great Britain
SOLUTION OF NONLINEAR BOUNDARY VALUE PROBLEMS-VIII. EVALUATION OF BRANCHING POINTS BASED ON SHOOTING METHOD AND GPM TECHNIQUE MILAN Department
of Chemical
KUBB?EK
Engineering,
and
VLADIMiR
HLAVACEK
Institute of Chemical Technology, Prague 6, Czechoslovakia
(Receiued
14 November
Suchbatorova
1903, 166 28
1973)
Abstract-General straightforward method for finding branching points of nonlinear boundary value problems is presented. The technique proposed utilizes the equations used in the GPM method. The method can be applied to problems arising in a number of physical applications, the technique is tailored to diffusion and heat conduction problems.
1. INTRODUCTION
subject to linear boundary
The nonlinear boundary value problems often may exhibit for a given set of parameters more than one solution. The kinds of problems we treat arise in wide variety of applications. In particular, some problems in chemical engineering as, for instance, axial mixing in tubular reactors, heat and mass transfer within and outside of a porous catalyst, boundary layer problems for exothermic reaction, equilibrium of neighbouring drops at different potentials etc. provide excellent motivation for our study. Of course, it is a difficult task to establish all possible solutions if more than one is expected. If the solutions are very close to each other the majority of straightforward methods fail to calculate the relevant profiles because of oscillation of a particular iteration procedure between both solutions. If both neighbouring profiles coincide a branching profile results. Because this profile corresponds in a functional space to a point the term branching point will be used in the following text. A great savings of effort can usually be made if it is possible to evaluate all branching points before calculation of solutions for particular values of governing parameters. It is the purpose of this paper to present a straightforward method for an evaluation of the branching points appearing in a two-point nonlinear boundary value problem. The method of solution will be outlined and employed to solve two practical problems. 2. DEVELOPMENT
@Y(a)
by’(b)
= 7%
(24
= ~1
(2b)
a function of the parameter initial condition
LY.On chasing the
as
y(a)
and supposing
= A
(3)
POZ 0, Eq. (2a) yields y’(a)=+-an&.
(4)
In accordance with our recent paper[l] a set of differential equations for the auxiliary variables
p(x)=w$
aY(x)
q(x)=aa
(5)
can be devloped[2]:
pfdfp +-$ p’
p(a)=
p’(a)=
1,
-2
af af q~~=-4+-rq~+-g. q(a)=q’(a)=O. ay a~ After rewritting
(7)
Eq. (2b) we have
F,(A, a) = a,y(b, A, a) + Pgy’(b, A, a) -
yl =
0 (8)
A necessary condition for a branching point yields: aF,
da zi -=---_& dA aF,
aa (1)
0
(6)
OF THE METHOD
Y, Y’, a)
+ by’(a)
a,y(b)+
2.1 One second order equation Consider the boundary value problem Y” = fk
conditions:
This equation 1695
gives, supposing
aF,laa # 0,
1696
MILAN
KUBftEK
and VLADIMIR HLAVXEK
F,(A,u)=$$=nip(b,A,n)+~,p’(b,A,a)=O Equations (8) and (9) are sufficient the unknown coordinates (A, &) of branching point. For instance, using method to solve (8) and (9) the Jacobian
(9) to establish a particular the Newton matrix r
r=
tions (3), (4) (6) (7) and (15) from x = a to x = b. (3) enumeration of F, and Ft according to (8) and (9) (4) enumeration of the elements of the Jacobian matrix (10) according to (11) and (16). (5) Evaluation of a new approximation of A and (Y using the Newton method
(10)
must be determined. The values of the first derivatives evaluated because
of F, can be
~=a,p(b,A,a)+Plp’(b,A,a) (11) $
where Ae(0, 1 > is a damping factor. (6) The calculation is stopped for ll(A, a)“‘” (A, CX)~‘~[~ < E (e is the prescribed tolerance), if not go to step 2. The method developed will be illustrated on a calculated example. Example 1. Heat and mass transfer in a porous catalyst of the plate shape is described by[5]: Y) $ =Sy exp [ 1+rP(JP(1- -Y) 1
= culq(b, A, a) + P,q’(b, A, a)
y$-0,
To calculate the first derivatives of F2 two additional auxiliary equations for the variables
(12) must be developed. These equations result after differentiation with respect to A and (Y as well:
y(l)=
1.
For y(0) = A and setting QI= 6 we can readily develop a set of differential Eqs. (6), (7), (13) and (14); for instance for r we have:
of (6) I”=
YP
6 [
(
YPY (l+P(l-y)P
(l+P(l-y))’
)
_rP(1+P(1-Y))+W2Y (1 +pu-
st&fpq +-&+P’I’+P%)+&& a’f a2
Apparently,
+-$st.
the initial conditions
are
r(u) = r’(u) = s(a) = s’(u) = 0. Now the first derivatives evaluated:
$ = alr(b,
(14)
(15)
of F2 can be easily
A, a)+ P,r’(b, A, a) (16)
g
= cx,s(b, A, a)+ P,s’(b, A, a).
A general algorithm for the calculation of Luand A is following: (1) initial guess of A and (Y (2) integration of 5 second order differential Eq. (l), (6), (7), (13) and (14) with the initial condi-
I
rP(l- Y) exp [ l+p(l-y)
I n* rP(l-Y) I exp [ l+p(l-y)
YPY
0 +-+J’+&p~+$s a(Y )
YP
+6 [ ‘-(l+p(l-y))Z
I r.
The relevant initial value problems have been integrated by virtue of Merson’s method[6]. Some results of our computations are shown in Table 1. The convergence properties of this procedure are very good. 2.2 Two second order equations We shall develop in this section an algorithm a set of two second order differential equations, general algorithm will be presented elsewhere However, the development for more equations be made in the same way. Consider a set of two differential equations: y:‘= fl(X,
Yl,
Y2,Yl, Yi, a)
y: = ji(X,
Yl,
Y2,Y’l, YL a)
for the [3]. can
(18)
Solution of nonlinear boundary value problems-VIII
1697
Table 1. Course of iteration for single second order equation (y = 20, p = 0.4, A = 1) Iteration
A
a=6
0 1 2 3
0.8000 0.7954 0.7929 0.7928
0
1 2 3 0 1 2 3 4
Ft
F,
0.1000 0.1367 0.1376 0.1376
- 0.04772 - 0~00104 0~00002
0.23304 0.01286 - owOO3
0.2000 0.2218 0.2272 0.2273
0.1000 0.0771 0.0779 0.0779
0.15636 - 0wKxO - OTlOO16
- 0.02447 ow470 0~00103
0.6000 0.8832 0.8024 0.7928 0.7928
0.1000 0.1451 0.1463 0.1372 0.1376
- 0.03936 0.01903 0.01065 - omo4o
- 0.31023 0.24825 - 0.01402 0.00165
subject to linear boundary
conditions
After inserting the values of y and y’ (which are the results of integration of (18), (20) and (21) into (19b) two nonlinear equations can be written R(A,, AZ, a) =
cukydb,AI, Azr a)
+ ~,~y;(b,A,,AZ, (Y)- ylk = 0
(24)
k = 1,2. The necessary condition for the branching point is K(A,,
ao’yi(a) + Poiy’r(a) = yoi
i = 1,2
(19a)
cul’yi(b)+ &‘y’,(b) = yli
i = 1,2.
Wb)
A~,(Y) = det WI,
here T is the Jacobian
AZ, a) ‘0,
(25)
matrix
Let us select the following two initial conditions y,(a)=A,, and supposing
yz(a)=Az
(20)
(26)
l-=
PO’Z 0, i = 1,2, we arrive at
y\(a)=L(yo’-c~o’A~) PO’
i=l,2.
(21)
In the same way as in the one-dimensional case a set of differential equations for auxiliary variables:
am)
Piik =
A!$
qi(x)
Pii T-j-&
a*Ydx) a~,a~k
qij@) =
(22)
zt.l&
This condition results from the implicit function theorem. Clearly, the branching point can be evaluated as a result of solving three nonlinear Eqs. (24) and (25) which yield the unknowns A,, A2 and (Y.If the Newton method is used to solve these nonlinear equations the terms aE/aA,, aTlao: i = 1,2,3; j = 1,2 must be evaluated. To evaluate all these terms the following relations can be developed:
I
i,j,k=l,2;
$=
kaj
can be developed [ 1,4] with the initial conditions = Pij(a)
Pii
=O,
i# j;
pii
= 1,
I
$=u,iq~(b)+Pl’q:(b) aF3 -=_
aAi
p:,(+-$ C&(U)= q:(a) = = qij(a)
pij!f(a)
= q:j(ll)
i, j, k = 1,2,
For instance,
for qij we have
=
P:jk(a)
= 0,
aFI
(23)
(27a)
i = 1,2 aF,
a2F2
a2F,
aA, m+%m-aA, ___aF2 a2F, aA, aA2aAj
cl
i, j = 1,2
cUl’pij(h)+ P,‘p:j(b)
Wb) aF,
a’F2
aA,aAj j=1,2
here, e.g. a2F2 -= aA2aAj
k 3 j.
and
a12p,,j(b) + fl1*pLj(b)r
(%a)
MILAN KUB~~EK and VLADIMIRHLAVACEK
1698
a& -=-
aa
aF, a2F2 aF, a2F2 aF, d2F, aA,ZGG+faA2aA,acvZaA,a~
---aF2 a2F, aA,
(28b)
akhaa'
here, e.g. we have
2
=a,‘q,,(b)+ P,*qL(b).
2
To illustrate the method for a set of two second order differential equations an example will be presented. Example 2. Axial heat and mass transfer in tubular reactors in described by two second-order differential equations [7,8] y - PeMy’ + PeMDaU - Y) em 0” - PeH8’ -
e
(- > 1 + le
=o
Pe&(O - 8,) +PeHDaB(l-y)exp
subject to boundary
(
&
>
=0
conditions
y’(1) = V(l) = 0 PeMy(0) - y’(0) = 0
e(l)=At,
0.1
a05
We denote y(l)=A,,
1
0.0
PeHe(0) - e’(o) = 0.
a=m
I 0.7!
Fig. 1. Dependences A,(a) and A*(a) for Pe, = PeM = 2, /3 =2.& =E =O,B = 12.
a=Da.
Table 2. Course of iteration for the values of parameters PeH=Pem=2,h=1. Iteration
A, = y(1)
A2 = e(l)
a = Da
0 1 2 3
0.2000 0.3759 0.3871 0.3871
1mOO 1.8061 1.8652 1.8652
0~1000 0.1330 0.1319 o-1319
0.00469 0~00003 0~0ooOo
0 1 2 3
0~7.500 0.7933 0.7860 0.7858
3-8000 4.0761 4.0313 4.02%
0.1150 0.1108 0.1115 0.1115
-0.01787 0.00424 OTIOO18
0 1 2 3 4 5 6
09000 0.9713 0.9533 0.9309 0.9207 0.9185 0.9184
4.0000 3.6593 36495 3.8366 3.8911 3.8990 3.8993
0.1400 0.1338 0.1206 0.1230 0.1216 0.1214 0.1214
0 1 2 3 4
0.9800 0.9677 0.970 1 0.9696 0.%96
3.8000 3.3720 3.3476 3.3527 3.3529
0.1150 0.1175 0.1154 0.1155 0.1155
F,
B = 12, p = 2, E = 8, = 0,
F*
F3
0.36868 OTIO138 0@0001
- 099197 -0aO444 -0aKKIo
- 0.02230 OdO815 OdOO56
0.45113 -0.11603 - oaO357
- 0.03548 0.23797 0.06945 0.02078 oaO393 0.00016
0.50492 244774 0.61371 0.22125 0.03630 ool145
- 1.89113 8.58409 3.94400 0.80386 0.11828 0.00421
0.47920 0.00328 - oaO113 - omOO3
3.97732 0.09019 - 0.01242 - 0.00038
20.8162 1.31692 - 040002 - 0~01564
Solution of nonlinear boundary
To find the branching points it is necessary to integrate 18 differential equations of second order (initial value problem). The course of the iteration process is reported in Table 2. The particular branching points are displayed in Fig. 1. where the dependences A,(a) and At(a) are drawn. It can be noted that the rate of convergence of the Newton iteration process is high. 3. DISCUSSION
AND CONCLUSIONS
There are numerous situations in which one is interested only in an evaluation of the branching points, i.e. the determination of boundaries where multiple solutions occur, not in the details of the whole dependences between some parameters and dependent variables. When such is the case there is no need to evaluate these dependences because the straightforward method proposed locates easily the parameters of the branching points in question. So far the location of a branching point was established only in an indirect way, i.e. it was necessary to explore the vicinity of a particular branching point by mapping or other methods (e.g. GPM technique [ 1,4,9]), however, the method suggested has proved effective because of its straight-forward character. Since obviously nonlinear elliptic partial differential equations can be converted by virtue of the ‘method of lines’ to a nonlinear boundary value problem for a set of ordinary differential equations the technique of finding the branching points can be also easily adopted. It should be
value problems-VIII
1699
noted, however, that there is a number of chemical engineering problems which cannot be integrated by shooting methods because of inherent instability of a relevant initial value problem. Clearly, the algorithm described can be also adapted to nonlinear algebraic equations; the boundary value problem can be then converted by making use of orthogonal collocation or high order difference scheme to a set of nonlinear algebraic equations which may be handled easily. The algorithm for algebraic equations will be presented elsewhere[3]. Further applications will include, probably, studies of thermal explosion equations, boundary layer structure, BCnard convection problems etc. REFERENCES
[l] Kubfcek M. and HlavaEek V., Chem. Engng Sci. 1972 27 743. [2] Coddington E. A. and Levinson N., Theory of Ordinarv Differential Eauations. McGraw-Hill, New York 1955. [3] Kubfcek M., .I. Inst. Math. Appfic., to be published. 141 Kubic’ek M. and Hlavacek V., J. Inst. Math. Applic. 1973. [5] HlavSlEek V., Marek M. and Kubfcek M., Chem. Engng Sci. I%8 23 1083. [6] Fox L., Numerical Solution of Ordinary and Partial Diflerential Equations. Pergamon Press, Oxford, (1%2). [7] HlavaEek V. and Hofmann H., Chem. Engng Sci. 1970 25 1517. [S] Hlavafek V., Hofmann H. and Kubicek M., Chem. 1973 12 287. [9] Kubicek M., Communs of the ACM, 1974.
Zusammenfassung-Eine allgkmeine direkte Methode zur Bestimmung der Verzweigungspunkte einer Die vorgeschlagene Technik verwendet die nichlinearen Randwertaufgabe wird angeftlhrt. Gleichungen, die fur die GPM Methode entwickelt wurden. Das Verfahren kann zur Liisung von vielen physikalischen Problemen verwendet werden, der Algorithmus ist zur Analyse der Diffusionsbezw. Warmeleitf5higkeitsprobleme geeignet.