- Email: [email protected]
Computing eigenvalues of Sturm-Liouville problems via optimal control theory
Mathl. Comput. Modelling Vol. 19, No. 10, pp. l-10, 1994 Copyright@1994 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0895-7177194 $7.00 + 0.00
08957177(94)00065-4
Computing Eigenvalues of Sturm-Liouville Problems via Optimal Control Theory C. J. GOH AND K. L. TEO Department of Mathematics The University of Western Australia, Australia R. P. AGARWAL Department of Mathematics National University of Singapore, Singapore (Received
and accepted
February 1994)
Abstract-In this paper, we shall address three problems arising in the computation of eigenvalues of Sturm-Liouville boundary value problems. We first consider a well-posed Sturm-Liouville problem with discrete and distinct spectrum. For this problem, we shall show that the eigenvalues can be computed by solving for the zeros of the boundary condition at the terminal point as a function of the eigenvalue. In the second problem, we shall consider the case where some coefficients and parameters in the differential equation are continuously adjustable. For this, the eigenvalues can be optimized with respect to these adjustable coefficients and parameters by reformulating the problem as a combined optimal control and optimal parameter selection problem. Subsequently, these optimized eigenvalues can be computed by using an existing optimal control software, MISER. The last problem extends the first to nonstandard boundary conditions such as periodic or interrelated boundary conditions. To illustrate the efficiency and the versatility of the proposed methods, several nontrivial numerical examples are included.
Keywords-Sturm-Liouville problems, Eigenvalues computation, Optimal parameter selection, Combined optimal control and optimal parameter selection, Optimal control software.
1. INTRODUCTION Finding eigenvalues and eigenfunctions of a regular Sturm-Liouville problem is an important question in mathematical physics and engineering. In fact, for the computation of eigenvalues of such problems, several successful algorithms are available, e.g., [l-lo]. In Section 2 of this paper, we shall show that the problem of finding the eigenvalues of regular Sturm-Liouville problems can be readily solved by computing for the zeros of the boundary condition at the terminal point as a function of the eigenvalue. Thus, any numerical scheme for solving the zeros of nonlinear functions can be used for this purpose. It is also known that gradient type methods are more efficient than those which require only the function evaluation, an example of this being the Newton method. Note that though the boundary condition is regarded as a function of the eigenvalue, A, the dependence on X is nonexplicit through the differential equation, and thus, the gradient has to be computed in a slightly round-about way. Nevertheless, such gradient formula can be easily derived, and hence, render the proposed approach of solving for eigenvalues a highly efficient and accurate one. It is also free of truncation error associated with some finite difference and finite element schemes. Typ=et by AM-W 1
2
C. J. GOH et al.
In Section 3, we shall consider the case where some coefficients and parameters in the differential equation are continuously adjustable. In this case, the eigenvalues are clearly dependent on how these coefficients and parameters are adjusted, the scope for optimization arises naturally. Analytical solution to this class of problems does exist occasionally, but more often than not, they admit only numerical solution. To compute for numerical solution, we reformulate the problem as a combined optimal control and optimal parameter selection problem. Subsequently, these optimized eigenvalues can be computed by using an existing optimal control software MISER [ll]. To illustrate the efficiency of the proposed method, we consider the celebrated problem of Keller [12]. In brief, the problem is to find the cross-sectional area of a column of constant length and volume such that the Euler buckling load (i.e., the lowest eigenvalue) is maximized. Finally, in Section 4, we shall extend the class of fixed boundary value problems of Section 2 to problems with periodic, semiperiodic or interrelated boundary conditions. Apart from additional constraints, the eigenvalues are readily computable as in the simpler cases.
2. REGULAR
STURM-LIOUVILLE
PROBLEMS
Consider a regular Sturm-Liouville boundary value problem
(p(t)y’)’+ q(t)y + WGY = 0,
(2.1)
UOY(O> + ~lY'(O> = 0,
u;+u::# 0,
(2.2a)
boy(l) + hY’(l)
b; + bT: # 0,
(2.2b)
= 0,
where q, T E C[O, l], p E C(l)[O, 11, and p(t) > 0, r(t) 2 0 in [O,l]. It is well known that for the problem (2.1), (2.2) there exists a countably infinite number of eigenvalues {X, : n = 1,2, . . . }. These eigenvalues can be arranged as a monotonic increasing sequence Xi < X2 < .+’ such that X, + 00 as n --+ 00. Furthermore, the eigenfunction &(t) (unique up to a multiplicative constant) corresponding to the nth eigenvalue X, has exactly (n- 1) zeros in the open interval (0,l). Setting y = xi and p(t)y’ = x2, the Sturm-Liouville problem can be expressed as (2.3a)
with the initial conditions cy and /3 being constrained by aoa + boxi
alp
+ bl &z(l) P(l)
=
0,
= 0,
u;+uf# 0, b; + b; # 0.
and
(2.4a) (2.4b)
Note that the eigenvalues of the problem (2.3) remain the same for different choices of the constants a and ,B as long as they are chosen such that the constraints (2.4a) are satisfied. Let us consider the case where al = 0 in (2.4a). This implies that (Y is necessarily 0, and p must be necessarily nonzero, since otherwise the problem admits only the trivial solution. Without loss of generality, we may arbitrarily choose p to be any nonzero value, say unity. This will also fix the corresponding eigenfunctions since the arbitrary multiplicative constant is prescribed through fixing p. Similarly, if ue = 0, p is necessarily zero and Q may be normalized to unity by the above argument. If we regard X as a parameter, the eigenvalues of (2.3), (2.4) are obviously the zeros of the terminal boundary condition g(X) = b,,xi(l)
1 + bl-x2(1) P(l)
= 0.
(2.5)
Sturm-LiouvilleProblems
3
It is clear that g(X) is a continuous function of X. This is, in essence, a variation of the shooting methods, e.g., [1,6,7]. It is well known (e.g., [S]) that there are countably infinite number of distinct zeros for the function g(X). To compute these zeros, one may use any of the existing numerical algorithms. It is known that algorithms that require gradient information are in general more efficient than those which require only function evaluation. Nevertheless, since the dependence of g on X is nonexplicit, one has to compute the gradient of g(X) in a somewhat round-about way, e.g., [1,6,7]. As this will come under the context of an optimal parameter selection problem, we shall summarize the specific gradient formula for the problem as follows; and leave the definition of an optimal parameter selection problem and the details of its gradient calculation to the Appendix. For a given X, the gradient of g(X) can be computed as follows: Solve the differential equations (2.3) from t = 0 to t = 1. Then, solve the costate differential equations $4(t) = Q(w2(t)
+
WW2($
(2.6a)
backward from t = 1 to t = 0. Finally, calculate the gradient of g(X) according to the formula g’(X) = 1’
-r(t)a(wz(qdt.
(2.7)
0
With this gradient formula, we can use any gradient type method (such as Newton’s method) to calculate the zeros of the function g(X). We note that this gradient formula is different from the usual one used for the shooting methods [1,6,7]. Alternatively, we may also reformulate the problem of finding zeros of g(X) to one that searches for local minima of (g(X))2. This is, in effect, an optimal parameter selection problem as discussed in the Appendix. To illustrate the efficiency of the proposed approach, six examples are solved using MISER [ll]. This is a FORTRAN program for solving constrained optimal control problem, which includes optimal parameter selection problems as special cases. The computed results for the first five eigenvalues of each problem are tabulated in Table 1. Some remarks should be made about the choice of the initial guess and search range for computing the particular kth eigenvalue &. Apart from being continuous and bounded, g(X) is also approximately periodic in 6. To substantiate this, we take the special case of p(t) = r(t) = 1, q(t) = 0, a0 = bo = 1, ai = br = 0. It can be easily shown that, by normalizing the initial condition to y’(O) = fi, the exact expression of g(X) = y(l) is g(X) = sin 45.
(2.3)
For the general problem, the known asymptotic formula
(2.9) can be used to establish a good initial guess xi as well as search range [&, A,] for the lc-th eigenvahie &, once &_.I is computed. This is done as follows:
Ak-1,
(2.10) (2.11) (2.12)
4
C. 3. GOH et al. Table 1. Six examples for eigenvaluea calculation.
Problem
Order A,
-
Initial
Number of
Guess
Iterations
A’
(I) p(t) = 1
1
40
q(t) = 0
2
181
215.76950 514.35637
r (t) = t(1 -t)
-
45.24288
3
485
a0 = bo = 1
4
914
940.97002
01 = bl = 0
5
1470
1495.59860
(II) PV) = 1
1
30
20.20543
q(t) = 0
2
80
87.85666
3
200
203.23746
ao=b=l
4
300
366.36550
al = bl = 0
5
450
577.21640
1
9
2 3
,lo=bo=l
21 = bl = 0
r(t) = sin
-
t
W) P(t) = 1 q(t) = 0 1*(t)
= cash
t
1 -
2
8.59065
26
2
33.89836
70
2
76.07736
4
150
3
135.10375
5
250
3
211.05166
W) p(t) = 1
1
9
2
7.63809
q(t) = 0
2
25
2
30.14484
r (t) = 1 +
3
65
2
67.60399
zQ=bo=l
4
115
2
120.04479
l1 = bl = 0
5
165
3
187.46795
30.45825
t2
1 1 -
W) 1
20
3
2
100
2
96.03750
(4 = *
3
200
1
199.71217
1;ro=bo=l
4
300
3
343.78217
5
450
3
528.78879
P(t) = 1
1
12.7
2
11.92434
Q(t) = -sint
2
45
2
47.33529
I
r(t) = cost
3
90
2
106.38230
1@l=bi)=l
4
160
3
189.09010
5
260
3
295.41930
P(t) = 1 ) = -(l+
qtt l.
-
1l1
t)4
= bl = 0 WI)
~1 = bl = 0
1
For the first eigenvalue, the available bound
(2.13) can be used to establish a good initial guess and search range for X1. Once this is obtained, the initial guess and search range for higher order eigenvalues can be computed recursively through equations (2.10)-(2.12). Note that though higher order bounds similar to (2.13) exist, they are of little help since they may be overlapping.
Sturm-Liouville
Problems
5
Unlike the traditional methods, the proposed method generates the eigenfunctions (unique up to a multiplicative constant) simultaneously with the eigenvalues through the solutions of the differential equation (2.3). Once the calculation has converged, the order of the eigenvalue can be confirmed easily by checking the number of zeros of the corresponding eigenfunction in the open interval (0,l). In obtaining the results of Table 1, we deliberately allow large variation in the initial guess of Xz (as oppose to what is suggested by (2.8)) to demonstrate the versatility and robustness of the method. In all six cases, no more than two or three iterations are required for convergence in spite of the remote initial guess used.
3. OPTIMIZATION
OF EIGENVALUES
Consider the generalized boundary value problem involving adjustable coefficients and parameters in its differential equation
Mt, 4% 5)Y')' + 4 (t74th 5)Y + XT-(t, u(t) ,5)Y = 0,
(3.la)
ooK)Y(O) + M)Y’(O)
= 0,
(3.lb)
bo(0Y(I)
= 0.
(3Sc)
+ WY’(I)
The functions p, q, T, a~, al, bo and br are assumed to satisfy the usual regularity conditions. A control function u( .) and a parameter < are to be chosen from U x R such that the lowest eigenvalue is maximized, where U denotes the class of all piecewise continuous real-valued functions defined in [0,11. As in the previous section, we rewrite (3.1) in the first order system form
(3.2a)
while the boundary conditions (3.lb),
(3.1~) are constrained by
91 = oo(<)o + Ul(S)P = 0,
(3.3a)
gz = boWi
(3.3b)
+
Note that CYand /I are arbitrary as long as they are chosen such that the boundary conditions (3.3) are satisfied without ending up with the trivial solution. As it was explained before, the nonunique eigenfunction may be fixed by normalizing either CYor p. Now the problem of maximizing the lowest eigenvalue may be posed as the following combined optimal parameter selection and optimal control problems: Subject to the system of differential system (3.2) with the boundary conditions (3.3), find a parameter < and a control u E U such that the cost functional go = --A
(3.4)
is minimized. To prevent X from converging to a higher order values, a suitable upper bound to this lowest eigenvalue should be imposed. Note also that here X is regarded as a decision parameter in addition to c. Sometimes, further restraint on the control function u(e) such as the isoperimetric constraint may often be required 1 u(t) dt = ii = constant. (3.5) I0
C. J. GOH et al.
6
The previous assumption that u(.) is piecewise continuous can be easily extended to functions of arbitrary smoothness. To solve the eigenvalue minimization problem, the optimal control software MISER [ll] is used. The underlying concept used is that of control parameterization [13-151. Essentially, the interval [O,l] is partitioned into a number of equal sub-interval, and a piecewise constant function defined on these sub-intervals is used to approximate the function u(.). As such, the optimal control problem is reduced to an optimal parameter selection problem with the constant function value in these subintervals as decision parameters. Alternatively, the eigenvalue minimization problem can also be solved by the gradient-restoration method [16-201. If r~(.) is required to be arbitrarily smooth, v additional differential equations can be appended to (3.2) via xi =
x4
x3(0)
x& =
x5
x4(0)
. . . x’l+v xi+,
=
. . . X2$-Y
=
v(t)
=
u
= il = . . .
~1+,(0)
i2
(3.6)
= i-l
x2+v(O)
= 6v,
where x2+i
(i-l)
i=1,2
I
)...,
Y.
(3.7)
It is clear that if v(t) is allowed to be piecewise constant, u(e) will be in C(“-‘)[O, 11. The initial conditions ii, i = 1,. . . , v, introduced in (3.6) are also regarded as decision parameters along E, a, P and A. To illustrate the applicability of the proposed approach, we column problem of Keller [12]. In brief, this problem considers volume. The aim is to find the variable cross sectional area (A(t) load, i.e., the lowest eigenvalue, is maximized. Mathematically, that the lowest eigenvalue X of the boundary value problem
(AW2y” + Xy
= 0,
t E
consider the celebrated strongest a column of constant length and say,) such that the Euler buckling the problem is to find A(t), such
[O,11,
Y(0) = Y(l) = 0
(3.8a) (3.8b)
is maximized. In this problem, we have normalized the length of the column to 1, and the isoperimetric constraint on the cross sectional area is given by 1
J
A(t) dt = vor
0
(3.9)
where vs is the volume of the column which may also be normalized to 1. Since A(t) is adjustable, it may be regarded as the control. In [12], Keller has shown that the analytical solution of optimum ;\ and A(t) are 4 A* = 3 sin2 0,
8 E [0, ~1,
t=i(,-isin20),
X’ = $2
2 13.1595.
(3.10) (3.11) (3.12)
Using the proposed approach, the strongest column problem of Keller may be restated as (3.13)
Sturm-Liouville
Problems
7
subject to the system of differential equations
x; = x2,
Xl(O)
x x; = -_-2x1,
=
(3.14a)
0,
Q(O) = 1 (normalized),
(3.14b)
x3(0)
(3.14c)
x3 xg = v,
= <,
and the constraints (3.15a)
0 5 A I 2& 91 =
Xl(l)
(3.15b)
= 0,
1 Q2 =
I0
x3(t)
dt - 1 = 0.
(3.15c)
Note that by choosing piecewise constant function for v, the cross-sectional area A(t) = x3(t) is, therefore, piecewise linear. To solve this combined optimal control and optimal parameter selection problem by MISER, we need to supply the number Np of partitions of the interval [O,l], and the initial guess for the control v and the parameters X and E. Using Np = 20, co = 1, v”(t) = 0 (and, hence, x!(t) = AO(t) = 1) and X0 = 10, the solution converges after 25 iterations to X* = 13.1551 (cf. exact value is 13.1595). The corresponding A*(t)
was found to be almost identical to the exact solution of Keller given
by (3.10), (3.11). As a straight forward extension to the above, we modify the boundary condition for (3.8b) to y’(0) = Y’(l) = 0.
(3.16)
It is easy to show that the solution for this modified problem yields the same X* as before, as confirmed by numerical calculation. However, the optimum A*(t) differs from the previous case (i.e., A*(t)) by a shifted symmetry
a*(t)
=
The above boundary conditions correspond to the hinged-hinged and clamped-clamped cases. Other variations of boundary condition such as hinged-free, hinged-clamped, clamped-free, etc., may also be easily computed. Interior point constraints, such as an internal hinge, can also be incorporated with slight modification.
4. NONSEPERATED
BOUNDARY
CONDITIONS
In the literature, boundary conditions of the Consider the differential equation (2.1). type (2.2a)-(2.2b) are referred to as separated boundary conditions. Problems with linear nonseparated boundary conditions of the type aov(O)
+ NY’(o)
= boy(l)
are slightly more difficult to solve in practice. conditions Y(0) = Y(I)
or
+ hY’(l)
(4.1)
This includes the class of periodic boundary Y’(O) = y’(l),
(4.2)
C. J. GOH et al.
8
or semi-periodic
boundary conditions
Y(O)= -Y(l)
or
y’(0) = -y’(l).
(4.3)
Correction techniques associated with both finite element and finite difference methods are available to improve on the accuracy of the computed eigenvalues [2,6-81. In this section, we consider a class of even more general nonseparated boundary conditions of the type fl2i (Y(O), Y’(O), Y(l)? Y’(i)) = 0,
i = 1,2,
(4.4)
where 0i are, in general, nonlinear functions of the arguments. Such nonlinear, nonseparated boundary conditions will unlikely be accessible to many of the existing methods in the literature. Nevertheless, it presents no added difficulty to our present formulation. As before, we assign the initial condition y(O) and g’(O) with suitable decision parameters cx and p, respectively. The nonseparated boundary condition (4.4) will thus reduce to the terminal state constraint G
Ccr, 07 51(l)r
22(1))
=
07
i =
1,2.
(4.5)
By arbitrarily fixing one of CYand p such that the resulting eigenfunction is nontrivial, (4.5) effectively reduces to two simultaneous nonlinear equation in two unknowns X and cy (or ,Q, Gradient formulae can be easily derived using the result of the Appendix, and hence, Newton’s method can again be used to solve for the zeros. The similar argument will still be valid if the differential equation involves continuously adjustable coefficients and parameters as in (3.la). For illustration, an example taken from [16] is computed using the present approach y” - 107r2 cos(27rt)y + 7r2xy = 0,
(4.6)
with the periodic boundary conditions Y(O)
= Y(l),
(4.7a)
Y’(O)
= Y’(l).
(4.7b)
Here, p(t) = 1, q(t) = -10n2cos(27rt), r(t) = x2. The eigenvalues &, k = 2,3,. . . , 10 computed by the present approach are, respectively, 2.0996, 7.4493, 16.6482, 17.0972, 36.3590, 36.3722, 64.1990, 64.2175, 100.1265 which are in close agreement with that reported in [3]. It is also noted that the convergence for even order eigenvalues are, in general, faster (requiring 2 to 3 iterations for reasonable stating guess), while convergence for odd eigenvalues are significantly slower (requiring 7 to 10 iterations instead).
5. CONCLUDING
REMARKS
We have presented an alternative approach to the computation of eigenvalues of Sturm-Liouville problems. Unlike conventional methods such as those using finite element or finite differences schemes, the eigenvalues are computed one at a time with great efficiency and accuracy. If desired, a higher order eigenvalue can also be determined directly (without having to generate all the lower order eigenvalues) if some knowledge of its approximate value is available. A further advantage is that the eigenfunctions are generated simultaneously with the eigenvalues. The versatility of the approach also extends to problems with continuously adjustable coefficients, and parameters. This results in a large class of eigenvalues optimization problem, of which many structural optimization problems are special cases. Lastly, the approach is readily applicable to problems with nonlinear and nonseparated boundary conditions, which could present difficulty for some existing methods.
Sturm-Liouville
Problems
9
APPENDIX A constrained
optimal
parameter selection problem is defined as follows
m;lnso(E) = 90 W),t3
+J1Tg(4QlC,t) 0
2' = fMq,t,t)
>
&
i!E [O,11,
(A-1)
(A.2a) (A.2b)
x(O) = z0(5), subject to the constraints (in canonical form) gi = 0,
i = l,...,Nh,
(A.3a)
gi L 0,
i=N/&+1,...,N,
(A.3b)
where S(5) = Qi(Z(l), r> +
I0
l ri(x(t),
5, t) &
(A.4)
< E RN< is the decision parameter vector to be determined; x(t) E RN= is the state vector of the system (A.2), f, x0, @‘i, Ti, i = 0,l , . . . , N are assumed to be continuously differentiable with respect to each of their respective arguments. It is easy to see that gi, i = 0, 1, . . . , N are functions of c. However, due to the nonexplicit dependence of gi on <, the gradients, and hence, first, order optimality condition cannot be computed directly 8s in the case of nonlinear mathematical programming problem. We shall state without proof the gradient formula as follows; its details can be found in [15]. N assign a costate vector fl E RN=. To compute the gradient, of gi Foreachgi, i=O,l,..., for a given [, solve the differential equation (A.2) forward from t = 0 to t = 1; then solve the corresponding costate differential equation backward from t = 1 to t = 0, namely
(A.5a) (+yl))T
=
a’%(X(l),t) ax(l)
.
(A.5b)
The gradient of gi with respect, to &, the jth component, of <, is given by
(A.61 REFERENCES 1. R.P. Agarwal and Y.M. Chow, Initial-value methods for computing eigenvalues of two-point boundary value problems, Arch. Mathematicum (Bmii) 23, 1-8, (1987). 2. R.S. Anderson and F.R. De Hoog, On the correction of finite difference eigenvalue approximations for SturmLiouville problems with general boundary conditions, BIT 24,401-412, (1984). Collection of finite difference eigenvalues of periodic Sturm-Liouville problems, Journal of 3. A.L.Andrews, Australian Mathematical Society, Series B 30, 46@-469, (1989). 4. M.M. Chawala and C.P. Katti, On Noumerov’s method for computing eigenvalues, BIT 20, 107-109, (1980). 5. M.M. Chawala and C.P. Katti, A new fourth order method for computing eigenvalues of two-point boundary value problems, BIT 20, 511-514, (1980). Boundary-Value Problems, Blaisdell Publ. Comp., Waltham, 6. H.B. Keller, Numerical Methods for Two-Point (1968). Value Problems, SIAM, Philadelphia, PA, (1976). 7. H.B. Keller, Numerical Solution of Two Point Boundary 8. J.W. Paine, F.R. De Hoog and R.S. Anderson, On the correction of finite difference eigenvalue approximations for Sturm-Liouville problems, Computing 26, 123-139, (1981). 9. R.A. Usmani and R.P. Agarwal, New symmetric finite difference methods for computing eigenvalues of a boundary value problem, Communications in Applied Numerical Methods 1,305-309,(1985).
10
C. J. GOH et al.
10. R.A. Usmani and R.P. Agarwal, Some higher order methods for computing eigenvaluee of two-point boundary value problems, Communications in Applied Numerical Methods 3, 5-9, (1987). 11. C.J. Goh and K.L. Teo, MISER: A FORTRAN program for solving optimal control problems, Adv. Engg. Software 10, 90-99, (1988). 12. J.B. Keller, The shape of the strongest column, Archive for Rational Mechanics and Analysis 5, 275-285, (1960). 13. C.J. Goh and K.L. Teo, Control parameterization: A unified approach to optimal control problem, Automatica 24, 3-18, (1988). 14. K.L. Teo and R.S. Womerslay, A control parameterization algorithm for optimal control problems involving linear systems and linear terminal inequality constraints, Numerical Functional Analysis and Optimization 6, 291-313, (1983). 15. K.L. Teo and C.J. Goh, A computational method for combined optimal parameter selection and optimal control problems with general constraints, Journal of Australian Mathematical Society, Series B 30, 350364, (1989). 16. S. Gonzalez and A. Miele, Sequential gradient-restoration algorithm for optimal control problems with general boundary conditions, Journal of Optimization Theory and Applications 26, 395-425, (1978). 17. A. Miele, Recent advances in gradient algorithms for optimal control problems, Journal of Optimization Theory and Applications 17, 361-430, (1975). 18. A. Miele, Gradient algorithms for the optimization of dynamic systems, In Control and Dynamic Systems: Advances in Theory and Applications 16, (Edited by C.T. Leondes), pp. l-52, Academic Press, New York, (1980). 19. A. Miele, J.N. Damoulakis, J.R. Cloutier and J.L. Tietze, Sequential gradient-restoration algorithm for optimal control problems with nondifferentiable constraints, Journal of Optimization Theory and Applications 13, 218-225, (1974). 20. A. Miele, R.E. Pritchard and J.N. Damoulakis, Sequential gradient-restoration algorithm for optimal control problems, Journal of Optimization Theory and Applications 5, 235-282, (197fl).
08957177(94)00065-4
Computing Eigenvalues of Sturm-Liouville Problems via Optimal Control Theory C. J. GOH AND K. L. TEO Department of Mathematics The University of Western Australia, Australia R. P. AGARWAL Department of Mathematics National University of Singapore, Singapore (Received
and accepted
February 1994)
Abstract-In this paper, we shall address three problems arising in the computation of eigenvalues of Sturm-Liouville boundary value problems. We first consider a well-posed Sturm-Liouville problem with discrete and distinct spectrum. For this problem, we shall show that the eigenvalues can be computed by solving for the zeros of the boundary condition at the terminal point as a function of the eigenvalue. In the second problem, we shall consider the case where some coefficients and parameters in the differential equation are continuously adjustable. For this, the eigenvalues can be optimized with respect to these adjustable coefficients and parameters by reformulating the problem as a combined optimal control and optimal parameter selection problem. Subsequently, these optimized eigenvalues can be computed by using an existing optimal control software, MISER. The last problem extends the first to nonstandard boundary conditions such as periodic or interrelated boundary conditions. To illustrate the efficiency and the versatility of the proposed methods, several nontrivial numerical examples are included.
Keywords-Sturm-Liouville problems, Eigenvalues computation, Optimal parameter selection, Combined optimal control and optimal parameter selection, Optimal control software.
1. INTRODUCTION Finding eigenvalues and eigenfunctions of a regular Sturm-Liouville problem is an important question in mathematical physics and engineering. In fact, for the computation of eigenvalues of such problems, several successful algorithms are available, e.g., [l-lo]. In Section 2 of this paper, we shall show that the problem of finding the eigenvalues of regular Sturm-Liouville problems can be readily solved by computing for the zeros of the boundary condition at the terminal point as a function of the eigenvalue. Thus, any numerical scheme for solving the zeros of nonlinear functions can be used for this purpose. It is also known that gradient type methods are more efficient than those which require only the function evaluation, an example of this being the Newton method. Note that though the boundary condition is regarded as a function of the eigenvalue, A, the dependence on X is nonexplicit through the differential equation, and thus, the gradient has to be computed in a slightly round-about way. Nevertheless, such gradient formula can be easily derived, and hence, render the proposed approach of solving for eigenvalues a highly efficient and accurate one. It is also free of truncation error associated with some finite difference and finite element schemes. Typ=et by AM-W 1
2
C. J. GOH et al.
In Section 3, we shall consider the case where some coefficients and parameters in the differential equation are continuously adjustable. In this case, the eigenvalues are clearly dependent on how these coefficients and parameters are adjusted, the scope for optimization arises naturally. Analytical solution to this class of problems does exist occasionally, but more often than not, they admit only numerical solution. To compute for numerical solution, we reformulate the problem as a combined optimal control and optimal parameter selection problem. Subsequently, these optimized eigenvalues can be computed by using an existing optimal control software MISER [ll]. To illustrate the efficiency of the proposed method, we consider the celebrated problem of Keller [12]. In brief, the problem is to find the cross-sectional area of a column of constant length and volume such that the Euler buckling load (i.e., the lowest eigenvalue) is maximized. Finally, in Section 4, we shall extend the class of fixed boundary value problems of Section 2 to problems with periodic, semiperiodic or interrelated boundary conditions. Apart from additional constraints, the eigenvalues are readily computable as in the simpler cases.
2. REGULAR
STURM-LIOUVILLE
PROBLEMS
Consider a regular Sturm-Liouville boundary value problem
(p(t)y’)’+ q(t)y + WGY = 0,
(2.1)
UOY(O> + ~lY'(O> = 0,
u;+u::# 0,
(2.2a)
boy(l) + hY’(l)
b; + bT: # 0,
(2.2b)
= 0,
where q, T E C[O, l], p E C(l)[O, 11, and p(t) > 0, r(t) 2 0 in [O,l]. It is well known that for the problem (2.1), (2.2) there exists a countably infinite number of eigenvalues {X, : n = 1,2, . . . }. These eigenvalues can be arranged as a monotonic increasing sequence Xi < X2 < .+’ such that X, + 00 as n --+ 00. Furthermore, the eigenfunction &(t) (unique up to a multiplicative constant) corresponding to the nth eigenvalue X, has exactly (n- 1) zeros in the open interval (0,l). Setting y = xi and p(t)y’ = x2, the Sturm-Liouville problem can be expressed as (2.3a)
with the initial conditions cy and /3 being constrained by aoa + boxi
alp
+ bl &z(l) P(l)
=
0,
= 0,
u;+uf# 0, b; + b; # 0.
and
(2.4a) (2.4b)
Note that the eigenvalues of the problem (2.3) remain the same for different choices of the constants a and ,B as long as they are chosen such that the constraints (2.4a) are satisfied. Let us consider the case where al = 0 in (2.4a). This implies that (Y is necessarily 0, and p must be necessarily nonzero, since otherwise the problem admits only the trivial solution. Without loss of generality, we may arbitrarily choose p to be any nonzero value, say unity. This will also fix the corresponding eigenfunctions since the arbitrary multiplicative constant is prescribed through fixing p. Similarly, if ue = 0, p is necessarily zero and Q may be normalized to unity by the above argument. If we regard X as a parameter, the eigenvalues of (2.3), (2.4) are obviously the zeros of the terminal boundary condition g(X) = b,,xi(l)
1 + bl-x2(1) P(l)
= 0.
(2.5)
Sturm-LiouvilleProblems
3
It is clear that g(X) is a continuous function of X. This is, in essence, a variation of the shooting methods, e.g., [1,6,7]. It is well known (e.g., [S]) that there are countably infinite number of distinct zeros for the function g(X). To compute these zeros, one may use any of the existing numerical algorithms. It is known that algorithms that require gradient information are in general more efficient than those which require only function evaluation. Nevertheless, since the dependence of g on X is nonexplicit, one has to compute the gradient of g(X) in a somewhat round-about way, e.g., [1,6,7]. As this will come under the context of an optimal parameter selection problem, we shall summarize the specific gradient formula for the problem as follows; and leave the definition of an optimal parameter selection problem and the details of its gradient calculation to the Appendix. For a given X, the gradient of g(X) can be computed as follows: Solve the differential equations (2.3) from t = 0 to t = 1. Then, solve the costate differential equations $4(t) = Q(w2(t)
+
WW2($
(2.6a)
backward from t = 1 to t = 0. Finally, calculate the gradient of g(X) according to the formula g’(X) = 1’
-r(t)a(wz(qdt.
(2.7)
0
With this gradient formula, we can use any gradient type method (such as Newton’s method) to calculate the zeros of the function g(X). We note that this gradient formula is different from the usual one used for the shooting methods [1,6,7]. Alternatively, we may also reformulate the problem of finding zeros of g(X) to one that searches for local minima of (g(X))2. This is, in effect, an optimal parameter selection problem as discussed in the Appendix. To illustrate the efficiency of the proposed approach, six examples are solved using MISER [ll]. This is a FORTRAN program for solving constrained optimal control problem, which includes optimal parameter selection problems as special cases. The computed results for the first five eigenvalues of each problem are tabulated in Table 1. Some remarks should be made about the choice of the initial guess and search range for computing the particular kth eigenvalue &. Apart from being continuous and bounded, g(X) is also approximately periodic in 6. To substantiate this, we take the special case of p(t) = r(t) = 1, q(t) = 0, a0 = bo = 1, ai = br = 0. It can be easily shown that, by normalizing the initial condition to y’(O) = fi, the exact expression of g(X) = y(l) is g(X) = sin 45.
(2.3)
For the general problem, the known asymptotic formula
(2.9) can be used to establish a good initial guess xi as well as search range [&, A,] for the lc-th eigenvahie &, once &_.I is computed. This is done as follows:
Ak-1,
(2.10) (2.11) (2.12)
4
C. 3. GOH et al. Table 1. Six examples for eigenvaluea calculation.
Problem
Order A,
-
Initial
Number of
Guess
Iterations
A’
(I) p(t) = 1
1
40
q(t) = 0
2
181
215.76950 514.35637
r (t) = t(1 -t)
-
45.24288
3
485
a0 = bo = 1
4
914
940.97002
01 = bl = 0
5
1470
1495.59860
(II) PV) = 1
1
30
20.20543
q(t) = 0
2
80
87.85666
3
200
203.23746
ao=b=l
4
300
366.36550
al = bl = 0
5
450
577.21640
1
9
2 3
,lo=bo=l
21 = bl = 0
r(t) = sin
-
t
W) P(t) = 1 q(t) = 0 1*(t)
= cash
t
1 -
2
8.59065
26
2
33.89836
70
2
76.07736
4
150
3
135.10375
5
250
3
211.05166
W) p(t) = 1
1
9
2
7.63809
q(t) = 0
2
25
2
30.14484
r (t) = 1 +
3
65
2
67.60399
zQ=bo=l
4
115
2
120.04479
l1 = bl = 0
5
165
3
187.46795
30.45825
t2
1 1 -
W) 1
20
3
2
100
2
96.03750
(4 = *
3
200
1
199.71217
1;ro=bo=l
4
300
3
343.78217
5
450
3
528.78879
P(t) = 1
1
12.7
2
11.92434
Q(t) = -sint
2
45
2
47.33529
I
r(t) = cost
3
90
2
106.38230
1@l=bi)=l
4
160
3
189.09010
5
260
3
295.41930
P(t) = 1 ) = -(l+
qtt l.
-
1l1
t)4
= bl = 0 WI)
~1 = bl = 0
1
For the first eigenvalue, the available bound
(2.13) can be used to establish a good initial guess and search range for X1. Once this is obtained, the initial guess and search range for higher order eigenvalues can be computed recursively through equations (2.10)-(2.12). Note that though higher order bounds similar to (2.13) exist, they are of little help since they may be overlapping.
Sturm-Liouville
Problems
5
Unlike the traditional methods, the proposed method generates the eigenfunctions (unique up to a multiplicative constant) simultaneously with the eigenvalues through the solutions of the differential equation (2.3). Once the calculation has converged, the order of the eigenvalue can be confirmed easily by checking the number of zeros of the corresponding eigenfunction in the open interval (0,l). In obtaining the results of Table 1, we deliberately allow large variation in the initial guess of Xz (as oppose to what is suggested by (2.8)) to demonstrate the versatility and robustness of the method. In all six cases, no more than two or three iterations are required for convergence in spite of the remote initial guess used.
3. OPTIMIZATION
OF EIGENVALUES
Consider the generalized boundary value problem involving adjustable coefficients and parameters in its differential equation
Mt, 4% 5)Y')' + 4 (t74th 5)Y + XT-(t, u(t) ,5)Y = 0,
(3.la)
ooK)Y(O) + M)Y’(O)
= 0,
(3.lb)
bo(0Y(I)
= 0.
(3Sc)
+ WY’(I)
The functions p, q, T, a~, al, bo and br are assumed to satisfy the usual regularity conditions. A control function u( .) and a parameter < are to be chosen from U x R such that the lowest eigenvalue is maximized, where U denotes the class of all piecewise continuous real-valued functions defined in [0,11. As in the previous section, we rewrite (3.1) in the first order system form
(3.2a)
while the boundary conditions (3.lb),
(3.1~) are constrained by
91 = oo(<)o + Ul(S)P = 0,
(3.3a)
gz = boWi
(3.3b)
+
Note that CYand /I are arbitrary as long as they are chosen such that the boundary conditions (3.3) are satisfied without ending up with the trivial solution. As it was explained before, the nonunique eigenfunction may be fixed by normalizing either CYor p. Now the problem of maximizing the lowest eigenvalue may be posed as the following combined optimal parameter selection and optimal control problems: Subject to the system of differential system (3.2) with the boundary conditions (3.3), find a parameter < and a control u E U such that the cost functional go = --A
(3.4)
is minimized. To prevent X from converging to a higher order values, a suitable upper bound to this lowest eigenvalue should be imposed. Note also that here X is regarded as a decision parameter in addition to c. Sometimes, further restraint on the control function u(e) such as the isoperimetric constraint may often be required 1 u(t) dt = ii = constant. (3.5) I0
C. J. GOH et al.
6
The previous assumption that u(.) is piecewise continuous can be easily extended to functions of arbitrary smoothness. To solve the eigenvalue minimization problem, the optimal control software MISER [ll] is used. The underlying concept used is that of control parameterization [13-151. Essentially, the interval [O,l] is partitioned into a number of equal sub-interval, and a piecewise constant function defined on these sub-intervals is used to approximate the function u(.). As such, the optimal control problem is reduced to an optimal parameter selection problem with the constant function value in these subintervals as decision parameters. Alternatively, the eigenvalue minimization problem can also be solved by the gradient-restoration method [16-201. If r~(.) is required to be arbitrarily smooth, v additional differential equations can be appended to (3.2) via xi =
x4
x3(0)
x& =
x5
x4(0)
. . . x’l+v xi+,
=
. . . X2$-Y
=
v(t)
=
u
= il = . . .
~1+,(0)
i2
(3.6)
= i-l
x2+v(O)
= 6v,
where x2+i
(i-l)
i=1,2
I
)...,
Y.
(3.7)
It is clear that if v(t) is allowed to be piecewise constant, u(e) will be in C(“-‘)[O, 11. The initial conditions ii, i = 1,. . . , v, introduced in (3.6) are also regarded as decision parameters along E, a, P and A. To illustrate the applicability of the proposed approach, we column problem of Keller [12]. In brief, this problem considers volume. The aim is to find the variable cross sectional area (A(t) load, i.e., the lowest eigenvalue, is maximized. Mathematically, that the lowest eigenvalue X of the boundary value problem
(AW2y” + Xy
= 0,
t E
consider the celebrated strongest a column of constant length and say,) such that the Euler buckling the problem is to find A(t), such
[O,11,
Y(0) = Y(l) = 0
(3.8a) (3.8b)
is maximized. In this problem, we have normalized the length of the column to 1, and the isoperimetric constraint on the cross sectional area is given by 1
J
A(t) dt = vor
0
(3.9)
where vs is the volume of the column which may also be normalized to 1. Since A(t) is adjustable, it may be regarded as the control. In [12], Keller has shown that the analytical solution of optimum ;\ and A(t) are 4 A* = 3 sin2 0,
8 E [0, ~1,
t=i(,-isin20),
X’ = $2
2 13.1595.
(3.10) (3.11) (3.12)
Using the proposed approach, the strongest column problem of Keller may be restated as (3.13)
Sturm-Liouville
Problems
7
subject to the system of differential equations
x; = x2,
Xl(O)
x x; = -_-2x1,
=
(3.14a)
0,
Q(O) = 1 (normalized),
(3.14b)
x3(0)
(3.14c)
x3 xg = v,
= <,
and the constraints (3.15a)
0 5 A I 2& 91 =
Xl(l)
(3.15b)
= 0,
1 Q2 =
I0
x3(t)
dt - 1 = 0.
(3.15c)
Note that by choosing piecewise constant function for v, the cross-sectional area A(t) = x3(t) is, therefore, piecewise linear. To solve this combined optimal control and optimal parameter selection problem by MISER, we need to supply the number Np of partitions of the interval [O,l], and the initial guess for the control v and the parameters X and E. Using Np = 20, co = 1, v”(t) = 0 (and, hence, x!(t) = AO(t) = 1) and X0 = 10, the solution converges after 25 iterations to X* = 13.1551 (cf. exact value is 13.1595). The corresponding A*(t)
was found to be almost identical to the exact solution of Keller given
by (3.10), (3.11). As a straight forward extension to the above, we modify the boundary condition for (3.8b) to y’(0) = Y’(l) = 0.
(3.16)
It is easy to show that the solution for this modified problem yields the same X* as before, as confirmed by numerical calculation. However, the optimum A*(t) differs from the previous case (i.e., A*(t)) by a shifted symmetry
a*(t)
=
The above boundary conditions correspond to the hinged-hinged and clamped-clamped cases. Other variations of boundary condition such as hinged-free, hinged-clamped, clamped-free, etc., may also be easily computed. Interior point constraints, such as an internal hinge, can also be incorporated with slight modification.
4. NONSEPERATED
BOUNDARY
CONDITIONS
In the literature, boundary conditions of the Consider the differential equation (2.1). type (2.2a)-(2.2b) are referred to as separated boundary conditions. Problems with linear nonseparated boundary conditions of the type aov(O)
+ NY’(o)
= boy(l)
are slightly more difficult to solve in practice. conditions Y(0) = Y(I)
or
+ hY’(l)
(4.1)
This includes the class of periodic boundary Y’(O) = y’(l),
(4.2)
C. J. GOH et al.
8
or semi-periodic
boundary conditions
Y(O)= -Y(l)
or
y’(0) = -y’(l).
(4.3)
Correction techniques associated with both finite element and finite difference methods are available to improve on the accuracy of the computed eigenvalues [2,6-81. In this section, we consider a class of even more general nonseparated boundary conditions of the type fl2i (Y(O), Y’(O), Y(l)? Y’(i)) = 0,
i = 1,2,
(4.4)
where 0i are, in general, nonlinear functions of the arguments. Such nonlinear, nonseparated boundary conditions will unlikely be accessible to many of the existing methods in the literature. Nevertheless, it presents no added difficulty to our present formulation. As before, we assign the initial condition y(O) and g’(O) with suitable decision parameters cx and p, respectively. The nonseparated boundary condition (4.4) will thus reduce to the terminal state constraint G
Ccr, 07 51(l)r
22(1))
=
07
i =
1,2.
(4.5)
By arbitrarily fixing one of CYand p such that the resulting eigenfunction is nontrivial, (4.5) effectively reduces to two simultaneous nonlinear equation in two unknowns X and cy (or ,Q, Gradient formulae can be easily derived using the result of the Appendix, and hence, Newton’s method can again be used to solve for the zeros. The similar argument will still be valid if the differential equation involves continuously adjustable coefficients and parameters as in (3.la). For illustration, an example taken from [16] is computed using the present approach y” - 107r2 cos(27rt)y + 7r2xy = 0,
(4.6)
with the periodic boundary conditions Y(O)
= Y(l),
(4.7a)
Y’(O)
= Y’(l).
(4.7b)
Here, p(t) = 1, q(t) = -10n2cos(27rt), r(t) = x2. The eigenvalues &, k = 2,3,. . . , 10 computed by the present approach are, respectively, 2.0996, 7.4493, 16.6482, 17.0972, 36.3590, 36.3722, 64.1990, 64.2175, 100.1265 which are in close agreement with that reported in [3]. It is also noted that the convergence for even order eigenvalues are, in general, faster (requiring 2 to 3 iterations for reasonable stating guess), while convergence for odd eigenvalues are significantly slower (requiring 7 to 10 iterations instead).
5. CONCLUDING
REMARKS
We have presented an alternative approach to the computation of eigenvalues of Sturm-Liouville problems. Unlike conventional methods such as those using finite element or finite differences schemes, the eigenvalues are computed one at a time with great efficiency and accuracy. If desired, a higher order eigenvalue can also be determined directly (without having to generate all the lower order eigenvalues) if some knowledge of its approximate value is available. A further advantage is that the eigenfunctions are generated simultaneously with the eigenvalues. The versatility of the approach also extends to problems with continuously adjustable coefficients, and parameters. This results in a large class of eigenvalues optimization problem, of which many structural optimization problems are special cases. Lastly, the approach is readily applicable to problems with nonlinear and nonseparated boundary conditions, which could present difficulty for some existing methods.
Sturm-Liouville
Problems
9
APPENDIX A constrained
optimal
parameter selection problem is defined as follows
m;lnso(E) = 90 W),t3
+J1Tg(4QlC,t) 0
2' = fMq,t,t)
>
&
i!E [O,11,
(A-1)
(A.2a) (A.2b)
x(O) = z0(5), subject to the constraints (in canonical form) gi = 0,
i = l,...,Nh,
(A.3a)
gi L 0,
i=N/&+1,...,N,
(A.3b)
where S(5) = Qi(Z(l), r> +
I0
l ri(x(t),
5, t) &
(A.4)
< E RN< is the decision parameter vector to be determined; x(t) E RN= is the state vector of the system (A.2), f, x0, @‘i, Ti, i = 0,l , . . . , N are assumed to be continuously differentiable with respect to each of their respective arguments. It is easy to see that gi, i = 0, 1, . . . , N are functions of c. However, due to the nonexplicit dependence of gi on <, the gradients, and hence, first, order optimality condition cannot be computed directly 8s in the case of nonlinear mathematical programming problem. We shall state without proof the gradient formula as follows; its details can be found in [15]. N assign a costate vector fl E RN=. To compute the gradient, of gi Foreachgi, i=O,l,..., for a given [, solve the differential equation (A.2) forward from t = 0 to t = 1; then solve the corresponding costate differential equation backward from t = 1 to t = 0, namely
(A.5a) (+yl))T
=
a’%(X(l),t) ax(l)
.
(A.5b)
The gradient of gi with respect, to &, the jth component, of <, is given by
(A.61 REFERENCES 1. R.P. Agarwal and Y.M. Chow, Initial-value methods for computing eigenvalues of two-point boundary value problems, Arch. Mathematicum (Bmii) 23, 1-8, (1987). 2. R.S. Anderson and F.R. De Hoog, On the correction of finite difference eigenvalue approximations for SturmLiouville problems with general boundary conditions, BIT 24,401-412, (1984). Collection of finite difference eigenvalues of periodic Sturm-Liouville problems, Journal of 3. A.L.Andrews, Australian Mathematical Society, Series B 30, 46@-469, (1989). 4. M.M. Chawala and C.P. Katti, On Noumerov’s method for computing eigenvalues, BIT 20, 107-109, (1980). 5. M.M. Chawala and C.P. Katti, A new fourth order method for computing eigenvalues of two-point boundary value problems, BIT 20, 511-514, (1980). Boundary-Value Problems, Blaisdell Publ. Comp., Waltham, 6. H.B. Keller, Numerical Methods for Two-Point (1968). Value Problems, SIAM, Philadelphia, PA, (1976). 7. H.B. Keller, Numerical Solution of Two Point Boundary 8. J.W. Paine, F.R. De Hoog and R.S. Anderson, On the correction of finite difference eigenvalue approximations for Sturm-Liouville problems, Computing 26, 123-139, (1981). 9. R.A. Usmani and R.P. Agarwal, New symmetric finite difference methods for computing eigenvalues of a boundary value problem, Communications in Applied Numerical Methods 1,305-309,(1985).
10
C. J. GOH et al.
10. R.A. Usmani and R.P. Agarwal, Some higher order methods for computing eigenvaluee of two-point boundary value problems, Communications in Applied Numerical Methods 3, 5-9, (1987). 11. C.J. Goh and K.L. Teo, MISER: A FORTRAN program for solving optimal control problems, Adv. Engg. Software 10, 90-99, (1988). 12. J.B. Keller, The shape of the strongest column, Archive for Rational Mechanics and Analysis 5, 275-285, (1960). 13. C.J. Goh and K.L. Teo, Control parameterization: A unified approach to optimal control problem, Automatica 24, 3-18, (1988). 14. K.L. Teo and R.S. Womerslay, A control parameterization algorithm for optimal control problems involving linear systems and linear terminal inequality constraints, Numerical Functional Analysis and Optimization 6, 291-313, (1983). 15. K.L. Teo and C.J. Goh, A computational method for combined optimal parameter selection and optimal control problems with general constraints, Journal of Australian Mathematical Society, Series B 30, 350364, (1989). 16. S. Gonzalez and A. Miele, Sequential gradient-restoration algorithm for optimal control problems with general boundary conditions, Journal of Optimization Theory and Applications 26, 395-425, (1978). 17. A. Miele, Recent advances in gradient algorithms for optimal control problems, Journal of Optimization Theory and Applications 17, 361-430, (1975). 18. A. Miele, Gradient algorithms for the optimization of dynamic systems, In Control and Dynamic Systems: Advances in Theory and Applications 16, (Edited by C.T. Leondes), pp. l-52, Academic Press, New York, (1980). 19. A. Miele, J.N. Damoulakis, J.R. Cloutier and J.L. Tietze, Sequential gradient-restoration algorithm for optimal control problems with nondifferentiable constraints, Journal of Optimization Theory and Applications 13, 218-225, (1974). 20. A. Miele, R.E. Pritchard and J.N. Damoulakis, Sequential gradient-restoration algorithm for optimal control problems, Journal of Optimization Theory and Applications 5, 235-282, (197fl).