Vibrational resonance: a study with high-order word-series averaging

We study a model problem describing vibrational resonance [12] by means of the high-order averaging technique introduced in the series of articles [5], [6], [7], [8] (see [15] for a summary).

In devices with vibrational resonance [12], [2], [16] the response of a system driven by a low-frequency forcing may be enhanced by the presence of high-frequency vibrations of suitable amplitude. Such devices feature in several current applications, including energy harvesting [9], [10], where the aim is to exploit the energy contained in the vibrations by converting it into electrical energy. The corresponding systems of differential equations are highly oscillatory in the sense that their solutions include periods much shorter than the time interval of interest. The numerical integration of such systems may be a very expensive task, as, typically, numerical integrators have to operate with step sizes significantly shorter than the shortest period present in the solution and the simulation has to be carried out over many periods. The literature contains many suggestions of nonstandard integrators to be applied in highly oscillatory situations, including the heterogeneous multiscale method [11], [19], methods with multiple time steps [18], multirevolution techniques [3], etc. However the use of such ad hoc integrators is prone to unexpected instabilities and inaccuracies [4].

Averaging techniques are very useful to deal with oscillatory problems; they may be applied in combination with numerical integration, to ‘soften’ the difficulties of the system to be simulated numerically, or in a purely analytical way. There are many techniques to carry out the computations required by the method of averaging (see e.g. [17]); here we use the word series [20] approach introduced in [5], [6], [7], [8]. With the technique applied here, the tasks of constructing the averaged system and the associated change of variables are divided into two parts. It is first necessary to build recursively a set of so-called word basis functions and, after that, all the required manipulations involve only scalar coefficients that are computed by means of simple recursions. As distinct from the situation with other approaches, with word-series high-order averaged systems may be derived without having to compute the associated change of variables. The use of word series techniques is not confined to averaging; word series may be applied to compute normal forms of discrete or continuous dynamical systems [13], [14], [15], find formal invariants of differential systems [14], analyze numerical integrators [13], [1], etc.

Section 2 provides a brief summary of the averaging technique we employ. It should be emphasized that while, for the sake of simplicity, this summary is restricted to the case of periodic forcing, the technique may be applied to quasiperiodic problems without essential changes [6], [15]. Section 3 studies a model problem taken from the original reference on vibrational resonance by Landa and McClintock [12]. We show how averaging provides insight into the mechanism causing vibrational resonance; this mechanism is related to the creation of an effective potential, akin e.g. to that responsible for stabilizing Kapitza’s inverted pendulum [19]. Furthermore, the construction of high-order averaged systems makes it possible to get very precise approximations to the true dynamics of the system under investigation.

The example here is based on the overdamped oscillator; the technique may be equally applied to under-damped models [9], [10]. Work to extend our study to realistic devices that exhibit vibrational resonance is under way.

After Fourier expansion, we may assume that the periodic problem to be studied is of the form:

ddtx=∑k=−∞∞exp(ikωt)fk(x),x(t0)=x0∈ℝD.$$\begin{array}{} \displaystyle \frac{d}{{dt}}x = \sum\limits_{k = - \infty }^\infty {{\rm{exp}}(ik\omega t){f_k}(x)} ,\qquad x\left( {{t_0}} \right) = {x_0} \in {\mathbb{R}^D}. \end{array}$$(1)

In the averaging technique used here, the elements (indices) k = 0,±1,±2,... in (1) are seen as the letters of an alphabet and each (possibly empty) string of letters k₁k₂ ...k_n, n = 0,1,..., is called a word. The symbol 𝒲 is used to denote the set of all words. With each word w ∈ 𝒲 we associate a word basis function f_w. For the empty word, f_∅(x) is the identity map x ↦ x in ℝ^D. For nonempty words the f_w(x) are constructed recursively from the f_k(x) that feature in (1). The recipe is

fk1…kn(x)=fk2…kn′(x)fk1(x),$$\begin{array}{} \displaystyle {f_{{k_1} \ldots {k_n}}}(x) = f_{{k_2} \ldots {k_n}}^\prime (x){f_{{k_1}}}(x), \end{array}$$

where fk2…kn′(x)$\begin{array}{} \displaystyle f_{{k_2} \ldots {k_n}}^\prime (x) \end{array}$ is the Jacobian matrix of f_k2...kn(x). We define ℂ^𝒲 as the set (vector space) of all mappings δ : 𝒲 → ℂ; if δ ∈ ℂ^𝒲 and w ∈ 𝒲, δ_w represents the complex number that δ associates with w. To each δ ∈ ℂ^𝒲 there corresponds a word series

Wδ(x)=∑w∈Wδwfw(x).$$\begin{array}{} \displaystyle {W_\delta }(x) = \mathop \sum \limits_{w \in {\mathscr W}} {\delta _w}{f_w}(x). \end{array}$$

With this terminology, the main result from [6] is as follows. Consider the averaged problem

ddtX=Wβ¯(t0)(X),X(t0)=x0$$\begin{array}{} \displaystyle \frac{d}{{dt}}X = {W_{\bar \beta ({t_0})}}(X),\qquad X({t_0}) = {x_0} \end{array}$$(2)

and the time-dependent change of variables

x=Wκ(tω,t0)(X),$$\begin{array}{} \displaystyle x = {W_{\kappa (t\omega ,{t_0})}}(X), \end{array}$$(3)

where the coefficients β¯w(t0),κw(tω,t0),w∈W$\begin{array}{} \displaystyle {\bar \beta _w}({t_0}),{\kappa _w}(t\omega ,{t_0}),w \in {\mathscr W} \end{array}$, necessary to write the word series may be computed explicitly (see below). Then the solution x(t) of (1) may be represented as

x(t)=Wκ(tω,t0)(X(t)),$$\begin{array}{} \displaystyle x(t) = {W_{\kappa (t\omega ,{t_0})}}(X(t)), \end{array}$$(4)

where X(t) is the solution of (2). The change of variables (3) is 2π/ω periodic in t and, in addition, at the stroboscopic times t_ℓ = t₀ + ℓ(2π/ω), x and X coincide. In general the series in the right hand-sides of (2) and (3) do not converge and have to be truncated as demonstrated in the next section. With the help of such truncations it is possible to approximate x(t) accurately without solving the oscillatory problem (1).

The coefficients (β¯w(t0))$\begin{array}{} \displaystyle ({\bar \beta _w}({t_0})) \end{array}$ are computed recursively by means of the formulas

β¯k(t0)=0,β¯0(t0) =1,β¯0r+1(t0) =0,β¯0rk(t0)=ikω(β¯0r−1k(t0)−β¯0r(t0)eikωt0),β¯kℓ1⋯ℓs(t0)=ikω(eikωt0β¯ℓ1⋯ℓs(t0)−β¯(k+ℓ1)ℓ2⋯ℓs(t0)),β¯0rkℓ1⋯ℓs(t0)=ikω(β¯0r−1kℓ1⋯ℓs(t0)−β¯0r(k+ℓ1)ℓ2⋯ℓs(t0)).$$\begin{array}{} \displaystyle \begin{array}{*{20}{c}} {{{\bar \beta }_k}({t_0})} \hfill & { = 0,} \hfill \\ {{{\bar \beta }_0}({t_0})} \hfill & { = 1,} \hfill \\ {{{\bar \beta }_{{0^{r + 1}}}}({t_0})} \hfill & { = 0,} \hfill \\ {{{\bar \beta }_{{0^r}k}}({t_0})} \hfill & { = \frac{i}{{k\omega }}({{\bar \beta }_{{0^{r - 1}}k}}({t_0}) - {{\bar \beta }_{{0^r}}}({t_0}){e^{ik\omega {t_0}}}),} \hfill \\ {{{\bar \beta }_{k{\ell _1} \cdots {\ell _s}}}({t_0})} \hfill & { = \frac{i}{{k\omega }}({e^{ik\omega {t_0}}}{{\bar \beta }_{{\ell _1} \cdots {\ell _s}}}({t_0}) - {{\bar \beta }_{(k + {\ell _1}){\ell _2} \cdots {\ell _s}}}({t_0})),} \hfill \\ {{{\bar \beta }_{{0^r}k{\ell _1} \cdots {\ell _s}}}({t_0})} \hfill & { = \frac{i}{{k\omega }}({{\bar \beta }_{{0^{r - 1}}k{\ell _1} \cdots {\ell _s}}}({t_0}) - {{\bar \beta }_{{0^r}(k + {\ell _1}){\ell _2} \cdots {\ell _s}}}({t_0})).} \hfill \\ \end{array} \end{array}$$

Here the integer k is ≠ 0 and 0^r, r > 0, denotes the word consisting of r zeros. Note that, from these formulas, (β¯w)$\begin{array}{} \displaystyle ({\bar \beta _w}) \end{array}$ is of size 𝒪(1/ωⁿ⁻¹) for n letter words w, n > 0. The, very similar formulas for the coefficients κ(tω,t₀) may be seen in [15].

The following one-dimensional overdamped double-well oscillator

dzdt=z−z3+Acosνt+Bωcosωt,z(0)=z0,$$\begin{array}{} \displaystyle \frac{{dz}}{{dt}} = z - {z^3} + A\cos \nu t + B\omega \cos \omega t,\qquad z(0) = {z_0}, \end{array}$$(5)

has been considered in [12] as a simple model to demonstrate vibrational resonance. Here 0 < v ≪ ω, the parameter A represents the amplitude of the applied driving force and B measures the size of the fast background vibration. We note that in [12] the vibrational term is written in the alternative format C cos ωt rather than as Bω cos ωt. Since no hypotheses are made in what follows as to the size of B, both formats are equivalent; when C and ω are given specific numerical values the problem considered in [12] may be cast in the form (5) by defining B = C/ω. However writing the amplitude of the vibration as Bω is more meaningful than writing it as C because if one sees ω as a parameter and keeps C constant then it is clear that for ω sufficiently large the vibrational term C cos ωt cannot be expected to exert any significant influence as it converges weakly to zero.

We begin by removing from (5) the 𝒪(ω) term by means of the preliminary change of variables z = x + B sin ωt (z and y coincide at the stroboscopic times ℓ(2π/ω)), which leads to

dydt=y−32B2y−y3+Acos⁡νt+B1−34B2−3y2sin⁡ωt+32B2ycos⁡2ωt+14B3sin⁡3ωt.$$\begin{array}{l} \frac{{dy}}{{dt}}{\rm{ }} = & {\rm{ }}y - \frac{3}{2}{B^2}y - {y^3} + A\cos \nu t\\ & + B\left( {1 - \frac{3}{4}{B^2} - 3{y^2}} \right)\sin \omega t + \frac{3}{2}{B^2}y\cos 2\omega t + \frac{1}{4}{B^3}\sin 3\omega t. \end{array}$$(6)

Note that the conservative force z − z³ in (5) has been ‘softened’ to (1 − 3B²/2)y − y³ in (6). In terms of the corresponding potentials, we have moved from −z²/2 + z⁴/4, with wells of depth 1/4 at z = ±1 to the new symmetric effective potential

(−12+34B2)y2+14y4.$$\begin{array}{} \displaystyle \left( { - \frac{1}{2} + \frac{3}{4}{B^2}} \right){y^2} + \frac{1}{4}{y^4}. \end{array}$$

As |B| increases from B = 0, the depth of the wells decreases and for (|B|≥2/3)$\begin{array}{} \displaystyle (|B| \ge \sqrt {2/3} ) \end{array}$ the two potential minima are merged into a single minimum at y = 0.

By considering an auxiliary variable ϕ with (d/dt)ϕ = v and the vector x = (ϕ,y)^T ∈ ℝ², (6) takes the form (1) studied in the preceding section with

f0(x)=(ν,y−32B2y−y3+Acosϕ)T,f1(x)=(0,−i2B+3i8B3+3i2By2)T,f2(x)= (0,34B2y)T,f3(x)= (0,−i8B3)T,$$\begin{array}{} \displaystyle \begin{array}{*{20}{c}} {{f_0}(x)} \hfill & = \hfill & {{{(\nu ,y - \frac{3}{2}{B^2}y - {y^3} + A\cos \phi )}^T},} \hfill \\ {{f_1}(x)} \hfill & = \hfill & {{{(0, - \frac{i}{2}B + \frac{{3i}}{8}{B^3} + \frac{{3i}}{2}B{y^2})}^T},} \hfill \\ {{f_2}(x)} \hfill & = \hfill & {{{(0,\frac{3}{4}{B^2}y)}^T},} \hfill \\ {{f_3}(x)} \hfill & = \hfill & {{{(0, - \frac{i}{8}{B^3})}^T},} \hfill \\ \end{array} \end{array}$$

f_k(x) equal to the complex conjugate of f_−k(x) for k = −1,−2,−3, and f_k(x) = 0 for |k| > 3.

Truncating the averaged system (2) to keep only contributions corresponding to words of one letter results in

ddt[ΦY]=[νY−32B2Y−Y3+AcosΦ],$$\begin{array}{} \displaystyle \frac{d}{{dt}}\left[ {\begin{array}{*{20}{c}} \Phi \\ Y \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} \nu \\ {Y - \frac{3}{2}{B^2}Y - {Y^3} + A\cos \Phi } \end{array}} \right], \end{array}$$

i.e. the effect of averaging is only to remove the terms of (6) that involve the background vibration. A simple computation of the required basis functions and coefficients that may be carried out by hand shows that truncating after two letter words yields (d/dt)Φ = v (as expected) and

ddtY=Y−32B2Y−Y3+Acos⁡Φ+1ωB−136B3+B5−52B3Y2+3BY4+6ABYcos⁡Φ.$$\begin{array}{l} \frac{d}{{dt}}Y{\rm{ }} = & Y - \frac{3}{2}{B^2}Y - {Y^3} + A\cos \Phi \\ &+ \frac{1}{\omega }\left( {B - \frac{{13}}{6}{B^3} + {B^5} - \frac{5}{2}{B^3}{Y^2} + 3B{Y^4} + 6ABY\cos \Phi } \right). \end{array}$$(7)

Thus the addition of two-letter words introduces corrections of size 𝒪(1/ω). The new effective potential

−12Y2+34B2Y2+14Y3+1ω(−BY+136B3Y−B5Y+56B3Y3−35BY5)$$\begin{array}{} \displaystyle - \frac{1}{2}{Y^2} + \frac{3}{4}{B^2}{Y^2} + \frac{1}{4}{Y^3} + \frac{1}{\omega }\left( { - BY + \frac{{13}}{6}{B^3}Y - {B^5}Y + \frac{5}{6}{B^3}{Y^3} - \frac{3}{5}B{Y^5}} \right) \end{array}$$

has lost the Y ↔ −Y symmetry and the amplitude A of the slow driving forcing has been incremented by an amount 6ABY/ω.

Due to the absence of fast vibrations, (7) may be integrated numerically with high accuracy with a negligible computational effort. We have done so for the values (taken from [12]) A = 0.2, v = 0.1, ω = 5. The initial condition is set to be −1, so that the oscillator starts at one of the minima of the potential well of (5). Figure 1 shows accurate numerical solution of (5) and (7) for two values of B. In each panel the oscillatory solution appears as a band due to its fast dynamics; the averaged solution varies at a much lower rate and is sufficient to describe the behaviour of the system. When B = 0.52, the motion is essentially confined to the basin of the potential minimum at −1; a small 2% increment in B is sufficient to let the moving particle visit the basin of the minimum at 1 without having to increase the amplitude A of the driving forcing, thus demonstrating the existence of vibrational resonance. Note that, once the stationary regime is attained, the averaged solution is periodic with period 2π/v = 20π and, as discussed above, does not have the Y ↔ −Y symmetry.

Fig. 1

As pointed out before, the solution Y of (7) approximates the true solution z at stroboscopic times. If approximations to z at non-stroboscopic times are also of interest, they may be easily obtained (without additional numerical integrations) by applying to the numerically computed Y the change of variables (3) (see (4)). If the change is truncated to exclude words with three or more letters, for the runs in Figure 1 the maximum on 0 ≤ t ≤ 400 of the magnitude of the discrepancy between the true z and the approximation obtained in this way is 0.080 for B = 0.52 and 0.481 for B = 0.53.

Additional corrections of sizes 𝒪(1/ω²), 𝒪(1/ω³), . . . may be added to (7) by considering words with three, four, . . . letters. Explicit formulas for the 𝒪(1/ω²) terms (three-letter words) are given in [15]. The number of words to be considered increases exponentially with the length of the word. The second column of Table 1 gives, for each n = 1,2,...,7, the number of words with n letters for which the associated basis function f_w is not identically zero. By using a computer algebra programme we have found the corrections of size 𝒪(1/ω²), . . . , 𝒪(1/ω⁶). The corresponding averaged systems are not reproduced here for obvious reasons.

Table 1

Higher-order averaging

For the parameter values used in Figure 1, we have integrated numerically the averaged system based on words with ≤ n letters n = 1,2,...,7. Also listed in Table 1 is the maximum error over 0 ≤ t ≤ 400 when z is approximated by applying the change of variables including words of 1, 2, . . . , n letters to the numerical solution of the averaged system employing the same words. It is apparent that the higher-order averaged systems, which, as emphasized before, may be easily integrated numerically, provide accurate approximations to the true oscillatory solution.