N u m e r . M a t h . 27, 4 4 9 — 4 6 9 (1977)
© b y S p r i n g e r - V e r l a g 1977
Comparing Routines for the Numerical Solution
of Initial Value Problems of Ordinary Differential Equations in Multiple Shooting
H . - J . Diekhoff, P . L o r y , H . J . Oberle, H . - J . Pesch, P . Rentrop, and R . Seydel
R e c e i v e d J u l y 26, 1976
Summary. T h e n u m e r i c a l Solution of t w o - p o i n t b o u n d a r y value problems a n d p r o b l e m s of o p t i m a l c o n t r o l b y shooting techniques requires i n t e g r a t i o n routines.
B y s o l v i n g 15 real-life problems four w e l l - k n o w n integrators are c o m p a r e d r e l a t i v e to r e l i a b i l i t y , fastness a n d precision. H i n t s are g i v e n , w h i c h routines c o u l d be used for a p r o b l e m .
1. Introduction
I n recent years a series of refined numerical methods has been developed for the Solution of i n i t i a l value problems i n ordinary differential equations. The large number of published integration routines necessitates a comparison of a l l these methods. E n r i g h t et al. [46], H u l l [23], and Davenport et al. [10] have tested a n d compared integration routines b y solving carefully selected i n i t i a l value problems.
I n applications many nonlinear two-point boundary value problems and prob- lems of optimal control arise. Convenient methods for the numerical Solution of these problems are shooting algorithms. Since these require a repeated compu- tation of i n i t i a l value problems, integration routines also play an important role.
I n this paper the authors have tested the suitability and the behaviour of certain integrators for a multiple shooting method. T h e y are m a i n l y interested i n a comparison of integration routines, which help a user select a method for a particular problem. Secondly several boundary value problems and o p t i m a l con t r o l problems are proposed and treated, which are realistic and representative of real-life applications. This m a y stimulate further comparisons of algorithms for solving boundary value problems.
I n order to compare the integrators, 15 two-point boundary value problems are solved, which arise i n different application areas: mathematics, physics, engineering, space science and economics. This report is concerned o n l y w i t h non-stiff problems. A l l the selected examples are solved b y means of a multiple shooting technique described, for example, i n Stoer and B u l i r s c h [40]. I n this algorithm a modification of Newton's method due to Deuflhard [11] is used.
It m a y be recalled that i n multiple shooting, sequences of i n i t i a l value prob- lems must be solved w i t h iteratively determined i n i t i a l values (for realistic prob-
lems about 100 i n i t i a l value problems i n every iteration i n v o l v i n g numerical differentiation). Because of the artificial character of these i n i t i a l values the integration routines have to overcome harder difficulties than i n the Solution of natural i n i t i a l value problems. I n addition, the singularities of the Solutions of nonlinear ordinary differential equations depend o n the initial values (movable singularities). That is w h y the integration routines must be reliable, fast a n d precise.
2. Integration Routines
The results of the investigations of 20 routines i n [16] suggested the following four integrators:
D I F S Y l : Bulirsch-Gragg-Stoer extrapolation method.
This is an extrapolation algorithm based on the midpoint rule w i t h Gragg stabili- zation [21]. I t was first published b y B u l i r s c h a n d Stoer [6]. I n Hussels [24] an improved stepsize control was implemented.
V O A S : Adams-Sedgwick method of variable order a n d variable step.
This is an implementation of a variable order variable step A d a m ' s method developed b y Sedgwick [39].
R K F 7 : Runge-Kutta-Fehlberg method of seventh order, R K F 4 : Runge-Kutta-Fehlberg method of fourth order.
These are R u n g e - K u t t a methods w i t h built-in estimators for the local error, see [16]. T h e y are based on the formulas of Orders seven a n d four, respectively, developed b y Fehlberg [18, 19] and E n g l a n d [17].
3. Multiple Shooting Method
A detailed description of the multiple shooting method for the Solution of a two-point boundary value problem
y' = f(t>y)\ y' [*, b]-+KN, /: [a, b] x R " - * R * (2.1.a)
r(y{a),y{b)) = 0; r: R ^ x R ^ I R * (2.1.b) m a y be found e.g. i n [40] a n d [7].
Here the interval [a, b] is suitably subdivided
a = t1<t2... tM_x<tM — b (M nodes).
Denoting y(t; tjt s7), j=\, M— 1 the Solution of the i n i t i a l value problem y'=i(t>y)\ y(tf)=s,-> <e [*y,</ + 1]
the A^-vectors s; have to be determined so, that the following N(M— 1) con- ditions h o l d :
continuity conditions (for M > 2)
Fj(sj9sj+1): =y(tj+1; tjts,)-sj+1 = 0 / = ! , . . . , AT —2 (2.2.a)
boundary conditions
*5if-i(si, SM-I) ' = r(svy{tM; tM_v sM_1)) = 0.
The conditions (2.2) define a System of N(M— 1) nonlinear equations
(2.2.b)
F{s): =
FM__2{SM-2> SM—l) FM-I(SI> SM~l)
-0 w i t h s: = (2-3)
This System is solved numerically b y the modified Newton-method:
sk+1=sk+XkAsk 0<Xhgi, As"=-DF (s*) _ 1.F(s*) {DF(s) denotes the Jacobian matrix).
(2.4.a) (2.4.b) A good strategy for choosing the Xk has been developed b y Deuflhard [12]. Start- ing the process (2.4) the following initial data must be available
6; 7 = 1, M — 1 .
W i t h the following abbreviations dy(ti+x;tj,sf)
A: 8 r
8 s,
B:
7 = 1, . . . , M - 1 ,
(2.4.b) can be written i n more detail (the Ä-index is omitted):
0 TM-2
0
3 As2
A sM_2 _A %/_!_
A V^M-l The zl s; allow a recursive determination b y
A s—G^A sH 1 + Ff^ j=2, . . . , M-1
£ : =^4 + BGM_XGM_2... Gx iteration matrix,
co : = — (FM_1+BGM_1FM_2-\ + BGM_1 ...G2F1).
(2.4.b')
w i t h
(2.5.a) (2.5.b) (2.6.a) (2.6.b) A, B and the G; are computed b y numerical differentiation. This requires the calculation of N trajectories. The use of Broyden's approximation formulas [3]
reduces Computing time up to 30%. It is apparent that the choice of a reliable, fast and precise integration routine plays an important role i n multiple shooting.
Remark. The piecewise constructed function
y(t):=y(t;tjfSi) for *€[*,-, / = ! , . . . . , M - l is called a trajectory. I n particular for s,=S/
y (t) is the starting trajectory.
4. Numerical Examples
T o compare the four integration methods, problems were chosen from the areas of mathematics, physics, technical science, economics and space science.
These examples include boundary value problems, integral equations, eigenvalue problems and problems of optimal control. The problems r ä n g e from simple ones of theoretical nature to highly complex ones.
A l l examples were transformed into two-point boundary value problems (ex- amples 1-11) or into boundary value problems w i t h switching functions (examples 12-15). T h e y were solved b y the Standard multiple shooting algorithm B O U N D -
S O L and its modified version especially trimmed to handle optimal control problems.
The experiments were run on the T R 440 of the Leibniz-Rechenzentrum der Bayerischen Akademie der Wissenschaften. The computations were performed in F O R T R A N single precision with a 38 bit mantissa (examples 1-11) and i n
F O R T R A N double precision w i t h an 84 bit mantissa (examples 12-15).
The statistics of the comparisons include the following parameters:
a) Parameters concerning the problem:
N : N u m b e r of differential equations.
M : N u m b e r of nodes for the multiple shooting method.
N T : N u m b e r of trajectories.
T O L : Tolerance for the i n i t i a l value methods.
O V H B : Overhead time for B O U N D S O L , i.e. the difference between the total time and the time spent i n solving i n i t i a l value problems.
E P S : Relative tolerance for B O U N D S O L ( 1 .1 0— 6 for a l l examples).
b) Parameters concerning the comparison of the methods:
T I M E : T o t a l time to solve the problem.
O V H M : Overhead time for M E T H O D , i.e. the difference between the time solving the initial value problems and the time spent i n evaluating the functions.
N F C : N u m b e r of function calls. N F C depends on the number and position of the nodes. A proper choice of the nodes, depending on a special inte- gration routine, can produce a 30% reduction i n Computing time. A n optimal choice for a l l four routines is not possible because of the dif- ferent kinds of stepsize control.
c) Parameter concerning the final assessment:
N C T : E x a c t time of solving an example b y one method i n per cent of the average time of the four methods.
Because of the different length of the problems a final assessment of the methods and their properties is difficult. T o overcome this difficulty,
" a normed Computing t i m e " N C T is introduced. T h i s s i m p l y calculated number has the advantage of being independent of the Computer used, and of the Computing time. F o r the single examples the N C T is not important, it is only used for the final comparison i n the con- clusion (see Chapter 5).
Example 2. A Standard test problem (N=2).
This linear two-point boundary value problem is discussed i n [40], p. 205 (see also Pereyra [30] and Daniel and M a r t i n [9]):
_y" + 400j> = — 400 cos2 {nt) — 2 n2 cos (2 nt) y(0)^y{\) = 0.
Single shooting technique w i l l fail to solve this problem [40]. Since the Solution has one exponentially increasing and one exponentially decreasing component, the following transformation is introduced:
y1=y
y2 = 20y+y'.
Nodes: 0, 0.25, 0.5, 0.75, 1.
Zero initial data
T I M E O V H M N F C N C T D I F S Y l 2.5 1.1 3901 39 V O A S 4.6 3.7 2222 72 R K F 7 3-2 1.8 4684 50 R K F 4 15-2 8.7 22635 238
Example 2. A singular bifurcation problem (N=2).
— u"+ [ K )3] ' = * u> u{0) = u{i) = 0.
This problem has been studied b y K ü p p e r [27]. It approximates the partial differential equation of the model of a continuous vibrating string considered by F e r m i et al. [20]. A s discussed i n [27], for Xe [8, n2] the problem has a positive Solution which is Symmetrie and bifurcates from the t r i v i a l Solution at X=n2. Outside of [8, n2] the problem has no positive Solution.
The Solution at 2 = 8. is calculated. F o r this parameter value the problem possesses a numerical singularity at t = 0. F o r that reason a series approximation i n the first subinterval is used.
M = 5
N T = 6 T O L = i .1 0- 6 O V H B - 0 . 3
Nodes: equidistant
I n i t i a l data: Solution at 2 = 8.5 (obtained b y continuation) T I M E O V H M N F C N C T
M = 1 0 M = 1 0
D I F S Y l 1.2 0.6 2023 80 N T = 8
V O A S 2.2 1.6 1667 147 T O L = 1 .1 0—6
R K F 7 1.6 0.9 2704 107 O V H B = 0.5 R K F 4 1.0 0.4 1300 67 A = 8.5->A = 8.
Example 3. An example due to Troesch (iV = 3)- The two-point boundary value problem is given b y :
y ^ A s i n h {Xy)t 2 ^ 0 ,
>>(0) = 0, ^ ( l ) = i .
T h i s problem had been discussed b y several authors (e.g. Scott and W a t t s [38], p. 74/75, where further references are listed). It describes the confinement of a plasma column b y radiation pressure. The exact Solution of the problem i n terms of Jacobian elliptic functions can be found i n [40], p. I69. F o r the numerical evaluation of the Jacobian elliptic functions it is recommended to employ the computationally economic, rapidly convergent method of the arithmetic-geo- metric mean [4]. Multiple shooting was used as a background procedure for comparison purposes only.
The problem was solved b y a modified continuation method [13], a t r i v i a l differential equation and a linear boundary condition for the homotopy parameter 2 is added.
The sensitivity of the problem is caused b y a logarithmic singularity at t = tx>\ of the Solution of the associated i n i t i a l value problem. One might expect numerical difficulties for the multiple shooting method because tx depends on the parameter X and approaches 1 + as X increases. Indeed, the norm of the iteration m a t r i x increased from 1.1 00 to 1.1 09 as X varied from 1 to 17.5. There- fore single precision computation had to be terminated.
Nodes: 0, 0.3, 0.5, 0.6, 0.7, 0.8, 0.85, 0.9, 0.92, 0.94, O.96, 0.97,0.98,0.99, 1.
I n i t i a l d a t a : Solution at 2 = 7.25.
T I M E O V H M N F C N C T
M ==15 M ==15
D I F S Y l 3.5 1.7 3393 48 N T = 8
V O A S 11.4 9.0 5224 156 T O L = 1 .1 0-•6
R K F 7 8.9 5 3 9114 122 O V H B - 0 . 8 R K F 4 5.5 3-2 5212 75 A = 7.25->A = 7-5
Remark. F o r small values of X V O A S or R u n g e - K u t t a methods are faster than D I F S Y l , but for all X = 4 D I F S Y l turns out to be best. F i r s t of all, V O A S faüs at 2 = 9-75 then R K F 7 at 2 = 1 1 . 5 and R K F 4 at 2 = 14.0. The Standard multiple shooting algorithm suggests a new node be inserted w h i c h coincides w i t h 1. i n the first four digits. The last successful step for D I F S Y l is 2 = 1 7 . 5 .
Example 4. Artificial boundary layer problem (N=2).
This linear problem was treated i n Daniel and M a r t i n [9], L e n t i n i and Pe- reyra [28], and Pereyra [30]:
jT&r: *<°>=
0' ^ > = V J T O T
At t=0 the problem has a boundary layer of thickness ]/2. One expects serious difficulties for the stepsize control of the four routines as X approaches zero.
Nodes: 0, | ] / l . i0- 7 , 0.01, 0.1.
I n i t i a l d a t a : Solution at X = 1.10— 5.
T I M E O V H M N F C N C T
M = 4 M = 4
D I F S Y l 0.7 0.4 1464 25 N T = 5
V O A S 7.0 6.6 3615 222 T O L = l .1 0- 6
R K F 7 1.4 1.0 2751 44 O V H B = 0 . 2
R K F 4 3.4 3.1 7116 108 2 = 1-io-" 5 - ^ = l .1 0
Remark. Solving this example D I F S Y l is the most reliable and V O A S the most susceptible routine. The last successful steps has been for
V O A S : A = l.1 0- 8 R K F 4 , R K F 7 : 2 = 1 .1 0- 1 0
D I F S Y l : 2 = 1 .1 0- 1 3 .
The stepsize control of the routines Orders a new node which is too close to a given one.
Example 5. An integral equation of the second kind (N=$).
67iq)(s)= f (3 sin r sin s + 2 sin 2s sin 2r) {cp{r)-\-(pz{T))dr.
0
This equation has three different Solutions i n addition to the t r i v i a l Solution (Pimbley [32]). E a c h of them can be obtained b y multiple shooting, after trans- forming the integral equation into a boundary value problem b y :
yi(t) = <p(t), t
y* W = / sin r (99 (T) + (pz ( T ) ) dr,
0
71
yz (t) = j sin r (y (r) + (p3 (r)) dr, t
t
t
31 Numer. Math., B d . 27
yA{t)=J sin (2T) (<p(r) + <f^{r))dr, o
y,(t)= / sin (2T) {<p{x) + ^(r))dx.
Nodes: 0, 0.5, n\29 2, n.
I n i t i a l d a t a : y=2 except for the given boundary values.
U s i n g the chosen data, multiple shooting yields the only positive Solution.
T I M E O V H M N F C N C T
D I F S Y l 7-5 3-5 5361 96 M = 5
V O A S 10.9 8.0 3411 140 N T = 1 8
R K F 7 5.8 3.0 3471 74 T O L = i .1 0- 6
R K F 4 7.0 3.6 4379 90 O V H B = 0 . 8
Example 6. Membrane theory (A7 = 2).
This problem arises i n the theory of stress distribution i n a spherical membrane having normal a n d tangential loads (see L e n t i n i a n d Pereyra [28], Russell and Shampine [36], a n d Scott [37]). The following differential equation holds:
y + (3 cot (0 + 2 tan (*)) / + 0.7y = 0f 3>(30°)=0 >'(60°) = 5.
The sharp rise i n y (t) a n d its derivatives m a y cause difficulties. The problem is scaled b y a factor 1.1 0—4.
Nodes: 30°, 31°, 60°.
I n i t i a l d a t a : y = 0, y' = TT/1800.
T I M E O V H M N F C N C T
D I F S Y l 7.9 4.4 15615 138 M = 3
V O A S 5.3 4.5 3097 93 N T = 6
R K F 7 4.3 2.6 6742 75 T O L - 1 .10
R K F 4 5-4 3-4 8796 94 O V H B = 0. 2
Example 7. Heat conduction (iV = 3).
The equations for one-dimensional heat conduction w i t h nonlinear heat gen- eration are treated i n N a and T a n g [29] a n d Stoer and Bulirsch [40]:
y'(0)=y(i)=0
where X: heat generation constant 0 < X = 0.8, y: temperature distribution.
The problem possesses a numerical singularity at t = Q. I n actual computation the differential equation must be separated (see [40] a n d [34]).
Nodes: 0., 0.1 (matching point), 1.
Zero i n i t i a l data.
T I M E O V H M N F C N C T
M = 3 M = 3
D I F S Y l 1.0 0.5 1573 82 N T = 7
V O A S 1.7 1.3 804 139 T O L = 1.10— 6
R K F 7 1.4 0.8 2132 114 O V H B= 0 . 2
R K F 4 0.8 0.4 1008 65 A=0.8
Example 8. Lubrication theory (N—2).
The following nonlinear eigenvalue problem was treated b y Cole et al. [8] and Keller [25]. It arises i n the theory of lubrication and concerns the flow of a viscous compressible fluid through a very narrow gap.
1 / sin^ t \
y'= — (sin2* - X — ^ — j : —n\2^t^ n\2 y{n\2)=y{-n\2) = \\
s is a parameter, X the required eigenvalue.
This boundary value problem was solved numerically for s = 1/100. F o r this parameter value the Solution y approaches zero at £ = 0.
The Single shooting method was performed i n the backward direction (inte- grating from TT/2 to — JT/2). In this case the norm of the iteration matrix was about 1.10 0. O n the other hand, during forward integration the norm increased to 0.5i010, indicating high sensitivity of the problem.
Nodes: — TZ/2, n\2 {M — 2, Single shooting).
I n i t i a l d a t a : y = 1, X = 4/3.
T I M E O V H M N F C N C T
D I F S Y l 15-4 9.0 22627 125 M = 2
V O A S 7-5 6.2 4017 61 N T = 5
R K F 7 10.6 6.7 13326 86 T O L = 1 .1 0— 6
R K F 4 15-7 10.1 19896 128 O V H B = 0 . 2
Example 9. Nonsymmetric bending of elliptic cylindrical tubes (N=6).
Recently Weinitschke [43, 44] studied the nonlinear problem of bending a thin-walled cylindrical tube w i t h elliptic cross section. N o assumptions are made of bending the tube about the minor or major axis. The differential equations are:
/ T = a2 / j/cos2! + £2 s i n2! ( s i n ß ( c o s ! + ^ sin !) — c o s ß ( y c o s ! — ^ s i n £)) q sin | cos |
— cos2 ! + £2 s i n2! ^ '
/ " = ] / c o s2! + e2s i n2! (cos ß (cos ! +y e sin !) + sin ß (y cos ! — e sin !)) q s i n ! cos !
~~~ c o s2 ! + *2s i n2! '
where:
! : angle 0 ^ ! ^ 2T Z ,
/ ( ! ) : dimensionless stress function, ß (!): dimensionless angular deflection,
b: minor axis, a: major axis,
Rx: projection of the deformed axis to the x — z plane, assumed to be radius of a circle,
Ryi projection to the y—z plane, D: flexural rigidity,
1 JA : extensional stiffness,
e = - , q=\-e\
Rx a*
Ry' fÄDRa = x •
S y m m e t r y of the problem leads to the boundary conditions:
ß(0)=ß(~) = arctany.
/'(0) = 0,/(f)=0.
A d d i t i o n a l l y , two integrals for the moments mx, my are computed:
jt/2
= - j f [((1 + C) cosß + 5 sinß) cos / - e ((1 + C) s i nß- S cos ß) sin t] dt, n
o
Wy = _ j ^ . j f[((i-C)sinß + Scosß)co$t+e((\-C)cosß-Ssinß)smt]dt o
where: S = sin (2 arctan y)t C = cos (2 arctan y).
I n the simpler case of a cylindrical tube, e = i, Solutions have been obtained b y means of a perturbation technique b y Reissner and Weinitschke [33]. Quasi- linearization has been successfully applied b y Thurston [42].
The problem was solved for the parameters:
e = 0.9, y = 0.25, a = 0.5.
Nodes: 0., n\\, n\2.
I n i t i a l d a t a :
ß = arctan y
ß' =7t2l}0 estimated from [33],
/ ( < » — 1 / ( T ) — «
/ ( f ) =
f t/'(0)=0 /'(f) = 2/B /'(f)=2M,
mx = 0,
T I M E O V H M N F C N C T
M = 3 N T = 12 T O L - 1 - 1 0 - 6 O V H B = 0.4
D I F S Y l 2.2 0.8 1317 93
V O A S 2.7 1.7 806 114
R K F 7 2.6 1.2 1387 109
R K F 4 2.0 0.9 1064 84
Example 10. A hydrodynamic problem (N=$).
The problem describes a laminar boundary layer produced b y the rotating flow of a viscous incompressible fluid over a stationary infinite disk. The external velocity should v a r y as some power of the radius: v~r~n; further the fluid is assumed to be interacting w i t h a magnetic field represented b y the parameter 5.
A complete discussion of this problem can be found i n K i n g and Lewellen [26].
The original Version leads to a two-point boundary value problem:
VV>" + n WY- i + &2-sy)' = 0,
~^y>&'+(n-i)y)'<P-s{0-- 1) = 0 boundary conditions:
y)(0)=:y)' (O) = 0(O)=O, y'(*/) = 0, 0(tf) = i where n and s are given parameters.
For several parameters n, s, tf numerical Solutions of this problem were ob- tained b y H o l t [22] using a combination of quasilinearization and finite difference methods a n d b y Roberts and Shipman [35] using a shooting method.
In [35] the Solution is given for the parameters
n=—0A, 5 = 0.2, tf=ii.}.
Nodes: equidistant.
Initial d a t a : Solution at tf=\\.}.
T I M E O V H M N F C N C T
M = 1 1 M = 1 1
D I F S Y l 10.0 6.4 13015 82 N T = 1 2
V O A S 16.7 14.1 7326 138 T O L = l. i o — 6
R K F 7 7-5 5-1 6876 62 O V H B = l. 3
R K F 4 14.3 10.4 15352 118 ^ = 1 1. 3 ^ = 1 2 .
Remark. Changing the direction of integration i n the multiple shooting algo- r i t h m a more favorable iteration matrix was obtained (see [13]).
Example 11. Re-entry problem for an Apollo orbiter type vehicle ( i V=7).
Of a l l the many difficult problems due to manned interplanetary travel per- haps the most critical is that of returning the vehicle from outer space to the earth's surface. The chief problem w i t h aerodynamic breaking is associated w i t h the severe heating effects experienced b y the vehicle. Therefore reentry trajec- tories are searched which minimize this quantity. This optimal control problem
leads to seven differential equations i n v o l v i n g the adjoint variables Xv, Ay, and the free t o t a l flight time T (cf. [40]).
P h y s i c a l differential equations:
SQV ' g' l \ ogs\n y r 2 m L w W ~ (1+f )2
SQV v cos y _ g cos y 2m <4 W + R ( \ * > ( l + f ) 2
(velocity),
T (flight path angle),
(normalized altitude).
H a m i l t o n i a n :
where
Cw(u) CA (u) sinu cosu a o R S/m
= 1.174 —0.9cos u,
= 0.6sinu (u is the control variable),
= - 0 . 6 Ay/a,
= — 0 . 9 ^ Av/ a ,
=209-, £ - 4 . 2 6 , e0= 2 . 7 0 41 0- 3 , g = 3 - 2 1 7 2l 0- 4 ,
= 53 200.
The boundary conditions are given b y : v{0) =0.36 v ( l ) = 0.27,
y(0) = - 8 . 1 ^ / 1 8 0 y(i) = 0.,
| ( 0 ) = 4/i? f ( l ) = 2.5/Ä,
^ U = o . Nodes: I n i t i a l d a t a :
V y T
K K
1. 0 . 2 7 0 0 E + 0 0 0.0 0.11961722488E-01 0. 2 3 0 E + 0 3 - 0 . 6 0 0 0 E + 0 0 0.47 500000000E -f-00 - 0 . 3 5 1 0 0 0 0 0 0 0 0 E + 02
0.762 0.2750E-fOO 0 . 7 0 6 8 5 8 3 4 7 0 5 E - 0 2 0. 1 1 7 2 2 4 8 8 0 3 8 E - 0 1 0. 2 3 0 E + 0 3 - 0. 3 9 0 0 E + 00 0.31500000000E+00 - 0. 6 5 700000000E + 01
0.496 0.2880E + 00 0.24085543677E-01 0.10574162680E-01 0.230E+03
— 0.2280E + 00 0.12600000000E + 00 - 0.34600000000E+01
0.26 0.3300E+00 - 0 . 2 1 2 9 3 0 1 6 8 7 4 E - 0 1 0 . 8 2 7 7 5 H 9 6 1 5 E - 0 2 0.230E+03 - 0 . 1 0 7 0 E + 00 -0.40000000000E - 0 1 0.0
0.097 0.3607E+00 -0.12217304764E + 00 0.14019138756E-01 0 . 2 3 0 E+03 - 0 . 5 0 0 0 E — 0 1 - 0 . 5 5000000000E - 01 —0.13000000000E— O l
0. 0.3600E + 00 -0.14137166941E+00 0.19138755981E—Ol 0.230E + 03
- 0 . 6 0 0 0 E —01 -0.60000000000E — Ol —0.20000000000E+00
T I M E O V H M N F C N C T
D I F S Y l 113.5 39.6 67224 60 M = 6
V O A S 194.4 143.1 46058 104 N T = 100
R K F 7 129.2 58.6 65142 69 T O L = l- i o- ö R K F 4 313.7 137.9 166992 167 O V H B = 3-4
Remark. I n order to prevent an overflow of the right hand side using R K F 7 , the argument of the exponential function must be limited.
Example 12. Optimal control of a mass production (iV = 3).
T h i s optimal control problem is due to Bauer, Neumann [2] and describes an automatic machine producing a mass article:
v: production velocity: v € [0,1 ], OLV : articles produced per unit-time: oc > 0, 1 — x: percentage of refuse: x € [0,1 ],
(xxv: returns of the production per unit-time,
k (v, t): production costs per unit-time: k (v, t) = \ v2 eT . The profit / is to be m a x i m i z e d :
T
I(x,v) = J [<xxv — k(v,t)]dt
0
under the restriction (depreciation):
x— —bxv b>0, x(0) = i
where a = 2., 6 = 0.08 and T= 1 0 .
The associated two-point boundary value problem is given b y p =v(0LX-±ve~T) p(0) = 0,
x=—bxv x(ö) = \,
K=v(K*>-*) **M=o w i t h the control variable
0, if v < 0 (this case does not occur i n actual computation) v, i f i? € [ 0, l ]
1, if 5 > 1 , t
v=e T[x(<x-Xxb)]
switching function: S =v — 1.
Nodes: equidistant.
Zero i n i t i a l data except for the given boundary values.
T I M E O V H M N F C N C T
D I F S Y l 6.1 1.9 2729 49 M = 3
V O A S 7-3 4.6 1420 59 N T = 15
R K F 7 17.0 1 1 3 4389 137 T O L = l - i o - 9
R K F 4 192 12.5 5147 155 O V H B = 0.7
Example 13. Optimal economic planning of a growing nation (N=7).
This example describes the economy of a growing nation which suffers from great unemployment due to the business structure. A policy leading to füll employment and stable growth is required. The data used held for Algeria i n the year 1961.
The economical model due to Stoleru [41] includes one constrained control variable and one first-order constrained State variable. U s i n g multiple shooting technique W i e k [45] has computed the associated two-point boundary value problem. The differential equations which consider two cases are given b y (x€ [0,1 ] is the independent variable):
Casel. O n unconstrained arcs (g : — maeyxT—z < 0 ) : y =xuyTt
z = a ( l — u)y T,
K= — (Xy<x.u + X2<x(l — u)) T,
\z=X=h=f==o.
Gase 2. O n constrained arcs (g = 0) the differential equations for Xy, Xz and Xr change t o :
X=(Xy-Xz)y*zT
w i t h the constants a = 0.25, o = }At A = B = 4, y = 0.125, and w = 0.45. The ad Joint variables Xz and Xr are discontinuous at the entry point xx of the con- strained arc:
X+\x^x=X7\^Xl + hmoye^T. The control variable u satisfies:
u =
if x < xx (entry point); g < 0 on the constrained arc; g = 0
if x > x2 (exit point); g < 0.
The boundary conditions are:
y{0)=i y{\)=AevT, z(0) = a z{\) = BevT, Ay(0) = l ,
#|*=o =
[
—^ + A y j / + A2i +A
rr],
= 0=
0 (Hamiltonian),# U = 0 .
Nodes: I n i t i a l data:
0. y= i. z= 3.1 =0.2 T = 20. /0 = 0.2 K= 1 0.5 15- z= 6. Ay= 0 . 0 5 = 0.005 T=20. l0=0.2 K= 1 1. y= 70. z = 7 0 . Xy=0. K-= 0.005 7 = 2 0 . /0= 0 . 2 K= 1
T I M E O V H M N F C N C T
D I F S Y l 29-0 190 11597 19 M = 3
V O A S 42.4 32.9 4893 27 N T = 28
R K F 7 94.9 84.8 14164 61 T O L = 9
R K F 4 460.6 410.9 74486 294 O V H B = 3.0
Example 14. Planar Earth—Mars transfer (N=9).
I n interplanetary space science i t is required to pilot an i o n rocket w i t h m i n i m a l mass-loss. I n this example a travel from the E a r t h i a n orbit to the M a r - tian orbit is treated. B o t h orbits are assumed to be circular and coplanar. T h e influence of the gravitational field of the E a r t h is neglected. This leads to an optimal control problem w i t h two control variables: the thrust ß € [0, ßm a x] (linear control variable) and the thrust angle ipE] — TC, TZ[ (nonlinear control variable).
F o r the (normed) physical quantities (distance from the sun r, radial-veloc- i t y (o, tangential-velocity vy massw) the four differential equations h o l d :
r=a>,
r r2 r m r cov
c = 1.872, Anax = 0.075, v = — Vß — cos w
m=-ß.
The boundary values are given b y : r(0) = l .
Ö>(0) = 0.
y(0) = l . w (0) = 1.
The problem w i t h free total flight time tf is yet unsolved r(</) = 1.525, o>(t,) = 0,
= 0.8098,
Nodes: equidistant.
Initial d a t a : Solution at t 5.0.
T I M E O V H M N F C N C T
M = 5
M = 5
D I F S Y l 46.4 20.2 12203 53 N T = 2 3
V O A S 54.8 37-9 5870 63 T O L = 1 .l 0- 9
R K F 7 82.2 61.9 8967 94 O V H B = 2 . 3
R K F 4 166.5 122.7 19470 190 // == 5 . 0 - ^ = 5 . 0 5
Example 15. Heating constraint crossrange maximization problem for a Space Shuttle orbiter—type vehicle (N—14).
Since the development of the Space Shuttle atmospheric entry of lifting ve- hicles is studied. These lifting vehicles should be capable of considerable lateral r ä n g e , which allows increased return frequencies from orbit to given landing sites.
F o r m a x i m u m lateral r ä n g e a heating constraint depending on velocity, altitude, and angle of attack is taken into account.
The underlying mathematical model of this o p t i m a l control problem is due to Dickmanns [14], who also prepared initial data for the unconstrained problem (no heating constraint) using an analytical approximate Solution. I n [31] the numerical Solution of the constrained problem has been treated using multiple shooting techniques especially trimmed to handle ill-conditioned two-point bound- ary value problems (cf. [11]). A n extensive presentation of the whole model would be beyond the scope of this paper, i t can be found i n [13]. The heating constraint was pushed down from 2 8 5 0 ° F to 1 7 0 0 ° F b y means of a continuation method. T h e special homotopy step reduces the level of the permitted skin temperature of the Space Shuttle from 2057 ° F to 2000 °F. Using the modified continuation method [13] the System comprises 14 differential equations.
The differential equations for the physical quantities (velocity v, heading angle %, flight path angle y, cross-range angle A, altitude h, down-range angle 6) are:
X=CA
' = CA
2m
at, s i n u e ßhv —
cos y 2m e ßhvcosju-
R+h go l R
cos y cos % tan A ,
\2 V R+hj ~ R+h cosy 3
h=vsmy,
6 cos ycos %
R+h cos A
where [i (aerodynamic bank angle) a n d cA (lift coefficient, constrained) are the two nonlinear control variables. The other quantities are constants. F o r a detailed discussion of the technical results see [15].
This example is one of the most complicated and most sensitive control Prob- lems to be solved so far b y multiple shooting algorithm.
Nodes and i n i t i a l data can be found i n [31].
T I M E O V H M N F C N C T
D I F S Y l 3850.6 795.9 346370 35 M = 11
V O A S 3071.9 1672.7 146689 28 N T = 95
R K F 7 7252.9 3970.7 384485 66 T O L = 1- i o -
R K F 4 F A I L 272 O V H B = 24.7
Remark. R K F 4 d i d not succeed i n 10800 sec (number of computed trajec- tories: N T= 3 0 , number of function evaluations: N C F = 582396). T o make the calculation of N C T for R K F 4 possible, the time was extrapolated to 30000.
5. Summary of Results
The 15 examples presented i n this paper may be considered as representative two-point boundary value problems. They cover different degrees of difficulties.
Some of the examples are rather simple and have been treated repeatedly i n the literature (expl. 1, 3, 4, 6), others are extremely difficult (expl. 11, 13, 14, 15).
The criteria for judging the four Integration routines are reliability, fastness and precision.
a) Reliability. A s far as reliability is concerned, D I F S Y l obtained the best results. D I F S Y l is not sensitive to the choice of nodes (expl. 3,4). It also runs well even i n sensitive problems described b y a large norm of the iteration m a t r i x
(cf. Chapter 3). A disadvantage of D I F S Y l is the increased Computing time, if D I F S Y l is used w i t h zero initial data i n the case of a homogeneous differential equation.
A s can be seen from the examples (3,4), V O A S was the most sensitive inte- grator. It does not allow large Integration intervals, and requires the greatest number of nodes. Nevertheless, V O A S proved to be more favourable i n appli- cations i n optimal control problems. The small integration steps of V O A S yielded a secure and fast determination of the switching points.
R K F 7 and R K F 4 are very reliable, except that the stepsize control of R K F 7 d i d not avoid an overflow i n the right hand side of example 11. B o t h of these integrators d i d not allow as large parameter values as D I F S Y l d i d i n the ex- amples 3 and 4.
b) Fastness. I n order to compare the fastness of the four routines, the total amount of Computing time is considered. The number of function calls, which has been measured i n [28] and [30], is not alone representative for fastness. Com- p u t i n g time must be paid, not the number of function calls. E v e n for complicated problems w i t h a large number of Operations i n the right hand side, the over- headtime of the routines can not be neglected (see Table 1).
A s is seen i n Table 1, the mean of the 0verheadtime for all problems (in per cent) requires a larger portion of time than the time spent i n the right hand side. O n l y i n example 15 T F C is dominant i n a l l routines. It should be mentioned that V O A S
Table l . Splitting of the total Computing time (100%) averaged over all examples D I F S Y l V O A S R K F 7 R K F 4
O V H M 47% 75% 61% 60%
T F C 40% 18% 2 9 % 30%
O V H B 13% 7% 10% 10%
O V H M == percentage for the overheadtime of the routines; T F C = percentage for the time spent i n the right hand side; O V H B = percentage for the overheadtime of Boundsol
Table 2. Averaged normed Computing time i n per cent of the R K F 7-average E x p l . l - i i E x p l . 12-15 E x p l . 1-15
D I F S Y l 95 44 81
V O A S 152 49 123
R K F 7 100 100 100
R K F 4 135 254 169
needs the largest, D I F S Y l the smallest part of O V H M . T h e overheadtime of B O U N D S O L , however, is negligible.
It is not possible to compare the four routines b y adding the total Computing times for each problem. The various difficult ies and sizes of the problems caused Computing times ranging between a few seconds (expl. 1) a n d several hours (expl. 15). The normed Computing times N C T introduced i n C h a p t e r 4 are com- parable quantities.
Table 2 includes the averaged N C T for each routine given i n per cent of the R K F 7-average i n each column.
This table illustrates the fastness of each routine for the examples 1-11 w i t h technical precision ( T O L = 1 .1 0— 6 ) , a n d for the sensitive o p t i m a l control prob- lems 12-15 w i t h higher precision ( T O L = 1.10—9) • Column three i n Table 2 contains the averaged values over a l l problems. A s is seen i n examples 1-11 D I F S Y l , and R K F 7 r u n similarly fast, whereas i n examples 12-15 D I F S Y l a n d V O A S are equivalent. R K F 4 is out of business. H a v i n g considered a l l problems, one would most likely prefer D I F S Y l . The above division of the examples into two groups, however, is unfavourable for R K F 4 , since R K F 4 was the fastest routine i n the examples 2, 7, 9.
c) Precision. I n shooting methods, the precision is prescribed ' a p r i o r i ' . A l l methods obtained the required precision.
Rule of Thumb. I n order to decide which routine should be applied for a par- ticular problem, the 15 examples are divided into four groups (according to their averaged total Computing time).
Simple examples: 2, 7, 9
Less difficult examples: 1, 3, 4, 5, 6, 8 Difficult examples: 10, 12, 14
Very difficult examples: 1 1,13, 1 5 .
F o r these four groups the averaged normed Computing times for each routine are compared i n Table 3. A g a i n they are given i n per cent of the R K F 7 - a v e r a g e i n each column.
Table 3. Averaged normed Computing time for the Classification i n per cent of the R K F 7-average
Simple Less difficult Difficult Very difficult
D I F S Y l 77 104 63 58
V O A S 121 165 89 81
R K F 7 100 100 100 100
R K F 4 65 163 158 374
F o r simple problems (column 1) R K F 4 is the fastest routine, whereas for less difficult problems (column 2) one would prefer R K F 7 or D I F S Y l . I n difficult and very difficult problems D I F S Y l and V O A S are the fastest routines. A m o n g very difficult problems, where the part of Computing time spent i n the right hand side dominates (see expl. 15), V O A S is faster than D I F S Y l . T h i s is caused b y the larger percentage of overheadtime of V O A S (see Table 1) (i.e. the smaller number of function calls). These Statements are summed up i n the Figure 1 w h i c h also includes the experience of the authors.
Good behaviour
Bad behaviour
VOAS,
— — 'DIFSYl
Simple
problems Complicated problems
F i g . 1. B e h a v i o u r of the i n t e g r a t i o n routines
The r e s u l t s of these i n v e s t i g a t i o n s v a l i d a t e t h e suggestions i n [16], t h a t t h e f o u r t e s t e d r o u t i n e s s h o u l d be i n c l u d e d i n a p r o g r a m l i b r a r y . The c h o i c e of a s u i t a b l e r o u t i n e for a p a r t i c u l a r p r o b l e m is m a d e easier b y t h e s k e t c h .
Acknowledgements. T h e authors w i s h to t h a n k Professor R . B u l i r s c h w h o s t i m u l a t e d a n d encouraged t h i s w o r k , t h e y are indebted to D r . D . Gerber for helpful comments.
References
1. A z i z , A . K . (ed.): N u m e r i c a l Solutions of b o u n d a r y v a l u e problems for o r d i n a r y differential equations. Proceedings of a S y m p o s i u m h e l d at the U n i v e r s i t y of M a r y l a n d , 1974. N e w Y o r k : A c a d e m i c Press 1975
2. B a u e r , H . , N e u m a n n , K . : B e r e c h n u n g o p t i m a l e r Steuerungen, M a x i m u m p r i n z i p u n d d y n a m i s c h e O p t i m i e r u n g . L e c t u r e Notes i n Operations Research a n d M a t h e - m a t i c a l Systems, V o l . 17. B e r l i n - H e i d e l b e r g - N e w Y o r k : Springer 1969
3. B r o y d e n , C . G . : A class of methods for s o l v i n g non-linear simultaneous equations.
M a t h . C o m p . 19, 577-593 (1965)
4. B u l i r s c h , R . : N u m e r i c a l c a l c u l a t i o n of e l l i p t i c integrals a n d e l l i p t i c functions, I, I I . N u m e r . M a t h . 7, 78-90, 353-354 (1965)
5. B u l i r s c h , R . , O e t t l i , W . , Stoer, J . (eds.): O p t i m i z a t i o n a n d o p t i m a l c o n t r o l . P r o - ceedings of a Conference h e l d at O b e r w o l f ach, 1974, L e c t u r e Notes, V o l . 477.
B e r l i n - H e i d e l b e r g - N e w Y o r k : S p r i n g e r 1975
6. B u l i r s c h , R . , Stoer, J . : N u m e r i c a l t r e a t m e n t of o r d i n a r y differential equations b y e x t r a p o l a t i o n methods. N u m e r . M a t h . 8, 1-13 (1966)
7. B u l i r s c h , R . , Stoer, J . , D e u f l h a r d , P . : N u m e r i c a l Solution of nonlinear two-point b o u n d a r y v a l u e problems, I. N u m e r . M a t h , (to be published)
8. Cole, J . D . , K e l l e r , H . B . , Saffmann, P . G . : T h e flow a viscous compressible f l u i d t h r o u g h a v e r y n a r r o w gap. S I A M J . A p p l . M a t h . 15, 605-617 (1967)
9. D a n i e l , J . W . , M a r t i n , A . J . : I m p l e m e n t i n g deferred corrections for N u m e r o v ' s difference m e t h o d for second-order t w o - p o i n t b o u n d a r y v a l u e problems. T e c h . R e p . C N A- 1 0 7 , Center for N u m e r i c a l A n a l y s i s , T h e U n i v e r s i t y of T e x a s at A u s t i n (1975) 10. D a v e n p o r t , S. M . , S h a m p i n e , L . F . , W a t t s , H . A . : C o m p a r i s o n of some codes for the i n i t i a l v a l u e p r o b l e m for o r d i n a r y differential equations. I n : [1], p p . 349-353 (1975)
11. D e u f l h a r d , P . : A m o d i f i e d N e w t o n m e t h o d for the Solution of ill-conditioned Systems of n o n l i n e a r equations w i t h a p p l i c a t i o n to m u l t i p l e shooting. N u m e r . M a t h . 22, 289-315 (1974)
12. D e u f l h a r d , P . : A r e l a x a t i o n s t r a t e g y for the modified N e w t o n m e t h o d . I n [5], p p . 59-73 (1975)
13- D e u f l h a r d , P . , Pesch, H . - J . , R e n t r o p , P . : A modified c o n t i n u a t i o n m e t h o d for the n u m e r i c a l Solution of n o n l i n e a r two-point b o u n d a r y value problems b y shooting techniques. N u m e r . M a t h . 26, 327-343 (1976)
14. D i c k m a n n s , E . D . : M a x i m u m r ä n g e , three-dimensional l i f t i n g p l a n e t a r y e n t r y . N A S A T e c h . R e p . R - M 199, M a r s h a l l Space F l i g h t Center, A l a b a m a (1972) 15. D i c k m a n n s , E . D . , Pesch, H . - J . : Influence of a r e r a d i a t i v e heating c o n s t r a i n t on
l i f t i n g e n t r y trajectories for m a x i m u m l a t e r a l r ä n g e . 11 t h I n t e r n a t i o n a l S y m - p o s i u m o n Space T e c h n o l o g y a n d Science, T o k y o , J u l y 1975
16. E n r i g h t , W . H . , B e d e t , R . , F a r k a s , I., H u l l , T . E . : Test results o n i n i t i a l v a l u e methods for non-stiff o r d i n a r y differential equations. T e c h . R e p . N o . 68, D e p a r t - m e n t of C o m p u t e r Science, U n i v e r s i t y of T o r o n t o (1974)
17. E n g l a n d , R . : E r r o r estimates for R u n g e - K u t t a t y p e Solutions to Systems of o r d i n a r y differential equations. C o m p . J . 12, 166-170 (1969)
18. F e h l b e r g , E . : K l a s s i s c h e R u n g e - K u t t a - F o r m e l n f ü n f t e r u n d siebenter O r d n u n g m i t S c h r i t t w e i t e n k o n t r o l l e . C o m p u t i n g 4, 93-106 (1969)
19. F e h l b e r g , E . : K l a s s i s c h e R u n g e - K u t t a - F o r m e l n f ü n f t e r u n d siebenter O r d n u n g m i t S c h r i t t w e i t e n k o n t r o l l e u n d ihre A n w e n d u n g auf W ä r m e l e i t u n g s p r o b l e m e . C o m p u t i n g 6, 61-71 (1970)
20. F e r m i , E . , P a s t a , J . R . , U l a m , S . : Studies of nonlinear problems, I. T e c h . R e p . N o . 1940, L o s A l a m o s (1955)
21. G r a g g , W . B . : O n e x t r a p o l a t i o n a l g o r i t h m s for o r d i n a r y i n i t i a l value problems.
S I A M J . N u m e r . A n a l . Ser. B 2, 384-403 (1965)
22. H o l t , J . F . : N u m e r i c a l Solution of n o n l i n e a r two-point b o u n d a r y v a l u e problems b y finite difference methods. C o m m . A C M . 7, 366-373 (1964)
23. H u l l , T . E . : N u m e r i c a l Solutions of i n i t i a l v a l u e problems for o r d i n a r y differential equations. I n [1], p p . 3-26 (1975)
24. Hussels, H . G . : S c h r i t t w e i t e n s t e u e r u n g bei der I n t e g r a t i o n g e w ö h n l i c h e r D i f - ferentialgleichungen m i t E x t r a p o l a t i o n s v e r f a h r e n . U n i v e r s i t ä t z u K ö l n , M a t h e - mathisches I n s t i t u t , D i p l o m a r b e i t (1973)
25. K e l l e r , H . B . : N u m e r i c a l m e t h o d s for t w o - p o i n t b o u n d a r y v a l u e problems.
C h a p t e r 6: P r a c t i c a l examples a n d c o m p u t a t i o n a l exercises. L o n d o n : B l a i s d e l l 1968
26. K i n g , W . S., L e w e l l e n , W . S.: B o u n d a r y - l a y e r s i m i l a r i t y Solutions for r o t a t i n g flows w i t h a n d w i t h o u t magnetic interaction. Ref. A T N- 6 3 (9227)-6, A e r o d y n a m i c s a n d P r o p u l s i o n Res. L a b . , Aerospace C o r p . , E l Segundo, Calif. (1963)
27. K ü p p e r , T . : A singular bifurcation p r o b l e m . M a t h . R e p . N o . 99, B a t t e l l e A d v a n c e d Studies Center, G e n e v a (1976)
28. L e n t i n i , M . , Pereyra, V . : B o u n d a r y p r o b l e m solvers for first order Systems based on deferred corrections. I n [1], p p . 293-315 (1975)
29. N a , T . Y . , T a n g , S. C . : A method for the Solution of c o n d u c t i o n heat transfer w i t h non-linear heat generation. Z A M M 49, 45-52 (1969)
30. P e r e y r a , V . : H i g h order finite difference Solution of differential equations. R e p o r t S T A N - C S-73-348, C o m p u t e r Science D e p a r t m e n t , Stanford U n i v e r s i t y (1973) 31. Pesch, H . - J . : Numerische Berechnung o p t i m a l e r Steuerungen m i t H i l f e der M e h r -
zielmethode dokumentiert a m P r o b l e m der R ü c k f ü h r u n g eines R a u m g l e i t e r s unter B e r ü c k s i c h t i g u n g v o n Aufheizungsbegrenzungen. U n i v e r s i t ä t zu K ö l n , M a t h e - matisches I n s t i t u t , D i p l o m a r b e i t (1973)
32. P i m b l e y , G . : E i g e n f u n c t i o n branches of nonlinear Operators a n d their bifur- cations. L e c t u r e Notes, V o l . 104. B e r l i n - H e i d e l b e r g - N e w Y o r k : Springer 1969 33. Reissner, E . , W e i n i t s c h k e , H . J . : F i n i t e pure bending of c i r c u l a r c y l i n d r i c a l tubes.
Q u a r t e r l y A p p l . M a t h . 20, 305-319 (1963)
34. R e n t r o p , P . : N u m e r i c a l Solution of the singular G i n z b u r g - L a n d a u equations b y multiple shooting. C o m p u t i n g 16, 61-67 (1976)
35- Roberts, S. M . , S h i p m a n , J . S.: T w o - p o i n t b o u n d a r y value p r o b l e m s : Shooting methods. N e w Y o r k : E l s e v i e r 1972
36. R ü s s e l , R . D . , Shampine, L . F . : A collocation m e t h o d for b o u n d a r y value prob- lems. N u m e r . M a t h . 19, 1-28 (1972)
37. Scott, M . R . : I n v a r i a n t i m b e d d i n g a n d its applications to o r d i n a r y differential equations: A n i n t r o d u c t i o n . R e a d i n g , Massachusetts: A d d i s o n - W e s l e y 1973 38. Scott, M . R . , W a t t s , H . A . : S U P O R T - A Computer code for two-point b o u n d a r y -
v a l u e problems viaorthonormalization. T e c h . R e p . S A N D-75-0198. S a n d i a L a b o - ratories, A l b u q u e r q u e , N e w M e x i c o (1975)
39. Sedgwick, A . E . : A n effective variable order variable step A d a m s method. P h . D . Thesis, T e c h . R e p . N o . 53, D e p a r t m e n t of C o m p u t e r Science, U n i v e r s i t y of T o r o n t o (1973)
40. Stoer, J . , B u l i r s c h , R . : E i n f ü h r u n g i n die Numerische M a t h e m a t i k , I I . H e i d e l - berger Taschenbuch, B d . 114. B e r l i n - H e i d e l b e r g - N e w Y o r k : Springer 1973 41. Stoleru, L . G . : A n o p t i m a l p o l i c y for economic g r o w t h . E c o n o m e t r i c a 33, 321-348
(1965)
42. T h u r s t o n , G . A . : N e w t o n ' s method applied to problems i n nonlinear mechanics.
J . A p p l . M e c h . 32, 383-388 (1965)
43. W e i n i t s c h k e , H . J . : D i e S t a b i l i t ä t elliptischer Zylinderschalen bei reiner B i e g u n g . Z A M M 50, 411-422 (1970)
44. W e i n i t s c h k e , H . J . : O n nonsymmetric bending i n s t a b i l i t y of e l l i p t i c c y l i n d r i c a l tubes. U n p u b l i s h e d manuscript, T U B e r l i n (to appear)
45. W i e k , R . : Numerische L ö s u n g volkswirtschaftlicher V a r i a t i o n s p r o b l e m e m i t Z u - s t a n d s b e s c h r ä n k u n g e n unter A n w e n d u n g der Mehrzielmethode. U n i v e r s i t ä t zu K ö l n , Mathematisches I n s t i t u t , D i p l o m a r b e i t (1973)
H . - J . D i e k h o f f P . L o r y
H . J . Oberle H . - J . Pesch P . R e n t r o p R . Seydel
I n s t i t u t 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 A r c i s s t r . 21
D-8000 M ü n c h e n 2
F e d e r a l R e p u b l i c of G e r m a n y
