Numerische Mathematik
© by Springer-Verlag 1980
Enlarging the Domain of Convergence
for Multiple Shooting by the Homotopy Method
Peter L o r y
Institut für M a t h e m a t i k der Technischen U n i v e r s i t ä t M ü n c h e n , Arcisstr. 21, D-8000 M ü n c h e n 2, G e r m a n y (Fed. Rep.)
Summary, The homotopy method is a frequently used technique i n overcom- ing the local convergence nature of multiple shooting. In this paper sufficient conditions are given that guarantee the homotopy process to be feasible. T h e results are applicable to a class of two-point boundary value problems.
F i n a l l y , the numerical solution of two practical problems arising in physiolo- gy is described.
Subject Classifications: A M S ( M O S ) : 6 5 L 1 0 ; C R : 5.17.
1. Introduction
M u l t i p l e shooting is a w e l l - k n o w n method for the numerical solution of nonlinear two-point boundary value problems. It has been thoroughly tested in numerous realistic applications (see e.g. Bulirsch [ 2 ] , Stoer and Bulirsch [25], D i e k h o f f e t a l . [ 8 ] , and K e l l e r [14]). D u e to the nonlinearity of the problem, starting values sufficiently close to the true solution have to be available in order to start the process. A characteristic local convergence theorem is given i n Weiss [28].
The present paper deals with a frequently used way in overcoming the local convergence nature o f the iterative process, namely with the homotopy method (or continuation method). T h i s technique takes advantage of the fact that most of the boundary value problems arising in applications depend o n a physical parameter i n a natural way. In general, this parameter t may appear i n the system of n ordinary differential equations
?=f{T,x9y); xefabl r e [0,1] ( l . l . a )
and i n the two-point boundary conditions
r(T,y(a\y(b)) = 0. (1.1.b) Here (1.1) is chosen so that for i = l (1.1) is equivalent to the original problem.
M o s t l y , (1.1) reduces to a somehow " f a m i l i a r " p r o b l e m for a certain value of the
0029-599X/80/0035/0231/S02.00
parameter, say T = 0. T h e n a solution y(0){x) of (1.1) for i = 0 is available. Starting with T0 = 0, a homotopy chain of subproblems is solved taking the computed solution for T = T¡__Í as an initial a p p r o x i m a t i o n for the ith subproblem. If Tn = 1 is reached (within a tolerable number N of homotopy steps), the original problem is solved.
The general idea of using the homotopy method i n actual computations of nonlinear equations seems to date back to the work of Lahaye [16]. Some authors (e.g. Roberts and S h i p m a n [22], W a c k e r [26]) present theoretical estimates of constant stepsizes J T = T¿ — TF_ p which guarantee the homotopy process to be successful. However, these estimates require the knowledge of computationally unavailable quantities and are by far too pessimistic in those parts of the t-interval where the solution depends on T in an undramatic manner.
H a v i n g an efficient stepsize c o n t r o l l i n g algorithm available it is sufficient to prove the pure existence of a partition 0 = T0 < i j <... < TN — 1 so that the computed solution of the (i — l )t h subproblem is contained in the d o m a i n of attraction for the next one. F o r nonlinear equations in R " this has been carried out by Ortega and R h e i n b o l d t [20], A v i l a [1], M e n z e l and Schwetlick [19]. The present paper treats the h o m o t o p y method for the problem (1.1) in connection with multiple shooting and gives sufficient conditions for its feasibility in terms of the original boundary value problem.
After some preliminaries in Sect. 2, the feasibility of the homotopy method is studied in Sect. 3. The results obtained there may be applied to a class of two- point boundary value problems (Sect. 4). This class contains as a special case a boundary value problem, which is used i n physiology to model salt and water transport i n epithelia. It is closely related to an extensive model of the renal counterflow system. B o t h problems were solved numerically, as described in Sect. 5. In these computations the stepsize control due to Deuflhard [6] was used successfully.
2. Preliminaries
The multiple shooting algorithm for the numerical solution of a two-point boundary value problem (1.1) is described in [25, 14] and realized i n the code in Bulirsch, Stoer, and Deuflhard [4]. Here the interval [a, 6] is suitably sub- divided
a = x{<x2< ..• <xm__ ¡ <xm = b (m nodes).
Let y(Tyx;xk,sk) (k = 1, . . . , m - 1) denote the solution of the initial value problem
y ' = f ( T , X , y ), y (xk) = sk, x € [ xf c, xf c + l] . (2.1) The n-vectors sk have to be determined so that the following n(m — 1) conditions
are satisfied:
continuity conditions (for m > 2 ) :
Fk(^ sk,sk+1) = = y ( t , xk+ ! ; xk, sk) - sk+ ! = 0, k = 1,..., m - 2;
boundary conditions:
i m - i ( ^si ^ m- i ) : = ' - ( T , s1, y ( T , xm; xm_1, sm_1) ) - 0 .
These conditions define a system of n ( m - l ) nonlinear equations Fl(x,sl,s2)
f ( T ^ U , 2 < ^ J " 0 W l t h S = ( / 1 ( 1 2 ) If
y
0 )(x)
is a solution of (1.1) for x = 0 and 40 ) ==y
( 0 )(x
J l),
the n ( m - l ) - v e c t o r s(0) satisfies F(0, .) = 0. Starting with this initial value, the homotopy F defines the following continuation process: Let the partition0 = T0< T ! < . . . < TN = 1 (2.3)
suitably subdivide the interval [0,1]. Then for each subproblem F ( T , . , . ) = 0 M , Newton-iterates are computed (according to the use of Newton's method i n [4, 14], and [25]):
sij+ i =si'J-lDsF{Ti,si'J)yl F(ThsiJ) j = 0,...,M-l i = l , . . . , i V - l .
(2.4) The starting values for these iterations are given by
s1>° = s<°\ s^^o = si,Mi ( 2 5 a )
("classical predictor") or
(2.5.b) s » • o = s<°) - T , • [ Ds F(0, s( 0 )) ] - 1 Dx F(0, s(0>)
s' > » - o = s« ' . " . - ( Ti.+ 1- Ti) . [ DsF ( Ti, si'A ,' ) ] -1ßti ' ( t1- , si'M' ) ("Euler predictor").
(2.4) and (2.5) yield an sN' °, which is used as starting value for the final iteration
¡¡NJ+ i =s"-J-[DsF(l,sN'Jy}-1 F(l,sN'j)
, = 0 , 1 , 2 , . . . . ( 1 6 )
O b v i o u s l y the following definition according to [20] and [1] is meaningful.
Definition. If a partition (2.3) and N-l integers M1,...9MN_i exist so that the process (2.4), (2.5.a) ((2.5.b) respectively) is well-defined and so that (2.6) con- verges to a solution of F(1,.) = 0, then the homotopy method (2.3)-(2.6) is called feasible with the classical predictor (the Euler predictor respectively).
Remark. The Euler predictor (2.5.b) represents an Euler step for the integration of D a v i d e n k o ' s differential equation (Davidenko [5]). A t a first glance, the simplest continuation procedure seems to be integrating this differential equa-
tion from T = 0 to T = 1. However, this approach w o u l d require the frequent evaluation of an explicitly not available right-hand side (cf. [6] and Feilmeier [11]). This argument and a comparison of the c o m p u t i n g times i n Feilmeier and W a c k e r [12] show that this method is very uneconomical.
3. The Feasibility of the Homotopy Method
In this section sufficient conditions for the feasibility of the homotopy method are given. T h e proof requires the following
Lemma. Let a: [0,1] - > R "( m^ l ) be a continuous function with <T(0) = S( O ). Define the stripe
S(Ö,<T): = { ( T , 5 ) | T e [ 0 , l ] , 5 G i ? "( w U, \\S-G(T)\\<0}
with some norm ||. || in Rn i m-l\ Assume that the homotopy F(i,s).
F: S(<9, o) Rn<<m-^ satisfies the following conditions:
F ( T , < X ( T ) ) = 0 for i e [ 0 , l ] . (3.1)
The functional matrix DSF{z,s) exists on S(&,a) and is continuous there.
DSF{T,(J(T)) is nonsingular for t e [ 0 , 1 ] . (3.2)
Then the homotopy method (2.3)-(2.6) is feasible with the classical predictor.
If additionally DTF(x,s) exists on S(<9, a) and is continuous there, then the homotopy method is feasible with the Euler predictor.
The proof slightly extends that of Theorem 10.4.2 in [20] and is omitted here.
It may be found in detail i n [17]. •
Let the functions / : [0,1 ] x [a, è ] x G -> Rn and r : [ 0 , l ] x G x G - > i ? " of the boundary value problem (1.1) be continuous. Here G may denote an open, convex and bounded subset of Rn and G its closure. Then
g(x) = y(x) - / ( T , x, y (x)), 4 ^ w = r(T,y(a),y(b))
defines a h o m o t o p y H ( i , y)
H: [ 0 , l ] x D - > Co[ a , f t ] x Äw (3.4.b) with D: = {yeC{ [a,b~]\y(x)eG for every x e [ a , 6 ] } .
A s usual CY[a,b] is equipped with the n o r m \\y\\{- = m a x ( | | y | |0, | | / | |0) , where II. II0 is the m a x i m u m n o r m .
Theorem. Suppose that in addition to the above assumptions the partial derivative Dyf{x,x,y) exists on [0,1] x \_a,b~] x G and is continuous there and that the same is valid for Du r ( i , w, v) and Dv r(r, u, v) on [0,1] x G x G . Assume that
if y0eD is a solution of i f ( i0, .) = 0, then the linearized problem y'(x) - Dyf{z0, x, y0(x)) • y (x) = 0
£>„ ' f r o , )>o(0), )><>(&)) * J>(*) + Dv r{T09y0(a), y0{b)) • y(b) - 0 ftas on/y í/ie trivial solution y = 0;
# ( T , . ) = 0 has no solution on the boundary of D for i e [ 0 , 1 ] , (3.6) Then, if y{0)eD is a solution of i f (0, .) = 0, there exists one and only one continuous
function Y]\ [ 0 , 1 ] with j/(0) = / 0 ) and if(i,*/(T)) = 0 for T E [ 0 , 1 ] . Additionally, the homotopy method (2.3)-(2.6) is feasible with the classical predictor and with s( 0 )= ( yo ,( x i ) , . . . , y, o )( xB 1.1) )r.
Sketch of the proof (for details see [17]). The implicit function theorem a n d (3.5) guarantee that a solution fj(z) with i f ( T , Í J Í ( T ) ) = 0 for x e [ 0 , a ] ( 0 ^ a < l ) may be continued locally. T h e so-called continuation property (see Rheinboldt [21]) can be proved by the Arzelà-Ascoli theorem and (3.6). Then a standard technique along the lines of Theorem 2.4 i n [21] demonstrates the first part of the statement.
In order to prove the second part, the above lemma is applied. T o this purpose the continuous function a\ [0,1] - >Rn { m~X ) is defined by
/ ax{x) \ / (r¡(T))(Xl) \
Evidently a(0) — s{0\ a n d a detailed investigation shows that J F ( T , S ) (see (2.2)) is defined o n a stripe S(0,<r). (3.1) is obvious and (3.2) follows from the fact, that the solution of an initial value problem depends o n its initial values i n a smooth manner (e.g. Walter [27]).
In order to demonstrate (3.3), let T = T0 and £ : =A + BGm_ i. . . Gi where
Gk: = DSky(T0,xk+x;xk,sk)
with sk: = ak(T0), fc = l , . . . , m - l ;
^ ^ ^ r i T o ^ ^ T o ) , ^ ^ ) ) K ( T ) - ( ^ ( T ) ) ( f e ) ) ;
B: = Dvr(T0,ax{z0l am{T0)).
The chain rule yields
E = A + B • Ds i y(t0, b; a, a{(i0)) and this matrix is proved to be nonsingular i n [28].
A s d e t ( £ ) = d e t ( DsF ( i0, ö r ( T0) ) ) (see [25]), (3.3) is shown. •
Corollary. Suppose that in addition to the hypotheses of the theorem Dzf(z,x,y) exists on [0,1] x [a,b] x G and is continuous there and that the same is valid for Dtr ( i , w, v) on [0,1] x G x G. Then the conclusions of the theorem are right for the Euler predictor, too.
Remark. F o r the sake of brevity the functions / , Dyf and DJ were assumed to be continuous i n x. A more detailed investigation shows that i n the above theorem and i n the corollary piecewise continuity would be sufficient (see [17]).
4. Application to a Class of Boundary Value Problems
In order to illustrate the results of the preceding section, consider the following nonlinear boundary value problem, which is used i n physiology to model salt and water transport i n epithelia (see Sect. 5):
-D-C"(x) + {v(x)-C{x))'=f(x\
v'{x) + 1 • J ( x , C(x)) = 0, 0 < x < 1, (4.1 ) C'(0) = 0, C ( l ) = û, v(0) = 0.
Here D and a are positive constants and the parameter r varies in [ 0 , 1 ] . The function f(x) is assumed to be continuous and nonnegative. J ( x , C) is con- tinuous and continuously differentiable with respect to the second variable and Jc( x , C ) < 0 .
T o apply the theorem, note that
C0 (x) = a + D 1 • \j\ (t) du v0(x) = 0, lj\ (x): = ]f(t)dt)
x \ 0 '
is the only solution of (4.1) for T = 0.
(3.5) is trivial for To = 0 and is demonstrated i n K e l l o g g [15] for T0> 0 . F o r the proof of (3.6), let CT(x), tT(x) be a solution of (4.1), where the subscript indicates the dependence on the parameter r. Then it may be shown analogously to [15] that 0 < CT( x ) < K with a constant K , which is independent of T . A S a consequence,
- T • ]j(t, 0) dt < vz(x) < - 1 • J J(u K)dt
0 0
and | i ?t( x ) | < M3 with a suitably defined constant M 3 . F i n a l l y , there exists a constant M 2 with \C'z(x)\<M2. T o see this, integrate the first differential equation i n (4.1) to obtain
Q x ) = D -,. [ Ct( x ) . i ;tM - /1( x ) ] . Defining
G : = { y G / ?3| 0 <y i< K , | y2| < M2, L v 3 | < M 3 } (3.6) follows immediately.
Conclusion. A p p l y i n g the homotopy method to the boundary value problem (4.1), the feasibility with both the classical and the Euler predictor is guaranteed.
5. Numerical Examples Arising in Physiology
The following two-point boundary value problems were solved by the program B S H O M . T h i s algorithm is based on the multiple shooting code i n [4] and is equipped with an automatic control of the homotopy stepsizes due to [6]. In
addition to the solution it computes a special norm of the matrix E (see Sect. 3).
This n o r m ( £ ) represents the sensitivity of the problem relative to the variation of si. The integration of the initial value problems (2.1) was performed by the routines D I F S Y 1 (Bulirsch and Stoer [3], stepsize control: Hussels [13]) and R K F 7 (Fehlberg [10]). A l l these programs were run on the computers of the Leibniz-Rechenzentrum der Bayerischen A k a d e m i e der Wissenschaften.
Example 1. Water and Solute Transport in Epithelia
M o s t epithelia absorb or secrete specific fluids, such as bile, gastric juice, sweat, and saliva. A t the ultrastructural level, they possess long, narrow channels open at one end and closed at the other. A c c o r d i n g to D i a m o n d and Bossert [7] the water and solute transport in such a channel may be described by the following boundary value problem:
- D • C"(x) + (v{x) • C(x))' = (2/r) • iV(x),
v\x) + {2/r)-P-(C0-C(x)) = 09 (5.1)
C'(0) = 0, C(L) = C0, v{0) = 0.
The dependent variables C(x) and v(x) represent the concentration and the velocity of the fluid flow at height x in a cylindrical channel of length L and radius r. The concentration C0 outside the channel, the water permeability P , and the diffusion coefficient D are positive constants. The rate of active solute transport N(x) is a nonnegative step function:
N0 for 0 ^ x < L / 1 0 , 0 for L / 1 0 ^ x ^ L .
Substituting T • P for P with r e [ 0 , 1 ] the water permeability is varied as a natural parameter. Considering the remark at the end of Sect. 3 (5.1) is contained i n the class of boundary value problems that was investigated i n Sect. 4. Therefore, the homotopy method is feasible.
The computations were performed on a T R 4 4 0 i n single precision arith- metic. The physiological values for N0 and P cover the following ranges:
l .1 0- l o g A r0^ l.1 0- 5 , 1 . 1 0- 6^ i ^ 2 .1 0-4.
In [7] solutions are computed for P = 2 .1 0_5 with N 0 = l .1 0-5, Co = 0.3, D
= l.10- 5 , L = 0.01, r = 5 .1 0-6. In the present paper P is raised to 2 .1 0-4. A s a
result the sensitivity of the problem increases significantly as indicated in the table.
Table 1. C o m p u t a t i o n of (5.1) by the homotopy method using the classical predictor
N H S = number of homotopy steps
^Tm a x / ^Tm i n= : r a ti0 ° f t h e maximal to the m i n i m a l homotopy stepsize
NHS ATmJAxmin n o r m ( £ )
2. i Q 5 5 11 8. ! 04
5. j 0 - s 7 23 2. J Q6
7 64 8.1 07
2. J Q - 4 10 149 6M 09
Remark. Performing the shooting method in backward direction n o r m ( £ ) may be reduced to 2 .1 0 4 (for P = 2 .1 0-4) , whereas the computing time remains nearly unchanged.
Example 2. Kidney Model
A c c o r d i n g zu Stephenson, Tewarson, and M e j i a [24] and Stephenson, M e j i a , and Tewarson [23] the solute and water movement in the kidney may be described by an extensive boundary value problem with the following differen- tial equations:
dFJdx=-Jiv, ] , 6
dFik/dx=-Jik, \ k = {2^ <5'2) Fik = Fiv-Cik~-Dik'dcik/dxJ
where
4 = hiv ' i(c6x-cn) + ( c6 2 - ci 2) ] ,
Jik = hik ' (<",•* - c*k) + aiJ( 1 + bjcikl
J3v = h3v ' KC0 I ~ C3 1 ) + (C02 ~~ C32Ï]*
J3k=h3k'(C3k-COk) + a3k/(i+bJc3k),
k = J Ik + ^2 k
+ 4
it + ^5 fc *This problem seems to be too complicated for a rigorous application of the theory of Sect. 3. However, (5.1) may be regarded as a strongly reduced version of (5.2). So it appears reasonable to attack the kidney model with the same homotopy that was used in E x a m p l e 1. Therefore, the three water permeability coefficients hl v, h3v, and h4v were chosen as natural homotopy parameters.
The numerical solution was performed on a C D C C y b e r 175 (double precision) using the parameter set given in F a r a h z a d and Tewarson [9]. T a k i n g advantage of some trivial parameters the number of differential equations may be reduced to 13 (see [17] and [18]). The problem was solved in 31 (classical) h o m o t o p y steps from permeabilities (0.,0.,0.) up to (10., 1., 10.). In the course of the computations, n o r m (£) increased from 3 .1 0 7 to 4 .1 0,2 indicating an extreme sensitivity of the problem. The homotopy stepsizes were spread by a factor of
about 530. The prescribed relative accuracy was e p s =1 0-5 for the intermediate subproblems and eps = 1 0_1 2 for the final problem. A s the algorithm performed an additional N e w t o n step, after it had reached this accuracy, the continuity and boundary conditions were satisfied with a precision of 1 0.3i for the final results.
A n important check is the conservation of mass which was computed to h o l d to a relative accuracy of 1 0.2 4. Further details are to be found i n [17] and [18].
Acknowledgement. The author wishes to thank R . Bulirsch who stimulated and encouraged this work.
References
1. A v i l a , J . H . : The feasibility of continuation methods for nonlinear equations. SI A M J. N u m e r . A n a l . 11, 102-122 (1974)
2. Bulirsch, R . : D i e Mehrzielmethode zur numerischen L ö s u n g von nichtlinearen Randwertproble- men und Aufgaben der optimalen Steuerung. Report der Carl-Cranz-Gesellschaft, 1971
3. Bulirsch, R., Stoer, L : N u m e r i c a l treatment of ordinary differential equations by extrapolation methods. N u m e r . M a t h . 8, 1-13 (1966)
4. Bulirsch, R., Stoer, J., Deuflhard, P.: (1980) N u m e r i c a l solution of nonlinear two-point boundary value problems I. Numer. M a t h . , H a n d b o o k Series A p p r o x i m a t i o n (in preparation)
5. D a v i d e n k o , D . : O n a new method of numerical solution of systems of nonlinear equations. D o k l . A k a d . N a u k S S S R 88, 601-602 (1953)
6. Deuflhard, P.: A stepsize control for continuation methods and its special application to multiple shooting techniques. Numer. M a t h . 33, 115-146 (1979)
7. D i a m o n d , J . M . , Bossert, W . H . : Standing-gradient osmotic flow. J. G e n . Physiol. 50, 2061-2083 (1967)
8. Diekhoff, H.-J., L o r y , P., Oberle, H J . , Pesch, H.-J., Rentrop, P., Seydel, R . : C o m p a r i n g routines for the numerical solution of initial value problems of ordinary differential equations i n multiple shooting. N u m e r . M a t h . 27, 449-469 (1977)
9. Farahzad, P., Tewarson, R . P . : A n efficient numerical method for solving the differential equation of renal counterflow systems. Comput. B i o l . M e d . 8, 57-64 (1978)
10. Fehlberg, E . : Klassische R u n g e - K u t t a F o r m e l n fünfter und siebenter O r d n u n g mit Schrittweiten- kontrolle. C o m p u t i n g 4, 93-106 (1969)
11. Feilmeier, M . : Numerische Aspekte bei der Einbettung nichtlinearer Probleme. C o m p u t i n g 9, 355-364 (1972)
12. Feilmeier, M . , Wacker, H J . : Z u r numerischen Praxis v o n Einbettungsmethoden. Z . Angew.
M a t h . M e c h . 52, T 2 0 0 - T 2 0 2 (1972)
13. Hussels, H . G . : Schrittweitensteuerung bei der Integration g e w ö h n l i c h e r Differentialgleichungen mit Extrapolationsverfahren. U n i v e r s i t ä t zu K ö l n , Mathematisches Institut, Diplomarbeit, 1973 14. Keller, H . B . : N u m e r i c a l methods for two-point boundary value problems. L o n d o n : Blaisdell
1968
15. K e l l o g g , R . B . : Uniqueness i n the Schauder fixed point theorem. Proc. A m e r . M a t h . Soc. 60, 207- 210 (1976)
16. Lahaye, E . : U n e m é t h o d e de résolution d'une catégorie d ' é q u a t i o n s transcendantes. C . R . A c a d . Sei. Paris 198, 1840-1842 (1934)
17. L o r y , P.: Homotopieverfahren und Anwendung der Mehrzielmethode auf mathematische M o - delle aus der Physiologie. Technische U n i v e r s i t ä t M ü n c h e n , Institut für Mathematik, Disserta- tion, 1978
18. L o r y , P.: N u m e r i c a l solution of a kidney model by multiple shooting. M a t h . Biosci. (in press, 1980)
19. M e n z e l , R., Schwetlick, H . : Z u r L ö s u n g p a r a m e t e r a b h ä n g i g e r nichtlinearer Gleichungen mit s i n g u l ä r e n Jacobi-Matrizen. Numer. M a t h . 30, 65-79 (1978)
20. Ortega, J . M . , Rheinboldt, W . C . : Iterative solution of nonlinear equations i n several variables.
N e w Y o r k : A c a d e m i c Press 1970
21. Rheinboldt, W . C . : L o c a l m a p p i n g relations and global implicit function theorems. Trans. A m e r . M a t h . Soc. 138, 183-198 (1969)
22. Roberts, S . M . , Shipman, J.S.: T w o - p o i n t boundary value problems: Shooting methods. N e w Y o r k : Elsevier 1972
23. Stephenson, J . L . , M e j i a , R., Tewarson, R . P . : M o d e l of solute and water movement in the kidney.
Proc. Nat. A c a d . Sei. U S A 73, 252-256 (1976)
24. Stephenson, J . L . , Tewarson, R . P . , M e j i a , R . : Quantitative analysis of mass and energy balance in non-ideal models of the renal counterflow system. Proc. Nat. A c a d . Sei. U S A 71, 1618-1622 (1974)
25. Stoer, J., Bulirsch, R . : Introduction to numerical analysis. Berlin-Heidelberg-New Y o r k : Springer 1980
26. Wacker, H J . : E i n Iterationsverfahren zur L ö s u n g spezieller nichtlinearer Randwertprobleme.
C o m p u t i n g 9, 275-291 (1972)
27. Walter, W . : G e w ö h n l i c h e Differentialgleichungen. Heidelberger Taschenbuch, Band 110. Berlin- Heidelberg-New Y o r k : Springer 1976
28. Weiss, R . : The convergence of shooting methods. B I T 13, 470-475 (1973)
Received December 21, 1979
