- Email: [email protected]
Solution of nonlinear boundary value problems—VII A novel method: differentiation with respect to boundary condition
ChemicalEngineering
Science. 1973, Vol. 28. pp. 1049-1052.
Pergmon
Press.
Printed inGreat
Britain
Solution of nonlinear boundary value problems-VII A novel method: differentiation with respect to boundary condition MILAN KUBICEK and VLADIMfR HLAVACEK Department of Chemical Engineering, Institute of Chemical Technology, Suchbatarova 1903, Prague 6, Czechoslovakia (Received 24 August 1972) Abstract-A simple numerical method for construction of the dependence of solutions to nonlinear boundary value problem on a parameter will be developed. The set of differential equations is diiferentiated with respect to the boundary condition chosen and the resulting partial differential equations are solved by a finite-difference method. The method is illustrated by an example of heat and mass transfer in a porous catalyst. 1. INTRODUCTION
Now supposing p0 # 0 we will choose
THERE ARE
few methods available for calculating the dependence of a solution of nonlinear boundary value problem on the given parameters so far[l-41. The basic technique used, general parameter mapping [2,3] and differentiation with respect to the actual parameter[ 1,4] transform the original boundary value problem into some initial value problem; in the first case the necessity of meeting boundary conditions requires solution of an initial value problem with an implicit right-hand side while the other method leads to a third order partial differential equation. The necessity of unique solution is the primary disadvantage of the last technique. For multiple solutions serious difficulties arise. It is the purpose of this article to develop a modified technique which may be readily used for handling problems with multiple solutions. 2. DEVELOPMENT
OF THE METHOD
We will assume, for simplicity, a single differential equation of second order Y” -“f-k
y, Y’, aY) = 0
where f is a nonlinear conditions are abY
function.
+PoY’(a)
%Y (b) +PIY
(3)
=A,
here A will be considered as a new parameter. Furthermore, we will assume that the function f is differentiable with respect to all arguments; on differentiation we get a partial differential equation for y(x, A)
aS_mY+2c?%+~dar_ axZaA_
ayaA ay’axaA
aadA
(4)
For a given A the boundary conditions are
cw(a,A) +Po
aY(GA) ax
=
Yo
(5)
cr,y(b,A)+P,~=y,
(6)
y(u,A)
(7)
= A.
Thus we have three boundary conditions. It seems that the problem is over-determined, however, in the partial differential equation (4) we have one unknown parameter
(1)
do C=J-J.
The boundary
= 3/o
(6) =
Y(U)
(2)
We note that the problem is completely mined. Now on choosing A=Ao:a=c+,,
Yl.
1049
y(x,A,)
=y,,(x)
deter(9)
MILAN
KUBiCEK
and VLADIMiR
HLAVACEK
we can integrate the differential equation (4); here y&) is the solution of Eq. (1) with boundary conditions (2) and an additional condition y,(a) = A. Frequently, the solution y,(x) may be readily calculated, sometimes for a = twoa trivial solution exists. 3. NUMERICAL
SOLUTION
ORDER
j+l
OF THE THIRD
EQUATION
For the numerical solution of Eq. (4) a finite difference method has been used. The approximation of the particular derivatives in the mesh point (i,j) is presented in both Table 1 and Fig. 1. The nonlinear terms @/lay, aflay’, and aflaa have to be computed at the old profile j. It is obvious that an evaluation of the new profile demands the solution of a set of linear equations having special almost tridiagonal matrix. For the new profile we must calculate the values of both y(+’ (i = 0, 1, . . ., N) and C. After dividing the region (a, b) into N parts following difference equations result
i-l
j
:
1
Y
Fig. 1. A sketch of the grid of the mesh points.
Note that the profile j has been calculated for Ai. This procedure yields a set of N+ 4 linear equations (lo)-(13) for N +4 variables yJ_:‘, 5+1 . ., yhyl, C. The set of equations has a 2 YIJ+’ Yo special foh so that following calculation procedure can be used:
1. Guess C = Cl. For i = 0 the values yl_:‘, j+l, Yo ylf+l can be calculated from Eqs. (lo), (11) &y&J + s;y2+1+ sj+,y;;; + rjc = 4: and (13). The recurrence relation (10) yields (10) . ., yi;;. In general, the values yi?i, Yi+l, Ys’+l i=O,l,...,N y,+!+l and >iyl do not meet Eq. (12) and a residue where the coefficients s:_, , s(, di+l, $, 4 me de- R (Cl ) depending on the estimation of C, C = Cl, pendent on the known values j$, = y,,’ (n = i- 1, results. i, i+ 1) and 6 = 01’. For the approximation of 2. Guess C = C, # Cl and repeat step 1. Beboundary conditions (5)-(7) we have cause the residue R(C) is a linear function of C the correct value of C can be calculated by linear Y1J+’ -ti’ = y 0 (11) interpolation. abYd+l+ PO 2h 3. For the correct value of C a new profile j+ 1 must be recalculated or interpolated from E-yKi_ Yl (12) two profiles evaluated for C = Cl and C = C,.. fflYN1+' +Pl y 4. This new profile yif+l may be employed to yoj+1 = Al+ k. (13) correct the coefficients sC1, st, s:+~, riJ, qrJ; these will be calculated for the average value y, deTable 1. Approximation of derivatives in fined as arithmetical mean, i.e. Eq. (4) 2h
-
ay _ Y,‘+‘-Y*’ aA k _apY _ YE: -Ye:
axaA _ a? _ ax2aA
yi = 4(Y2+Y(+‘), -Yd+1 +Ycl 2kh
Y/3 -2Y“+1+Yi=:-Y:-,
kh2
+St’-Yt?l
ar = aJ+&x.
(14)
For the corrected coefficients repeat steps 1 and 2. The recomputed profile is considered to be the new profilej + 1. 5. Calculate a new value Aj+‘; A’+’ = A’+ k and the corresponding value of 1050
Solution of nonlinear boundary value problems- VII
af+l=
(15)
cd+ kc.
Thus the new profile is completely determined. 4. EXAMPLE
OF APPLICATION
Consider the nonlinear boundary value problem which arises in chemical reaction engineering
g2+fg =6y exp
rP(1 -Y)
Table 2. The effect of the grid spacing on the accuracy of calculated parameter 8. The case without correction
1 (16)
0.05 0.02 0.01 0.005
[ l+P(l-Y)
n=O.dy=() ‘dx
x= l:y=
(17)
1.
Differentiation of Eq. (16) with respect to A yields the following third order partial differential equation
a9
a --a2y
YPY
ax2aA+ iiaxaA 6 [ 1-(l+j3(l-_y))2 YN-Y)
xexp [ l+p(l-y) =yexp
yP(l-Y)
l+p(l-y)
1 1
ay
1
ax
5
10
20
40
80
1.62077 1.69264 1.71812 1.73117
1.61490 168582 1*71095 1.72381
1.61344 168411 l-70916 1.72198
1.61307 168369 1.70871 1.72152
1.61298 168358 l-70860
y(0) = 0.5, a = 0, /3 = 0. The exact value 6 = 1.73438. Table 3. The effect of the grid spacing on the accuracy of calculated parameter 6. The case with correction. y(0) = 0.5, a = 0, 6 = 0. The exact value S = 1.73438.
0.05 0.02 0.01 0.005
5
10
20
40
80
l-74056 1.74378 1a74426 1a74439
1.73310 l-73626 1.73673 1.73685
1.73124 1.73438 1.73485 1.734%
1.73077 1.73391 1.73438 1.73450
1.73066 1.73379 1a73426
z=
d8 dA
(18)
subject to boundary conditions
aY(o4 = 0; y(l,A)
For the nonisothermal case (fl f 0) the results are in Table 4, the exact values are calculated by the GPM technique. Again the correction essentially improves the accuracy attained, the results
Table 4. The effect of the grid spacing on the accuracy of calculated parameter 6. y(O) = 0.5, (I = 0, y = 20, p = 0.05. The exact value 6 = 1.15462. Correction
= 1; y(O,A) =A.
no
cash ti=-
1 Y(O) -
(19)
The results calculated without correction (step 4 in algorithm) are presented in Table 2, the corrected results are in Table 3. It is obvious that the correction yields more accurate results.
N
0.05
(19)
For numerical integration of this equation we must choose an appropriate initial condition, for 6 = 0 a trivial profile y(x, 1) = 1 results. A simple comparison of the numerical accuracy can be carried out for the isothermal case (p = 0) and plate geometry (a = 0) where the analytical solution of (16) and (17) is simple. The dependence between 6 and y(0) can be evaluated from the transcendental equation
t
\ ._
0.02 0.01 0.005 0.05
yes
0.02 0.01 0.005
5
10
20
40
1.12763 1.14437 1.15028 1.15331
1.12659 1.14317 1.14903 1.15202
1.12634 1.14288 1.14872 1.15170
1.12628 1.14281 1.14864 1.15162
l-15542 1.15622 1.15633 1.15636
1.15412 1.15490 1.15501 l-15504
l-15380 1.15458 1.15469 1.15472
1.15372 l-15450 1.15461 1.15464
calculated for the event without correction and halved step k/2, which requires approximately the same computational effort as the algorithm with correction, are worse than that computed with correction. In practice, the relatively sparse grid (iV = 10, k = O-02) and the correction procedure yields satisfactory results. On using this method multiple solutions do not complicate the computation as far as y(x, A) and cw(A)are con-
1051
MILAN
KUBfCEK
and VLADIMfR
Table 5. The effect of the 8rid spacing on the accuracy of calculated parameter 8. The case with multiple solutions. Algorithm with correction.y(0) = 05, a = 2, y = 40, p = 0.2. -. N
k\ 0.05 0.02 0.01 0.005
5 0.406% 0.41242 0.41311 0.41328
0.41279 0.41805 0.41871 0.41887
0.41414 0.41934 0.42000 O-42015
0.41447 O-41%5 O-42031 0.42047
tinuously differentiable. The case which has been solved by the differentiation with respect to an actual parameter with serious difficulties [ 11, can be calculated by this method in the straightforward way. The results for multiple solutions are reported in Table 5. 5. DISCUSSION
AND
CONCLUSION
For a set of differential equation, e.g.
Y” =f(y, T”
= g(y,
Y’, T, T’,
a)
y’,
a)
T, T’,
with two-point boundary conditions the preceding algorithm can be used in an analogous way. Here a seven-diagonal band matrix results. The technique can be considered as an altemative to the differentiation with respect to an actual parameter. For the case of multiple solutions we can take use of this new method, on the other hand, for unique solutions the differentiation with respect to an actual parameter can be adopted. The results of calculation show that a relatively sparse grid of mesh point yields satisfactory accuracy. Thus, for example, branching points can be readily found. Here the parameter C is zero. As previously pointed out the method can be extended to other classes of nonlinear boundary value problems.
REFERENCES [ll KUBfCEKM. and HLAVACEKV., [2] KUBfCEK M. and HLAVACEK V., [31 KUBfCEK M. and HLAVACEK V., [4] KUBfCEK M., Sci. Pap. Inst. Chem.
HLAVACEK
Gem. EngngSci. 197126705. Chem. Engng Sci. 1972 27 743.
Chem. Engng Sci. 1972 27 2095. Technol. Prague, in press.
1052
Science. 1973, Vol. 28. pp. 1049-1052.
Pergmon
Press.
Printed inGreat
Britain
Solution of nonlinear boundary value problems-VII A novel method: differentiation with respect to boundary condition MILAN KUBICEK and VLADIMfR HLAVACEK Department of Chemical Engineering, Institute of Chemical Technology, Suchbatarova 1903, Prague 6, Czechoslovakia (Received 24 August 1972) Abstract-A simple numerical method for construction of the dependence of solutions to nonlinear boundary value problem on a parameter will be developed. The set of differential equations is diiferentiated with respect to the boundary condition chosen and the resulting partial differential equations are solved by a finite-difference method. The method is illustrated by an example of heat and mass transfer in a porous catalyst. 1. INTRODUCTION
Now supposing p0 # 0 we will choose
THERE ARE
few methods available for calculating the dependence of a solution of nonlinear boundary value problem on the given parameters so far[l-41. The basic technique used, general parameter mapping [2,3] and differentiation with respect to the actual parameter[ 1,4] transform the original boundary value problem into some initial value problem; in the first case the necessity of meeting boundary conditions requires solution of an initial value problem with an implicit right-hand side while the other method leads to a third order partial differential equation. The necessity of unique solution is the primary disadvantage of the last technique. For multiple solutions serious difficulties arise. It is the purpose of this article to develop a modified technique which may be readily used for handling problems with multiple solutions. 2. DEVELOPMENT
OF THE METHOD
We will assume, for simplicity, a single differential equation of second order Y” -“f-k
y, Y’, aY) = 0
where f is a nonlinear conditions are abY
function.
+PoY’(a)
%Y (b) +PIY
(3)
=A,
here A will be considered as a new parameter. Furthermore, we will assume that the function f is differentiable with respect to all arguments; on differentiation we get a partial differential equation for y(x, A)
aS_mY+2c?%+~dar_ axZaA_
ayaA ay’axaA
aadA
(4)
For a given A the boundary conditions are
cw(a,A) +Po
aY(GA) ax
=
Yo
(5)
cr,y(b,A)+P,~=y,
(6)
y(u,A)
(7)
= A.
Thus we have three boundary conditions. It seems that the problem is over-determined, however, in the partial differential equation (4) we have one unknown parameter
(1)
do C=J-J.
The boundary
= 3/o
(6) =
Y(U)
(2)
We note that the problem is completely mined. Now on choosing A=Ao:a=c+,,
Yl.
1049
y(x,A,)
=y,,(x)
deter(9)
MILAN
KUBiCEK
and VLADIMiR
HLAVACEK
we can integrate the differential equation (4); here y&) is the solution of Eq. (1) with boundary conditions (2) and an additional condition y,(a) = A. Frequently, the solution y,(x) may be readily calculated, sometimes for a = twoa trivial solution exists. 3. NUMERICAL
SOLUTION
ORDER
j+l
OF THE THIRD
EQUATION
For the numerical solution of Eq. (4) a finite difference method has been used. The approximation of the particular derivatives in the mesh point (i,j) is presented in both Table 1 and Fig. 1. The nonlinear terms @/lay, aflay’, and aflaa have to be computed at the old profile j. It is obvious that an evaluation of the new profile demands the solution of a set of linear equations having special almost tridiagonal matrix. For the new profile we must calculate the values of both y(+’ (i = 0, 1, . . ., N) and C. After dividing the region (a, b) into N parts following difference equations result
i-l
j
:
1
Y
Fig. 1. A sketch of the grid of the mesh points.
Note that the profile j has been calculated for Ai. This procedure yields a set of N+ 4 linear equations (lo)-(13) for N +4 variables yJ_:‘, 5+1 . ., yhyl, C. The set of equations has a 2 YIJ+’ Yo special foh so that following calculation procedure can be used:
1. Guess C = Cl. For i = 0 the values yl_:‘, j+l, Yo ylf+l can be calculated from Eqs. (lo), (11) &y&J + s;y2+1+ sj+,y;;; + rjc = 4: and (13). The recurrence relation (10) yields (10) . ., yi;;. In general, the values yi?i, Yi+l, Ys’+l i=O,l,...,N y,+!+l and >iyl do not meet Eq. (12) and a residue where the coefficients s:_, , s(, di+l, $, 4 me de- R (Cl ) depending on the estimation of C, C = Cl, pendent on the known values j$, = y,,’ (n = i- 1, results. i, i+ 1) and 6 = 01’. For the approximation of 2. Guess C = C, # Cl and repeat step 1. Beboundary conditions (5)-(7) we have cause the residue R(C) is a linear function of C the correct value of C can be calculated by linear Y1J+’ -ti’ = y 0 (11) interpolation. abYd+l+ PO 2h 3. For the correct value of C a new profile j+ 1 must be recalculated or interpolated from E-yKi_ Yl (12) two profiles evaluated for C = Cl and C = C,.. fflYN1+' +Pl y 4. This new profile yif+l may be employed to yoj+1 = Al+ k. (13) correct the coefficients sC1, st, s:+~, riJ, qrJ; these will be calculated for the average value y, deTable 1. Approximation of derivatives in fined as arithmetical mean, i.e. Eq. (4) 2h
-
ay _ Y,‘+‘-Y*’ aA k _apY _ YE: -Ye:
axaA _ a? _ ax2aA
yi = 4(Y2+Y(+‘), -Yd+1 +Ycl 2kh
Y/3 -2Y“+1+Yi=:-Y:-,
kh2
+St’-Yt?l
ar = aJ+&x.
(14)
For the corrected coefficients repeat steps 1 and 2. The recomputed profile is considered to be the new profilej + 1. 5. Calculate a new value Aj+‘; A’+’ = A’+ k and the corresponding value of 1050
Solution of nonlinear boundary value problems- VII
af+l=
(15)
cd+ kc.
Thus the new profile is completely determined. 4. EXAMPLE
OF APPLICATION
Consider the nonlinear boundary value problem which arises in chemical reaction engineering
g2+fg =6y exp
rP(1 -Y)
Table 2. The effect of the grid spacing on the accuracy of calculated parameter 8. The case without correction
1 (16)
0.05 0.02 0.01 0.005
[ l+P(l-Y)
n=O.dy=() ‘dx
x= l:y=
(17)
1.
Differentiation of Eq. (16) with respect to A yields the following third order partial differential equation
a9
a --a2y
YPY
ax2aA+ iiaxaA 6 [ 1-(l+j3(l-_y))2 YN-Y)
xexp [ l+p(l-y) =yexp
yP(l-Y)
l+p(l-y)
1 1
ay
1
ax
5
10
20
40
80
1.62077 1.69264 1.71812 1.73117
1.61490 168582 1*71095 1.72381
1.61344 168411 l-70916 1.72198
1.61307 168369 1.70871 1.72152
1.61298 168358 l-70860
y(0) = 0.5, a = 0, /3 = 0. The exact value 6 = 1.73438. Table 3. The effect of the grid spacing on the accuracy of calculated parameter 6. The case with correction. y(0) = 0.5, a = 0, 6 = 0. The exact value S = 1.73438.
0.05 0.02 0.01 0.005
5
10
20
40
80
l-74056 1.74378 1a74426 1a74439
1.73310 l-73626 1.73673 1.73685
1.73124 1.73438 1.73485 1.734%
1.73077 1.73391 1.73438 1.73450
1.73066 1.73379 1a73426
z=
d8 dA
(18)
subject to boundary conditions
aY(o4 = 0; y(l,A)
For the nonisothermal case (fl f 0) the results are in Table 4, the exact values are calculated by the GPM technique. Again the correction essentially improves the accuracy attained, the results
Table 4. The effect of the grid spacing on the accuracy of calculated parameter 6. y(O) = 0.5, (I = 0, y = 20, p = 0.05. The exact value 6 = 1.15462. Correction
= 1; y(O,A) =A.
no
cash ti=-
1 Y(O) -
(19)
The results calculated without correction (step 4 in algorithm) are presented in Table 2, the corrected results are in Table 3. It is obvious that the correction yields more accurate results.
N
0.05
(19)
For numerical integration of this equation we must choose an appropriate initial condition, for 6 = 0 a trivial profile y(x, 1) = 1 results. A simple comparison of the numerical accuracy can be carried out for the isothermal case (p = 0) and plate geometry (a = 0) where the analytical solution of (16) and (17) is simple. The dependence between 6 and y(0) can be evaluated from the transcendental equation
t
\ ._
0.02 0.01 0.005 0.05
yes
0.02 0.01 0.005
5
10
20
40
1.12763 1.14437 1.15028 1.15331
1.12659 1.14317 1.14903 1.15202
1.12634 1.14288 1.14872 1.15170
1.12628 1.14281 1.14864 1.15162
l-15542 1.15622 1.15633 1.15636
1.15412 1.15490 1.15501 l-15504
l-15380 1.15458 1.15469 1.15472
1.15372 l-15450 1.15461 1.15464
calculated for the event without correction and halved step k/2, which requires approximately the same computational effort as the algorithm with correction, are worse than that computed with correction. In practice, the relatively sparse grid (iV = 10, k = O-02) and the correction procedure yields satisfactory results. On using this method multiple solutions do not complicate the computation as far as y(x, A) and cw(A)are con-
1051
MILAN
KUBfCEK
and VLADIMfR
Table 5. The effect of the 8rid spacing on the accuracy of calculated parameter 8. The case with multiple solutions. Algorithm with correction.y(0) = 05, a = 2, y = 40, p = 0.2. -. N
k\ 0.05 0.02 0.01 0.005
5 0.406% 0.41242 0.41311 0.41328
0.41279 0.41805 0.41871 0.41887
0.41414 0.41934 0.42000 O-42015
0.41447 O-41%5 O-42031 0.42047
tinuously differentiable. The case which has been solved by the differentiation with respect to an actual parameter with serious difficulties [ 11, can be calculated by this method in the straightforward way. The results for multiple solutions are reported in Table 5. 5. DISCUSSION
AND
CONCLUSION
For a set of differential equation, e.g.
Y” =f(y, T”
= g(y,
Y’, T, T’,
a)
y’,
a)
T, T’,
with two-point boundary conditions the preceding algorithm can be used in an analogous way. Here a seven-diagonal band matrix results. The technique can be considered as an altemative to the differentiation with respect to an actual parameter. For the case of multiple solutions we can take use of this new method, on the other hand, for unique solutions the differentiation with respect to an actual parameter can be adopted. The results of calculation show that a relatively sparse grid of mesh point yields satisfactory accuracy. Thus, for example, branching points can be readily found. Here the parameter C is zero. As previously pointed out the method can be extended to other classes of nonlinear boundary value problems.
REFERENCES [ll KUBfCEKM. and HLAVACEKV., [2] KUBfCEK M. and HLAVACEK V., [31 KUBfCEK M. and HLAVACEK V., [4] KUBfCEK M., Sci. Pap. Inst. Chem.
HLAVACEK
Gem. EngngSci. 197126705. Chem. Engng Sci. 1972 27 743.
Chem. Engng Sci. 1972 27 2095. Technol. Prague, in press.
1052