- Email: [email protected]
Numerical experiments on the Stodola-Vianello method
Journal of Sound and Vibration (1995) 183(5), 925–927
LETTER TO THE EDITOR NUMERICAL EXPERIMENTS ON THE STODOLA–VIANELLO METHOD P. A. A. L Institute of Applied Mechanics (CONICET) and Department of Engineering, Universidad Nacional del Sur, 8000 Bahı´ a Blanca, Argentina L. E Mechanical Systems Analysis Group, Facultad Regional Bahı´ a Blanca (U.T.N.) , Argentina G. S´ S Department of Physics, School of Engineering, Universidad Nacional de Buenos Aires, Argentina (Received 6 April 1994, and in final form 26 September 1994)
1. As stated by Hildebrand in his classical textbook [1], ‘‘the method of Stodola and Vianello is a useful iterative procedure which allows for the approximate determination of the characteristic numbers and functions of a boundary value problem’’. Hildebrand illustrated the procedure by considering the problem of a vibrating string fixed at its ends and governed by the differential system d2y/dx 2 = −ly,
y(0) = y(L) = 0,
(1a–c)
where l is the eigenvalue under investigation. Replacing y on the right side of equation (1a) by a first approximation y1 (x) which satisfies equations (1b, c), one is able to integrate equation (1a) and obtain a second approximation which can be written in the form y = lf1 (x).
(2)
Requiring now that y1 (x) and lf1 (x) agree as well as possible (in some sense), for instance integrating them over the interval (0, L) one is able to obtain an approximate value of l. One may improve further the result by repeating the procedure: e.g., substituting f1 (x) in equation (1a) and obtaining y = lf2 (x). As shown by Hildebrand convergence of the procedure is achieved in some instances and upper bounds are obtained when polynomials with integer powers of the variable are used. The present note describes numerical experiments in which the approximate co-ordinate functions are simple polynomials which contain a minimization exponential parameter g. By judiciously selecting g one is able to improve the calculated eigenvalue determined in the first cycle, although one cannot guarantee the existence of a lower or upper bound since non-integer powers of the independent variable are used. Three numerical examples are presented in the following section. 2. In the case of the differential system (1) one takes y1 (x) = x − L(x/L)g.
(3)
925 0022–460X/95/250925 + 03 $08.00/0
7 1995 Academic Press Limited
926
Substituting this into the right side of equation (1a), integrating twice and using the boundary conditions (1b, c) one obtains
$
y(x) = l
1
L
g−1
0
1%
xg + 2 x3 1 1 − + L2 − x . (g + 1)(g + 2) 6 6 (g + 1)(g + 2)
(4)
Integrating equations (3) and (4) over the interval (0, L) and equating the results, one obtains l1(1)L 2 =
$
1 1 − 2 g+1
%>$
%
1 1 1 I − + = 1. (g + 1)(g + 2)(g + 3) 2(g + 1)(g + 2) 24 I2
(5)
A second cycle of the interaction is performed by substituting equation (4) in equation (1a). An improved value of y(x) is obtained, and repeating the previous procedure one obtains g1(2)L 2 = I3 /I4 ,
(6)
where
0 1
1 1 1 1 1 I3 = − + − , S 24 6 R 2 R = (g + 1)(g + 2),
I4 = −
S = R(g + 3),
$
0 1
%
1 1 1 1 1 1 − W , − + − 2 U 720 24 6 R T = S(g + 4),
U = T(g + 5),
1 1 1 1 W= − + − . T 6R 36 120 Minimizing expression (5) with respect to g, one obtains l1(1)L 2 = 9·983, while the same procedure applied to expression (6) yields g2(2)L 2 = 9·88. The exact eigenvalue is p 2 2 9·87, and the approximate values obtained in reference [1] are g1(1)L 2 = 10 and g1(2)L 2 = 9·882 (which correspond to taking a fixed value of g = 2). It is observed that a slight improvement has been achieved. Next, consider the buckling of a simply supported beam, the flexural rigidity of which varies according to the functional relation Cx/L. The differential system is defined by [1] d2y/dx 2 = −(m 2/L)(y/x),
y(0) = y(L) = 0,
(7a–c)
where m 2 = Pcr L 2/C. Taking y1 (x) = x − L(x/L)g
(8)
and following the previously explained procedure one obtains, as a first approximation, m1(1) =
6$
1 1 − 2 g+1
%>$
1 1 1 + − g(g + 1)(g + 2) 2g(g + 1) 12
%7
1/2
,
(9)
This functional relation is a monotonically increasing function of g. Obviously, the value g = 1 is not valid, since equation (8) will be identically zero. On the other hand, values of g very close to unity, say g = 1·01, will correspond to extremely flat deflection shapes. It seems reasonable to choose a value of g which differs considerably from g = 2 (which yields m1(1) = 2 [1]) but which still maintains a reasonable shape for the deflected beam configuration. For instance, for g = 1·10 one obtains l1(1) = 1·90. On the other hand, for g = 1·01 a value of g1(1) = 1·898 results, and for g = 1·40 one obtains g1(1) = 1·94. Accordingly, a pragmatic selection of g (for instance, g = 1·10) seems to be a reasonable one, but using other values will not introduce considerable differences, the exact value being g1 = 1·916 [1].
927
Now consider the case of axisymmetric vibrations of a circular membrane of unit radius. The differential system is now (1/r)(d/dr)(r du/dr) = −l 2u,
u(1) = 0,
(10a, b)
taking u1 = 1 − r g
(11)
and following the usual procedure one obtains l1(1) =
6$
1 1 − 2 g+2
%>$
1 1 1 + − (g + 2)2(g + 4) 2(g + 2)2 16
%7
1/2
.
(12)
With g = 2 one obtains l1(1) = 2·45, while for g = 1·01 the value l1(1) = 2·39 is determined. Following the reasoning previously stated, one selects g = 1 = 1·10 and determines l1(1) = 2·40, which is in excellent agreement with the exact eigenvalue: l1 = 2·4048. 3. Summarizing, one can conclude the following. The classical Stodola–Vianello iterative procedure converges in many situations of mathematical physics of practical interest, and upper bounds of the exact eigenvalues are determined when co-ordinate functions with integer powers of the independent variable are used, as proved by Hildebrand in his well known treatise [1]. A large number of iterations is possible nowadays by using a computer program, such as MACSYMA, to avoid the tedious algebra, as the complexities of the function increase from one iteration to the next. The minimization procedure proposed in this note leads, obviously, to the use of polynomials with non-integer powers and does not yield upper bounds. However, from a practical viewpoint, and as shown in the three examples presented here, considerable improvement is achieved in the determined eigenvalue when only the first cycle of the Stodola–Vianello method is employed. It does certainly require judicious selection of the non-integer optimization parameter. The technique seems convenient when used in connection with more complex eigenvalue problems. The present study has been sponsored by Secretarı´ a de Ciencia y Technologı´ a of Universidad Nacional del Sur (Project 1994–1995; Program Director Professor R. E. Rossi). The authors are indebted to Professor P. E. Doak and to the referee of the present note for their useful comments. 1. F. B. H 1962 Advanced Calculus for Applications. Englewood Cliffs, New Jersey, Prentice-Hall. See pp. 200–206.
LETTER TO THE EDITOR NUMERICAL EXPERIMENTS ON THE STODOLA–VIANELLO METHOD P. A. A. L Institute of Applied Mechanics (CONICET) and Department of Engineering, Universidad Nacional del Sur, 8000 Bahı´ a Blanca, Argentina L. E Mechanical Systems Analysis Group, Facultad Regional Bahı´ a Blanca (U.T.N.) , Argentina G. S´ S Department of Physics, School of Engineering, Universidad Nacional de Buenos Aires, Argentina (Received 6 April 1994, and in final form 26 September 1994)
1. As stated by Hildebrand in his classical textbook [1], ‘‘the method of Stodola and Vianello is a useful iterative procedure which allows for the approximate determination of the characteristic numbers and functions of a boundary value problem’’. Hildebrand illustrated the procedure by considering the problem of a vibrating string fixed at its ends and governed by the differential system d2y/dx 2 = −ly,
y(0) = y(L) = 0,
(1a–c)
where l is the eigenvalue under investigation. Replacing y on the right side of equation (1a) by a first approximation y1 (x) which satisfies equations (1b, c), one is able to integrate equation (1a) and obtain a second approximation which can be written in the form y = lf1 (x).
(2)
Requiring now that y1 (x) and lf1 (x) agree as well as possible (in some sense), for instance integrating them over the interval (0, L) one is able to obtain an approximate value of l. One may improve further the result by repeating the procedure: e.g., substituting f1 (x) in equation (1a) and obtaining y = lf2 (x). As shown by Hildebrand convergence of the procedure is achieved in some instances and upper bounds are obtained when polynomials with integer powers of the variable are used. The present note describes numerical experiments in which the approximate co-ordinate functions are simple polynomials which contain a minimization exponential parameter g. By judiciously selecting g one is able to improve the calculated eigenvalue determined in the first cycle, although one cannot guarantee the existence of a lower or upper bound since non-integer powers of the independent variable are used. Three numerical examples are presented in the following section. 2. In the case of the differential system (1) one takes y1 (x) = x − L(x/L)g.
(3)
925 0022–460X/95/250925 + 03 $08.00/0
7 1995 Academic Press Limited
926
Substituting this into the right side of equation (1a), integrating twice and using the boundary conditions (1b, c) one obtains
$
y(x) = l
1
L
g−1
0
1%
xg + 2 x3 1 1 − + L2 − x . (g + 1)(g + 2) 6 6 (g + 1)(g + 2)
(4)
Integrating equations (3) and (4) over the interval (0, L) and equating the results, one obtains l1(1)L 2 =
$
1 1 − 2 g+1
%>$
%
1 1 1 I − + = 1. (g + 1)(g + 2)(g + 3) 2(g + 1)(g + 2) 24 I2
(5)
A second cycle of the interaction is performed by substituting equation (4) in equation (1a). An improved value of y(x) is obtained, and repeating the previous procedure one obtains g1(2)L 2 = I3 /I4 ,
(6)
where
0 1
1 1 1 1 1 I3 = − + − , S 24 6 R 2 R = (g + 1)(g + 2),
I4 = −
S = R(g + 3),
$
0 1
%
1 1 1 1 1 1 − W , − + − 2 U 720 24 6 R T = S(g + 4),
U = T(g + 5),
1 1 1 1 W= − + − . T 6R 36 120 Minimizing expression (5) with respect to g, one obtains l1(1)L 2 = 9·983, while the same procedure applied to expression (6) yields g2(2)L 2 = 9·88. The exact eigenvalue is p 2 2 9·87, and the approximate values obtained in reference [1] are g1(1)L 2 = 10 and g1(2)L 2 = 9·882 (which correspond to taking a fixed value of g = 2). It is observed that a slight improvement has been achieved. Next, consider the buckling of a simply supported beam, the flexural rigidity of which varies according to the functional relation Cx/L. The differential system is defined by [1] d2y/dx 2 = −(m 2/L)(y/x),
y(0) = y(L) = 0,
(7a–c)
where m 2 = Pcr L 2/C. Taking y1 (x) = x − L(x/L)g
(8)
and following the previously explained procedure one obtains, as a first approximation, m1(1) =
6$
1 1 − 2 g+1
%>$
1 1 1 + − g(g + 1)(g + 2) 2g(g + 1) 12
%7
1/2
,
(9)
This functional relation is a monotonically increasing function of g. Obviously, the value g = 1 is not valid, since equation (8) will be identically zero. On the other hand, values of g very close to unity, say g = 1·01, will correspond to extremely flat deflection shapes. It seems reasonable to choose a value of g which differs considerably from g = 2 (which yields m1(1) = 2 [1]) but which still maintains a reasonable shape for the deflected beam configuration. For instance, for g = 1·10 one obtains l1(1) = 1·90. On the other hand, for g = 1·01 a value of g1(1) = 1·898 results, and for g = 1·40 one obtains g1(1) = 1·94. Accordingly, a pragmatic selection of g (for instance, g = 1·10) seems to be a reasonable one, but using other values will not introduce considerable differences, the exact value being g1 = 1·916 [1].
927
Now consider the case of axisymmetric vibrations of a circular membrane of unit radius. The differential system is now (1/r)(d/dr)(r du/dr) = −l 2u,
u(1) = 0,
(10a, b)
taking u1 = 1 − r g
(11)
and following the usual procedure one obtains l1(1) =
6$
1 1 − 2 g+2
%>$
1 1 1 + − (g + 2)2(g + 4) 2(g + 2)2 16
%7
1/2
.
(12)
With g = 2 one obtains l1(1) = 2·45, while for g = 1·01 the value l1(1) = 2·39 is determined. Following the reasoning previously stated, one selects g = 1 = 1·10 and determines l1(1) = 2·40, which is in excellent agreement with the exact eigenvalue: l1 = 2·4048. 3. Summarizing, one can conclude the following. The classical Stodola–Vianello iterative procedure converges in many situations of mathematical physics of practical interest, and upper bounds of the exact eigenvalues are determined when co-ordinate functions with integer powers of the independent variable are used, as proved by Hildebrand in his well known treatise [1]. A large number of iterations is possible nowadays by using a computer program, such as MACSYMA, to avoid the tedious algebra, as the complexities of the function increase from one iteration to the next. The minimization procedure proposed in this note leads, obviously, to the use of polynomials with non-integer powers and does not yield upper bounds. However, from a practical viewpoint, and as shown in the three examples presented here, considerable improvement is achieved in the determined eigenvalue when only the first cycle of the Stodola–Vianello method is employed. It does certainly require judicious selection of the non-integer optimization parameter. The technique seems convenient when used in connection with more complex eigenvalue problems. The present study has been sponsored by Secretarı´ a de Ciencia y Technologı´ a of Universidad Nacional del Sur (Project 1994–1995; Program Director Professor R. E. Rossi). The authors are indebted to Professor P. E. Doak and to the referee of the present note for their useful comments. 1. F. B. H 1962 Advanced Calculus for Applications. Englewood Cliffs, New Jersey, Prentice-Hall. See pp. 200–206.