- Email: [email protected]
Fast galerkin method and its application to determine periodic solutions of non-linear oscillators
IN. 1. Non-Luwar Mechanics. Vol. 22. No. 2. pp. 89-98. 1987 Printed in Grwt Britain
0
co2a-7462187 93.00 + .oo 1987 Pcrgamon Journals Ltd.
FAST GALERKIN METHOD AND ITS APPLICATION TO DETERMINE PERIODIC SOLUTIONS OF NON-LINEAR OSCILLATORS F. H. Department
of Engineering
LING
Mechanics,
(Received 9 August
and X. X. WV
Shanghai Jiao Tong University, P.R. China
1985; received for publication
23 September
Shanghai
200030,
1986)
Abstract-The Galerkin method is an approximate method which finds wide application in solving differential and integral equations. But a large amount of computation is needed in order to get a high order approximation by using the method. Applying the FFT technique to form a so-called fast Galerkin method, we can reduce the computation work greatly, when taking trigonometric functions as characteristic functions. Taking the periodic solution of non-linear oscillators as an example, we illustrate the procedure and the efficiency of the method. Moreover, with some modifications we extend the applicability of the method, so that not only periodic solutions with known periods, but also those with unknown periods, as well as subharmonics, combination tones, etc., can be treated with the method. Some techniques are described which can be used to simplify the computation. I. INTRODUCTION
The Galerkin method is an approximate method which is often applied to solve differential equations and integral equations numerically. The basic idea is to seek a solution in terms of a series of characteristic functions which satisfy the boundary conditions. For partial differential equations this procedure leads to a system of ordinary differential equations, which can be treated with, say, the Galerkin method once again. As for ordinary differential equations or integral equations, it leads to a system of algebraic equations which the coefficients of the series must satisfy. Solving these equations, usually by a numerical method, we get a Galerkin approximate solution of the equations. In this paper, we regard the periodic solution of non-linear oscillators as a boundary value problem of ordinary differential equations (as in [l-3]), but the method described here can also be applied to other boundary value problems as well as integral equations. The problem of determining a periodic solution can be put forward as follows, Ti = f(x, t),
XER',
f: UP+1- R,
f(x, t + 2x) = f(x),
x(0) - x(279 = 0,
(1)
where f(x, t) is periodic in t with a period of 271without losing generality. Choosing the following trigonometric functions: {1, sin nt, cos nt>, n = 1,2,. . . , as characteristic functions, we obtain a Galerkin approximation of order m:
x,(0 = a0 + f (a2”_ 1sin nt n=1
+ a2,cos nt),
(2)
where x,(t), ao, a+ I, azn (n = 1, 2,. . . , m) are all I-dimensional vectors. Urabe [4,5] has investigated the above problem. His main results are as follows. Assume that f(x, t) and its derivatives with respect to the x-coordinate are continuously differentiable with respect to the same x-coordinate and t in a closed bounded region of the x-space and the axis of t. If there is an isolated periodic solution x = n(t), there exists a Galerkin approximation x = x,,,(t), which converges uniformly as m + co to the exact solution x = n(t). On the other hand, if the Galerkin approximation x = x,(t) satisfies certain conditions, there must be an exact isolated periodic solution lying in the neighborhood of xl&). Urabe applied the classic Fourier transform (the FFT was not discovered at that time), which requires a large amount of computation, especially for obtaining a higher order 89
F. H. LINCand X. X. Wu
90
approximation, so as to restrict its usage greatly. In this paper we introduce the FFT to form the so-called fast Galerkin method, in order to reduce the amount of computation and hence to raise its applicability. Besides, with some modifications the method can be used to solve not only the periodic solutions with known periods but also those with unknown periods, as well as to get subharmonics and combination tones, etc. We illustrate the efficiency in some numerical examples. Some techniques which can be used to reduce the amount of computation are also listed here. 2. PROCEDURE
FOR
DETERMINING PERIODIC GALERKIN METHOD
SOLUTIONS
BY USING
THE
Substituting the expression (2) into equation (1) and equating the coefficients of sin nt and cosnt on both sides of the equation, we have
FO(4
fcxm(s),
= &
slds = 0,
s 2n
F2,-l(4
= i
Rx,,,(s), s]sin
ns ds + na2,, = 0,
0
fCx,(s),s]cosns
F2,(a) = a
ds - nazn_ i = 0,
(n = 1,2,...,m)
(3)
where a = [a$, a:, . . . , ai,,J’ is a (2m + 1) x I-dimensional vector, which gathers all unknown coefficients. The equations (3) are non-linear algebraic equations for determining a and can be solved with, say, a Newton-Raphson iterative method. In this case a system of linear equations for the correctors Aa of the unknown coefficients results: JAa + F = 0,
(4)
where F = [Fc, FT,. . . , FT,,,]’ is a (2m + 1) x I-dimensional vector, and J = dF/t?a is the first order derivative matrix, whose elements are definite integrals of the products of the elements of the derivative $ = af/a x with cosnt and/or sinnt and can be regarded as the Fourier coefficients of the elements of +. Some other methods, for example, the Broyden method [6], can also be used to solve equations (3); the corresponding iterative scheme takes the form of p(il = =
p(i),
yW =
F(i+
A,(i)
H(‘+
_ H(i)F(i)
1) = H”’
,
1) _
F(i),
_
(i) ,, (i) _
(H
(il
(i)rH(i)/p(i)TH(i)y(i).
P )P
(5)
Thereby, we sum up the procedure of using the Galerkin method to determine a periodic solution as follows: (1) Estimate the initial value a(‘) of the coefficient vector a in expression (2); (2) Select 2m, equidistant nodes ti, i = 1, 2,. . . , 2me, in the time interval (0,2n), then repeat the following steps with k = 0, l,-2,. . .: (a) Evaluate x$‘(ti) from a (kl by using the inverse Fourier transform; (b) Calculate fCX~‘(fi),ti], and when using the Newton-Raphson method, also calculate
+lIx!tt’(ti)9 lil; (c)Calculate J
from $[X$‘(tJ, ti] with the Fourier transform, when using the NewtonRaphson method; (d) Evaluate F from fCX~‘(Ci),ti] with the Fourier transform according to equations (3); (e) Solve equation (4) when using the Newton-Raphson method, or equation (5) when
Fast Galerkin method
91
using the Broyden method to determine correctors Aa”‘; (f) End the iteration, if r is small enough, otherwise set z(~+I) = zCk)-I- A,(‘) and return to step (a). 3. FORMING
THE FAST GALERKIN METHOD TECHNIQUE
BY INTRODUCING
THE FFT
From the procedure stated above, we see that the computational work in the Galerkin method comes from three source% the solution of the algebraic equations [step2(e)], the calculation of f and $ at nodes ti [step2(b)], and the Fourier transforms [steps 2(a), (c) and (d)]. For the first one, there are many effective methods, such as the Broyden method and the super relaxation iterative method, which both require 0(4m*I*) times of multiplications, where m is the order of the harmonics to be found. For the second one, generally snaking, the number of multiplications is of O(Zm,& which varies with the complexity of the equation, where m, is the order of the Fourier transform. We must have m, > 2m for a reason to be stated later. For the third one, the classical discrete Fourier transform requires a large amount of computation. In the steps 2(a), (c)and (d), the numbers of multiplications are 1(2m,)2, f2(2m,)2 and I(2mJ2 respectively, and their sum is (I f 2)t(2mJ2. Introducing the FFT technique (see [?‘J, for example) into the Galerkin method may reduce the number of multiplications to l/Z(I + 2~I(2m~)log~(2m~) with the FFT algo~thm of index 2, i.e. only 18% of the original computer work for m, = 4, 12.5% for m, = 8, and 7.8% for m, = 16. Now we explain the reason of taking m, 3 2m,where m is the order of the harmonics to be found and m, is the order of the Fourier transform. As is well known, the coefficients of higher harmonics evaluated with the Fourier transform are no longer valid, and usually only the first half of the m, coefficients can be used. On the other hand, in order to reach a high degree of accuracy, the value of r, which is equal or greater than )/P(t) - f@(t), t)ll, must be kept small enough, where S(t) is an approximate solution. Expanding (k(t) - f@(t), t)) into a Fourier series, we have a(t) - f@(r), t) = A,, f
i
(Azn_ I sin nt + Al,, cosat),
&?=1
which is equal approximately
to
k(t) - f@(t), t) = A0 + 3 (Azn-. 1sin nt + Aln cos d. n-1
From equation (3), we find that (i = O,l,... ,2m).
Ai =z: Fi, Let
EB = m,ax A,, + 2 (A,,_ I sin nt + Az,cosnt) II
EE = rn,“x f$+ li
,
n-1 1(AZ,--I sin nt -I- A2*cos nt) , 1/
then we can take r = EB + EE. The practical calculation shows that the control of EE is very important. See next section for details. This requires that m, 3 2m, otherwise EE cannot be determined. In this way we see that the Fourier transform is a time consuming work in the method and the application of the FFT technique reduces the amount of computation greatly. Taking the method with the Broyden method as an example, in a two-degree-of-freedom
F. H, LING and X. X. WV
92
system (I = 4), we find that to determine the harmonics through the 8th order (m = 8, taking m, = 16), the above three kinds of computation need 4 x 8’ x 42 = 4096,2 x 16 x 4 = 128, and 2 x 4 x 322 = 8192 times of multiplications respectively, and altogether 12,416 times. After adopting the FFT, the third one decreases to 2 x 4 x 32 x log,32 = 1280 times, and the sum reduces to 5504 times, which is only 44.3% of the original. In compiling the program we take m, = 2m, and EE and EB are enlarged to
J
A; + f (At,It=1
EB =
EE =
J
z
n=m+ 1
(A;,-,
1 +
Ain)
+ A:,)
respectively. 4. SOME
REMARKS
4.1. Choice of the initial value a(” of the coejkient vector a This is an important problem which has a strong effect on the convergence speed and even the convergence itself. The commonly used methods are: (a) Make an extrapolation of the solutions obtained to the equations with slightly different coefficients to give the initial value of the solution which is sought. For example, we can make use of the solution with neighbouring forcing frequencies to determine the response curves. For weak non-linear systems we can make use of the solution of a corresponding linearized system. (b) Use a lower-order Galerkin approximation of the same system. A lower-order Galerkin approximation is denoted by x = a,* + z (a&_, sinnt + ar,cos nt). PI=1
Thus the initial value may be chosen as (0) =
%
a? 1
a!O) I = 0.
I
(i = 0,1,...,2m,) (i = 2m, + 1,. . . ,2m,).
4.2. Necessity of controlling the residue EE The well-known harmonic balance method may be regarded as a low-order Galerkin method, in which only the conditions Fi = 0 (i = 1,. . . , m), i.e. EB = 0 are demanded. But according to our experience gathered in the computation we establish that this is not sufficient. It may happen that an approximate solution with a small EB but a large EE is far away from the exact solution. So we must check the condition of EE < EPS after the condition of EB < EPS is satisfied. If the condition EE < EPS is not satisfied, we should increase the order m, and if the error EE is still unacceptable even with a large m, then there seems to be no proper periodic solution nearby. 4.3. A comparison between the Newton-Raphson here
method and the Broyden method used
Both methods are tested in -the computation. With the Newton-Raphson method 0(9m313)times of multiplications are required for solving the linear algebraic equation, and the derivative matrix J has to be evaluated. So the program is complicated, and the numerical stability in the computation is poor, i.e. it requires a fairly good initial guess, otherwise it would diverge. The advantage of the quick convergence of the NewtonRaphson method cannot compensate for these disadvantages. With the Broyden method the computation of J is avoidable, and we obtain a concise program. Although it converges more slowly (in about 20 iterations), the total amount of computation is less than that
93
Fast Galerkin method
with the Newton-Raphson method, because it needs only 0(4m21*) times of multiplications for every iteration. Moreover, the numerical stability of the Broyden method is better than that of Newton-Raphson method, i.e. the convergence of the method is nearly always guaranteed even with a poor initial guess. 4.4. Treatment of a second order equation Most non-linear oscillation equations are of second order and need not be transformed into the standard form (1). We substitute the expression (2) directly into the equations, by noticing that k,(t) = “cl n(a2”- l cos nt - a2” sin nt)
In this way, the number of unknown coefficients is only half of that for the transformed first-order equations, so the amount of computation is reduced to nearly one-fourth. This idea can be extended to a more general case, in which the equations are of higher order. We even transform a set of first-order equations into a set of higher-order equations in order to reduce the number of the equations when possible. But this makes the usersupplied program somewhat complicated. 4.5. Treatment of the linear equation If there are linear equations in the system, we can derive some simple relations among the coefficients of the variables so as to decrease the total number of equations, see Example 3. It is notable that this is the usual case in non-linear oscillations. 4.6. Use of symmetry In practice, we often encounter equations with a kind of symmetry, for example, the solutions of which do not include the harmonics of odd or even order, or of the form of sin nt or cos nt. Keeping these features in mind, we can reduce the unknown coefficients of the Galerkin approximation greatly. 4.7. Equations in implicit form. In some numerical methods it is necessary to express the equations in an explicit form (I), say, in the shooting type direct numerical method. In the Galerkin method the equations in implicit form g(ir, x, t) = 0 are treated directly. Then instead of equation (3) we have gC%,(s),x,(s),
slds= 0,
g[$,,(s), x,(s), s]sin ns ds = 0,
F2”_.i =;
s
2X F 2n = ;1 g[&,(s), x,(s), s]cos ns ds = 0. 0 In fact, higher-order derivatives can also be included in the function g. 5. APPLICATION
OF THE GALERKIN CALCULATIONS
METHOD TO AUTONOMOUS OF SUBHARMONICS
SYSTEMS
AND
In the above sections we described the application of the Galerkin method to equations with periodical excitation. As for an autonomous system t = f(x),
(6)
it is impossible to write the Fourier expansion of the periodic solution because the period
94
F. H. LIHG and X. X. Wu
is unknown. There are two different cases: the limit cycle, and the dense periodic solutions with various amplitudes and periods, which fill up the phase plane. A self-excited system, say, the van der Pol oscillator, belongs to the first category, and a conservative system, such as an undamped pendulum, belongs to the second category. In the first case, we have one extra unknown T, or equally, o = 2n/T. After a time transformation at = r, equation (6) becomes wx’ = f(x),
(7)
where the prime represents differentiation with respect to T. Then the periodic solution with an unknown frequency o in the time scale t becomes that with the known frequencyone in the time scale T. Equation (7) has a similar form to equation (2), but it contains an unknown parameter w, so the number of equations in equation (3) is one less than that of the unknown variables. But equation (6) represents an autonomous system whose solution is not distinguishable to the initial time, so we can always assume that there is an unknown coefficient a, which equals zero. Then the number of unknown variables is decreased by one, so that the problem becomes well defined. In the second case, we can always find a periodic solution corresponding to the given w in the same way as in the above section. But the solution is not isolated, so the abovementioned theoretical results are no longer applicable. However, very good results are still obtained with this method as shown in examples. In many non-linear oscillation systems, there are subharmonics of frequency w/n or super-subharmonics of frequency mm/n. In these cases, we make a time transformation wt = nt/m to transform the harmonics of frequency ma/n in the time scale t into basic periodic harmonics of frequency one in the time scale r. Thus, with the Galerkin method we can also treat subharmonics, see Example 4. Burton [9] suggested a time transformation approach in which an autonomous equation is transformed into a periodic one with known periods. But the method can be applied practically only to single-degree-of-freedom systems. 6. NUMERICAL
EXAMPLES
We have calculated many examples with the fast Galerkin method which show the efficiency of the method. We list some of them as follows. Example 1. Combination tones of a two-degree-of-freedom system 2, + 0.0021, + x1 + 0.002x,x2 + 0.003x: = 3.5~0~20~ + 6.5 cos 30t, 2, + 0.001512 + 0.81x, + 0.003x,x,
+ 0.002x: = 4cos20t
+ 13cos 3wt.
This problem has been investigated in [lo] with the perturbation method, and three Galerkin approximations for w = 0.92 are also given, one of which was revealed to be false in [l] with the direct numerical method. With m = 6 we calculated the coefficients of cos t, cos 2t and cos 3t for some w near 0.91 and depicted them in Figs 1 and 2. The other coefficients are small and are not presented here. The results we obtained agree very well with these given in [l]. Example 2. Limit cycles of van der Pol equation
This is an autonomous system which can be treated with the method described in the last section. Because of symmetry, there are no harmonics of even order. Taking a Galerkin approximation through the seventh order with an accuracy of lo-‘, we compared the results acquired with two methods-the fast Galerkin method (FGM) and the direct numerical method (DNM) Cl] in Table 1. It is seen that the accuracy is satisfactory. Example 3. Forced oscillations of an oscillator with two-degree-of-freedom. 22, + 3.477x, + 0.2x: + (x1 - x*) = 0.2cosor, 2, - (x1 - x*) = 0.
Fast Galerkin method
95
-
coszwt
---
cor3wt
ofx,
4-
2I-) 1
(-)
--.-_. t-
I 0.90
I
I
1
I
,
L
1
0.92
!
0.94
w
Fig. 1. Fourier coefficients of the Galerkin approximation of the solution x, of the equations f, + 0.0021, + x, + 0.002x,x, + 0.003x: = 3Scos2ot + 6Scos3wt 2, +0.0015x, + 0.81~~ + 0.003x,x, +0.002.x; = 4cos2rct + 13~0~3~~
0-
coswt
---cos2wt -
OfW* cos3wt
----
6-
4-
0.90
0.92
0.94 IJJ
Fig. 2. Fourier coefficients of the Galerkin approximation of the solution x2 of the equations 2, + 0.0@2& + x, + 0.002x,x, + 0.003~~ = 3.5c0s2ot + 6Scos3wt f2 + 0.00151, + 0.81~~ + 0.003x,x, +0.002x; = 4cosZot + 13cos3wt.
Table 1. Frequencies of the limit cycles of the van der Pol equation FGM DNM
0.1 0.2 0.3 0.4 0.5 0.999312 0.997510 0.994421 0.990143 0.984723 0.999316 0.997509 0.994420 0.990142 0.984721 0.6 0.7 0.8 0.9 1.0
FGM DNM
0.978219 0.970704 0.962260 0.952978 0.942960 0.978217 0.970101 0.962256 0.952975 0.942956
According to the fifth item in the last section, taking X1
=
,zliliCOS(2i -
~2 = izt
1)0X,
bicos(2i - l)~t,
F. H. LING and X. X. WU
96
from the second equation we have - (2i - l)‘fx’]
bi = aJ[l
(i = 1,2.. . .).
Then we solve the first equation, which is now a single degree-of-freedom system. By calculating only two coefficients a, and a2 a high accuracy is reached. The coefficients aI and bl are shown in Figs 3 and 4, in which the curves starting at U* = 0.678 and c_? = 2.56 are response curves for free oscillations of the unforced system 2f, + 3.477x, + 0.2x; + (x1 - x*) = 0, iz - (x1 - x2) = 0. The treatment is based on that described for the second case in the last section. Example 4. Subharmonics of the Duffing equation i + x + 0.1x3 = cos WC.
a1 -x1 -e-x*
6-
4-
0.76
0.68
0.60
Fig. 3. Resonance 2f,
curves near wz = 0.68 of the equations
+ 3.477x,
+ 0.2x: + (1, - XL) = 0.2cosor ~2-(.x,-.Yz)=o.
I
---x2
3 t
2.50
2.66
Fig. 4. Resonance ZP, + 3.477x,
curves
2.82
near (o* = 2.56 of the equations
+ 0.2.x: + (x, - x,) = 0.2coswt f, - (YI - x2) = 0.
20,
Fast Galerkin method
97
The subharmonic of l/3 order in the equation has been given in [l l] by using the perturbation method. When using the Galerkin method, we set wt = 3~, and transform the equation into the following form:
(1 W
T Y + x + 0.1x3 = cos 35.
Here the prime represents differentiation with respect to 5. Then the subharmonic of l/3 order transforms into the basic harmonic in the time scale ‘I. Let x
=
aI,3
cos
T +
a3,3
cos
35 + a5,3 cos 5r + a7,3 cos 75,
the subharmonic appears for w > 3, as shown in Fig. 5.
Fig. 5. Resonance curves of the subharmonic of l/3 order of the equation B + x + 0.1x’ = COSuJf.
7. CONCLUSIONS
There exists a strict theory for the existence and error-estimation of the Galerkin method. Introducing the FFT into it to form the so-called fast Galerkin method can greatly reduce the amount of computation. In addition, we extend the method to solve for the limit cycles in autonomous systems and subharmonics, etc. Besides we show that for a system with symmetry or with linear equations in it, the amount of computation can be reduced further. It is also mentioned that the Galerkin method can be used to treat implicit equations directly. A number of examples show that the fast Galerkin method is an effective tool to determine periodic solutions of a non-linear oscillator. It is remarkable that in order to find a satisfactory approximation, not only the residue EB corresponding to the harmonics, which is to be estimated, but also the EE corresponding to higher harmonics must be controlled. One of the authors has developed another numerical method, the shooting type direct numerical method Cl-33 for the same problem. The numerical examples indicate that both methods are applicable to complicated problems in non-linear oscillations, and can achieve a high accuracy. The advantages of the fast Galerkin method are the completeness in the theory, the clearness in physical meaning and the directness in error estimation. Especially for simpler problems or problems requiring only low accuracy, the amount of computation by the fast Galerkin method is less than that of the direct numerical method. REFERENCES 1. F. H. Ling, Numerische Berechnung periodischer Liisungen einiger nichtlinearer Schwingungssysteme. Dissertation, Uni. Stuttgart (1981). 2. F. H. Ling, Numerische Berechnung period&her Liisungen nichtlinearer Schwingungssysteme. 2. angew. Math. Mech. 62, T55-58 (1982).
98
F. H. LING and X. X. WV
3. F. H. Ling, A numerical treatment of the periodic solutions of non-linear vibration systems. Appl. Math. Mech.4, 525-546 (1983). 4. M. Urabe, Galerkin’s procedure for non-linear periodic systems. Arch. Rat. Mech. Anal. 20, 120-152 (1965). 5. M. Urabe and A. Reiter, Numerical computation of non-linear forced oscillations by Galerkin’s procedure. J. Math. Analysis Applic. 14, 107-140 (1966). 6. C. G. Broyden, A new method of solving non-linear simultaneous equations. Computer J. 12, 94-99 (1969). 7. E. 0. Brigham, The Fast Fourier Transform. Prentice-Hall, Englewood Cliffs, NJ (1974). 8. W. Gautsche, Attenuation factor in practical Fourier analysis. Num. Marh. 18, 373-400 (1972). 9. T. D. Burton, Non-linear oscillator limit cycle analysis using a time transformation approach. Int. J. Nonlinear Mech. 17, 7-19 (1982). 10. R. van Dooren, Differential tones in a damped mechanical system with quadratic and cubic non-linearities. Inr. J. Non-linear Mech. 8, 575-583 (1973). 11. J. J. Stoker, Non-linear Vibrarions in Mechanical and Electrical Systems. Interscience, New York (1950).
0
co2a-7462187 93.00 + .oo 1987 Pcrgamon Journals Ltd.
FAST GALERKIN METHOD AND ITS APPLICATION TO DETERMINE PERIODIC SOLUTIONS OF NON-LINEAR OSCILLATORS F. H. Department
of Engineering
LING
Mechanics,
(Received 9 August
and X. X. WV
Shanghai Jiao Tong University, P.R. China
1985; received for publication
23 September
Shanghai
200030,
1986)
Abstract-The Galerkin method is an approximate method which finds wide application in solving differential and integral equations. But a large amount of computation is needed in order to get a high order approximation by using the method. Applying the FFT technique to form a so-called fast Galerkin method, we can reduce the computation work greatly, when taking trigonometric functions as characteristic functions. Taking the periodic solution of non-linear oscillators as an example, we illustrate the procedure and the efficiency of the method. Moreover, with some modifications we extend the applicability of the method, so that not only periodic solutions with known periods, but also those with unknown periods, as well as subharmonics, combination tones, etc., can be treated with the method. Some techniques are described which can be used to simplify the computation. I. INTRODUCTION
The Galerkin method is an approximate method which is often applied to solve differential equations and integral equations numerically. The basic idea is to seek a solution in terms of a series of characteristic functions which satisfy the boundary conditions. For partial differential equations this procedure leads to a system of ordinary differential equations, which can be treated with, say, the Galerkin method once again. As for ordinary differential equations or integral equations, it leads to a system of algebraic equations which the coefficients of the series must satisfy. Solving these equations, usually by a numerical method, we get a Galerkin approximate solution of the equations. In this paper, we regard the periodic solution of non-linear oscillators as a boundary value problem of ordinary differential equations (as in [l-3]), but the method described here can also be applied to other boundary value problems as well as integral equations. The problem of determining a periodic solution can be put forward as follows, Ti = f(x, t),
XER',
f: UP+1- R,
f(x, t + 2x) = f(x),
x(0) - x(279 = 0,
(1)
where f(x, t) is periodic in t with a period of 271without losing generality. Choosing the following trigonometric functions: {1, sin nt, cos nt>, n = 1,2,. . . , as characteristic functions, we obtain a Galerkin approximation of order m:
x,(0 = a0 + f (a2”_ 1sin nt n=1
+ a2,cos nt),
(2)
where x,(t), ao, a+ I, azn (n = 1, 2,. . . , m) are all I-dimensional vectors. Urabe [4,5] has investigated the above problem. His main results are as follows. Assume that f(x, t) and its derivatives with respect to the x-coordinate are continuously differentiable with respect to the same x-coordinate and t in a closed bounded region of the x-space and the axis of t. If there is an isolated periodic solution x = n(t), there exists a Galerkin approximation x = x,,,(t), which converges uniformly as m + co to the exact solution x = n(t). On the other hand, if the Galerkin approximation x = x,(t) satisfies certain conditions, there must be an exact isolated periodic solution lying in the neighborhood of xl&). Urabe applied the classic Fourier transform (the FFT was not discovered at that time), which requires a large amount of computation, especially for obtaining a higher order 89
F. H. LINCand X. X. Wu
90
approximation, so as to restrict its usage greatly. In this paper we introduce the FFT to form the so-called fast Galerkin method, in order to reduce the amount of computation and hence to raise its applicability. Besides, with some modifications the method can be used to solve not only the periodic solutions with known periods but also those with unknown periods, as well as to get subharmonics and combination tones, etc. We illustrate the efficiency in some numerical examples. Some techniques which can be used to reduce the amount of computation are also listed here. 2. PROCEDURE
FOR
DETERMINING PERIODIC GALERKIN METHOD
SOLUTIONS
BY USING
THE
Substituting the expression (2) into equation (1) and equating the coefficients of sin nt and cosnt on both sides of the equation, we have
FO(4
fcxm(s),
= &
slds = 0,
s 2n
F2,-l(4
= i
Rx,,,(s), s]sin
ns ds + na2,, = 0,
0
fCx,(s),s]cosns
F2,(a) = a
ds - nazn_ i = 0,
(n = 1,2,...,m)
(3)
where a = [a$, a:, . . . , ai,,J’ is a (2m + 1) x I-dimensional vector, which gathers all unknown coefficients. The equations (3) are non-linear algebraic equations for determining a and can be solved with, say, a Newton-Raphson iterative method. In this case a system of linear equations for the correctors Aa of the unknown coefficients results: JAa + F = 0,
(4)
where F = [Fc, FT,. . . , FT,,,]’ is a (2m + 1) x I-dimensional vector, and J = dF/t?a is the first order derivative matrix, whose elements are definite integrals of the products of the elements of the derivative $ = af/a x with cosnt and/or sinnt and can be regarded as the Fourier coefficients of the elements of +. Some other methods, for example, the Broyden method [6], can also be used to solve equations (3); the corresponding iterative scheme takes the form of p(il = =
p(i),
yW =
F(i+
A,(i)
H(‘+
_ H(i)F(i)
1) = H”’
,
1) _
F(i),
_
(i) ,, (i) _
(H
(il
(i)rH(i)/p(i)TH(i)y(i).
P )P
(5)
Thereby, we sum up the procedure of using the Galerkin method to determine a periodic solution as follows: (1) Estimate the initial value a(‘) of the coefficient vector a in expression (2); (2) Select 2m, equidistant nodes ti, i = 1, 2,. . . , 2me, in the time interval (0,2n), then repeat the following steps with k = 0, l,-2,. . .: (a) Evaluate x$‘(ti) from a (kl by using the inverse Fourier transform; (b) Calculate fCX~‘(fi),ti], and when using the Newton-Raphson method, also calculate
+lIx!tt’(ti)9 lil; (c)Calculate J
from $[X$‘(tJ, ti] with the Fourier transform, when using the NewtonRaphson method; (d) Evaluate F from fCX~‘(Ci),ti] with the Fourier transform according to equations (3); (e) Solve equation (4) when using the Newton-Raphson method, or equation (5) when
Fast Galerkin method
91
using the Broyden method to determine correctors Aa”‘; (f) End the iteration, if r is small enough, otherwise set z(~+I) = zCk)-I- A,(‘) and return to step (a). 3. FORMING
THE FAST GALERKIN METHOD TECHNIQUE
BY INTRODUCING
THE FFT
From the procedure stated above, we see that the computational work in the Galerkin method comes from three source% the solution of the algebraic equations [step2(e)], the calculation of f and $ at nodes ti [step2(b)], and the Fourier transforms [steps 2(a), (c) and (d)]. For the first one, there are many effective methods, such as the Broyden method and the super relaxation iterative method, which both require 0(4m*I*) times of multiplications, where m is the order of the harmonics to be found. For the second one, generally snaking, the number of multiplications is of O(Zm,& which varies with the complexity of the equation, where m, is the order of the Fourier transform. We must have m, > 2m for a reason to be stated later. For the third one, the classical discrete Fourier transform requires a large amount of computation. In the steps 2(a), (c)and (d), the numbers of multiplications are 1(2m,)2, f2(2m,)2 and I(2mJ2 respectively, and their sum is (I f 2)t(2mJ2. Introducing the FFT technique (see [?‘J, for example) into the Galerkin method may reduce the number of multiplications to l/Z(I + 2~I(2m~)log~(2m~) with the FFT algo~thm of index 2, i.e. only 18% of the original computer work for m, = 4, 12.5% for m, = 8, and 7.8% for m, = 16. Now we explain the reason of taking m, 3 2m,where m is the order of the harmonics to be found and m, is the order of the Fourier transform. As is well known, the coefficients of higher harmonics evaluated with the Fourier transform are no longer valid, and usually only the first half of the m, coefficients can be used. On the other hand, in order to reach a high degree of accuracy, the value of r, which is equal or greater than )/P(t) - f@(t), t)ll, must be kept small enough, where S(t) is an approximate solution. Expanding (k(t) - f@(t), t)) into a Fourier series, we have a(t) - f@(r), t) = A,, f
i
(Azn_ I sin nt + Al,, cosat),
&?=1
which is equal approximately
to
k(t) - f@(t), t) = A0 + 3 (Azn-. 1sin nt + Aln cos d. n-1
From equation (3), we find that (i = O,l,... ,2m).
Ai =z: Fi, Let
EB = m,ax A,, + 2 (A,,_ I sin nt + Az,cosnt) II
EE = rn,“x f$+ li
,
n-1 1(AZ,--I sin nt -I- A2*cos nt) , 1/
then we can take r = EB + EE. The practical calculation shows that the control of EE is very important. See next section for details. This requires that m, 3 2m, otherwise EE cannot be determined. In this way we see that the Fourier transform is a time consuming work in the method and the application of the FFT technique reduces the amount of computation greatly. Taking the method with the Broyden method as an example, in a two-degree-of-freedom
F. H, LING and X. X. WV
92
system (I = 4), we find that to determine the harmonics through the 8th order (m = 8, taking m, = 16), the above three kinds of computation need 4 x 8’ x 42 = 4096,2 x 16 x 4 = 128, and 2 x 4 x 322 = 8192 times of multiplications respectively, and altogether 12,416 times. After adopting the FFT, the third one decreases to 2 x 4 x 32 x log,32 = 1280 times, and the sum reduces to 5504 times, which is only 44.3% of the original. In compiling the program we take m, = 2m, and EE and EB are enlarged to
J
A; + f (At,It=1
EB =
EE =
J
z
n=m+ 1
(A;,-,
1 +
Ain)
+ A:,)
respectively. 4. SOME
REMARKS
4.1. Choice of the initial value a(” of the coejkient vector a This is an important problem which has a strong effect on the convergence speed and even the convergence itself. The commonly used methods are: (a) Make an extrapolation of the solutions obtained to the equations with slightly different coefficients to give the initial value of the solution which is sought. For example, we can make use of the solution with neighbouring forcing frequencies to determine the response curves. For weak non-linear systems we can make use of the solution of a corresponding linearized system. (b) Use a lower-order Galerkin approximation of the same system. A lower-order Galerkin approximation is denoted by x = a,* + z (a&_, sinnt + ar,cos nt). PI=1
Thus the initial value may be chosen as (0) =
%
a? 1
a!O) I = 0.
I
(i = 0,1,...,2m,) (i = 2m, + 1,. . . ,2m,).
4.2. Necessity of controlling the residue EE The well-known harmonic balance method may be regarded as a low-order Galerkin method, in which only the conditions Fi = 0 (i = 1,. . . , m), i.e. EB = 0 are demanded. But according to our experience gathered in the computation we establish that this is not sufficient. It may happen that an approximate solution with a small EB but a large EE is far away from the exact solution. So we must check the condition of EE < EPS after the condition of EB < EPS is satisfied. If the condition EE < EPS is not satisfied, we should increase the order m, and if the error EE is still unacceptable even with a large m, then there seems to be no proper periodic solution nearby. 4.3. A comparison between the Newton-Raphson here
method and the Broyden method used
Both methods are tested in -the computation. With the Newton-Raphson method 0(9m313)times of multiplications are required for solving the linear algebraic equation, and the derivative matrix J has to be evaluated. So the program is complicated, and the numerical stability in the computation is poor, i.e. it requires a fairly good initial guess, otherwise it would diverge. The advantage of the quick convergence of the NewtonRaphson method cannot compensate for these disadvantages. With the Broyden method the computation of J is avoidable, and we obtain a concise program. Although it converges more slowly (in about 20 iterations), the total amount of computation is less than that
93
Fast Galerkin method
with the Newton-Raphson method, because it needs only 0(4m21*) times of multiplications for every iteration. Moreover, the numerical stability of the Broyden method is better than that of Newton-Raphson method, i.e. the convergence of the method is nearly always guaranteed even with a poor initial guess. 4.4. Treatment of a second order equation Most non-linear oscillation equations are of second order and need not be transformed into the standard form (1). We substitute the expression (2) directly into the equations, by noticing that k,(t) = “cl n(a2”- l cos nt - a2” sin nt)
In this way, the number of unknown coefficients is only half of that for the transformed first-order equations, so the amount of computation is reduced to nearly one-fourth. This idea can be extended to a more general case, in which the equations are of higher order. We even transform a set of first-order equations into a set of higher-order equations in order to reduce the number of the equations when possible. But this makes the usersupplied program somewhat complicated. 4.5. Treatment of the linear equation If there are linear equations in the system, we can derive some simple relations among the coefficients of the variables so as to decrease the total number of equations, see Example 3. It is notable that this is the usual case in non-linear oscillations. 4.6. Use of symmetry In practice, we often encounter equations with a kind of symmetry, for example, the solutions of which do not include the harmonics of odd or even order, or of the form of sin nt or cos nt. Keeping these features in mind, we can reduce the unknown coefficients of the Galerkin approximation greatly. 4.7. Equations in implicit form. In some numerical methods it is necessary to express the equations in an explicit form (I), say, in the shooting type direct numerical method. In the Galerkin method the equations in implicit form g(ir, x, t) = 0 are treated directly. Then instead of equation (3) we have gC%,(s),x,(s),
slds= 0,
g[$,,(s), x,(s), s]sin ns ds = 0,
F2”_.i =;
s
2X F 2n = ;1 g[&,(s), x,(s), s]cos ns ds = 0. 0 In fact, higher-order derivatives can also be included in the function g. 5. APPLICATION
OF THE GALERKIN CALCULATIONS
METHOD TO AUTONOMOUS OF SUBHARMONICS
SYSTEMS
AND
In the above sections we described the application of the Galerkin method to equations with periodical excitation. As for an autonomous system t = f(x),
(6)
it is impossible to write the Fourier expansion of the periodic solution because the period
94
F. H. LIHG and X. X. Wu
is unknown. There are two different cases: the limit cycle, and the dense periodic solutions with various amplitudes and periods, which fill up the phase plane. A self-excited system, say, the van der Pol oscillator, belongs to the first category, and a conservative system, such as an undamped pendulum, belongs to the second category. In the first case, we have one extra unknown T, or equally, o = 2n/T. After a time transformation at = r, equation (6) becomes wx’ = f(x),
(7)
where the prime represents differentiation with respect to T. Then the periodic solution with an unknown frequency o in the time scale t becomes that with the known frequencyone in the time scale T. Equation (7) has a similar form to equation (2), but it contains an unknown parameter w, so the number of equations in equation (3) is one less than that of the unknown variables. But equation (6) represents an autonomous system whose solution is not distinguishable to the initial time, so we can always assume that there is an unknown coefficient a, which equals zero. Then the number of unknown variables is decreased by one, so that the problem becomes well defined. In the second case, we can always find a periodic solution corresponding to the given w in the same way as in the above section. But the solution is not isolated, so the abovementioned theoretical results are no longer applicable. However, very good results are still obtained with this method as shown in examples. In many non-linear oscillation systems, there are subharmonics of frequency w/n or super-subharmonics of frequency mm/n. In these cases, we make a time transformation wt = nt/m to transform the harmonics of frequency ma/n in the time scale t into basic periodic harmonics of frequency one in the time scale r. Thus, with the Galerkin method we can also treat subharmonics, see Example 4. Burton [9] suggested a time transformation approach in which an autonomous equation is transformed into a periodic one with known periods. But the method can be applied practically only to single-degree-of-freedom systems. 6. NUMERICAL
EXAMPLES
We have calculated many examples with the fast Galerkin method which show the efficiency of the method. We list some of them as follows. Example 1. Combination tones of a two-degree-of-freedom system 2, + 0.0021, + x1 + 0.002x,x2 + 0.003x: = 3.5~0~20~ + 6.5 cos 30t, 2, + 0.001512 + 0.81x, + 0.003x,x,
+ 0.002x: = 4cos20t
+ 13cos 3wt.
This problem has been investigated in [lo] with the perturbation method, and three Galerkin approximations for w = 0.92 are also given, one of which was revealed to be false in [l] with the direct numerical method. With m = 6 we calculated the coefficients of cos t, cos 2t and cos 3t for some w near 0.91 and depicted them in Figs 1 and 2. The other coefficients are small and are not presented here. The results we obtained agree very well with these given in [l]. Example 2. Limit cycles of van der Pol equation
This is an autonomous system which can be treated with the method described in the last section. Because of symmetry, there are no harmonics of even order. Taking a Galerkin approximation through the seventh order with an accuracy of lo-‘, we compared the results acquired with two methods-the fast Galerkin method (FGM) and the direct numerical method (DNM) Cl] in Table 1. It is seen that the accuracy is satisfactory. Example 3. Forced oscillations of an oscillator with two-degree-of-freedom. 22, + 3.477x, + 0.2x: + (x1 - x*) = 0.2cosor, 2, - (x1 - x*) = 0.
Fast Galerkin method
95
-
coszwt
---
cor3wt
ofx,
4-
2I-) 1
(-)
--.-_. t-
I 0.90
I
I
1
I
,
L
1
0.92
!
0.94
w
Fig. 1. Fourier coefficients of the Galerkin approximation of the solution x, of the equations f, + 0.0021, + x, + 0.002x,x, + 0.003x: = 3Scos2ot + 6Scos3wt 2, +0.0015x, + 0.81~~ + 0.003x,x, +0.002.x; = 4cos2rct + 13~0~3~~
0-
coswt
---cos2wt -
OfW* cos3wt
----
6-
4-
0.90
0.92
0.94 IJJ
Fig. 2. Fourier coefficients of the Galerkin approximation of the solution x2 of the equations 2, + 0.0@2& + x, + 0.002x,x, + 0.003~~ = 3.5c0s2ot + 6Scos3wt f2 + 0.00151, + 0.81~~ + 0.003x,x, +0.002x; = 4cosZot + 13cos3wt.
Table 1. Frequencies of the limit cycles of the van der Pol equation FGM DNM
0.1 0.2 0.3 0.4 0.5 0.999312 0.997510 0.994421 0.990143 0.984723 0.999316 0.997509 0.994420 0.990142 0.984721 0.6 0.7 0.8 0.9 1.0
FGM DNM
0.978219 0.970704 0.962260 0.952978 0.942960 0.978217 0.970101 0.962256 0.952975 0.942956
According to the fifth item in the last section, taking X1
=
,zliliCOS(2i -
~2 = izt
1)0X,
bicos(2i - l)~t,
F. H. LING and X. X. WU
96
from the second equation we have - (2i - l)‘fx’]
bi = aJ[l
(i = 1,2.. . .).
Then we solve the first equation, which is now a single degree-of-freedom system. By calculating only two coefficients a, and a2 a high accuracy is reached. The coefficients aI and bl are shown in Figs 3 and 4, in which the curves starting at U* = 0.678 and c_? = 2.56 are response curves for free oscillations of the unforced system 2f, + 3.477x, + 0.2x; + (x1 - x*) = 0, iz - (x1 - x2) = 0. The treatment is based on that described for the second case in the last section. Example 4. Subharmonics of the Duffing equation i + x + 0.1x3 = cos WC.
a1 -x1 -e-x*
6-
4-
0.76
0.68
0.60
Fig. 3. Resonance 2f,
curves near wz = 0.68 of the equations
+ 3.477x,
+ 0.2x: + (1, - XL) = 0.2cosor ~2-(.x,-.Yz)=o.
I
---x2
3 t
2.50
2.66
Fig. 4. Resonance ZP, + 3.477x,
curves
2.82
near (o* = 2.56 of the equations
+ 0.2.x: + (x, - x,) = 0.2coswt f, - (YI - x2) = 0.
20,
Fast Galerkin method
97
The subharmonic of l/3 order in the equation has been given in [l l] by using the perturbation method. When using the Galerkin method, we set wt = 3~, and transform the equation into the following form:
(1 W
T Y + x + 0.1x3 = cos 35.
Here the prime represents differentiation with respect to 5. Then the subharmonic of l/3 order transforms into the basic harmonic in the time scale ‘I. Let x
=
aI,3
cos
T +
a3,3
cos
35 + a5,3 cos 5r + a7,3 cos 75,
the subharmonic appears for w > 3, as shown in Fig. 5.
Fig. 5. Resonance curves of the subharmonic of l/3 order of the equation B + x + 0.1x’ = COSuJf.
7. CONCLUSIONS
There exists a strict theory for the existence and error-estimation of the Galerkin method. Introducing the FFT into it to form the so-called fast Galerkin method can greatly reduce the amount of computation. In addition, we extend the method to solve for the limit cycles in autonomous systems and subharmonics, etc. Besides we show that for a system with symmetry or with linear equations in it, the amount of computation can be reduced further. It is also mentioned that the Galerkin method can be used to treat implicit equations directly. A number of examples show that the fast Galerkin method is an effective tool to determine periodic solutions of a non-linear oscillator. It is remarkable that in order to find a satisfactory approximation, not only the residue EB corresponding to the harmonics, which is to be estimated, but also the EE corresponding to higher harmonics must be controlled. One of the authors has developed another numerical method, the shooting type direct numerical method Cl-33 for the same problem. The numerical examples indicate that both methods are applicable to complicated problems in non-linear oscillations, and can achieve a high accuracy. The advantages of the fast Galerkin method are the completeness in the theory, the clearness in physical meaning and the directness in error estimation. Especially for simpler problems or problems requiring only low accuracy, the amount of computation by the fast Galerkin method is less than that of the direct numerical method. REFERENCES 1. F. H. Ling, Numerische Berechnung periodischer Liisungen einiger nichtlinearer Schwingungssysteme. Dissertation, Uni. Stuttgart (1981). 2. F. H. Ling, Numerische Berechnung period&her Liisungen nichtlinearer Schwingungssysteme. 2. angew. Math. Mech. 62, T55-58 (1982).
98
F. H. LING and X. X. WV
3. F. H. Ling, A numerical treatment of the periodic solutions of non-linear vibration systems. Appl. Math. Mech.4, 525-546 (1983). 4. M. Urabe, Galerkin’s procedure for non-linear periodic systems. Arch. Rat. Mech. Anal. 20, 120-152 (1965). 5. M. Urabe and A. Reiter, Numerical computation of non-linear forced oscillations by Galerkin’s procedure. J. Math. Analysis Applic. 14, 107-140 (1966). 6. C. G. Broyden, A new method of solving non-linear simultaneous equations. Computer J. 12, 94-99 (1969). 7. E. 0. Brigham, The Fast Fourier Transform. Prentice-Hall, Englewood Cliffs, NJ (1974). 8. W. Gautsche, Attenuation factor in practical Fourier analysis. Num. Marh. 18, 373-400 (1972). 9. T. D. Burton, Non-linear oscillator limit cycle analysis using a time transformation approach. Int. J. Nonlinear Mech. 17, 7-19 (1982). 10. R. van Dooren, Differential tones in a damped mechanical system with quadratic and cubic non-linearities. Inr. J. Non-linear Mech. 8, 575-583 (1973). 11. J. J. Stoker, Non-linear Vibrarions in Mechanical and Electrical Systems. Interscience, New York (1950).