Residue-regulating homotopy method for strongly nonlinear oscillators

It is highly desired yet challenging to obtain analytical approximate solutions to strongly nonlinear oscillators accurately and efficiently. Here we propose a new approach, which combines the homotopy concept with a “residue-regulating” technique to construct a continuous homotopy from an initial guess solution to a high-accuracy analytical approximation of the nonlinear problem, namely the residue-regulating homotopy method (RRHM). In our method, the analytical expression of each order homotopy-series solution is associated with a set of base functions which are preselected or generated during the previous order of approximations, while the corresponding coefficients are solved from deformation equations specified by the nonlinear equation itself and the auxiliary residue functions. The convergence region, rate, and final accuracy of the homotopy are controlled by a residue-regulating vector and an expansion threshold. General procedures of implementing RRHM are demonstrated using the Duffing and Van der Pol–Duffing oscillators, where approximate solutions containing abundant frequency components are successfully obtained, yielding significantly better convergence rate and performance stability compared to the other conventional homotopy-based methods.

Abstract It is highly desired yet challenging to obtain analytical approximate solutions to strongly nonlinear oscillators accurately and efficiently. Here we propose a new approach, which combines the homotopy concept with a ''residue-regulating'' technique to construct a continuous homotopy from an initial guess solution to a high-accuracy analytical approximation of the nonlinear problem, namely the residue-regulating homotopy method (RRHM). In our method, the analytical expression of each order homotopy-series solution is associated with a set of base functions which are preselected or generated during the previous order of approximations, while the corresponding coefficients are solved from deformation equations specified by the nonlinear equation itself and the auxiliary residue functions. The convergence region, rate, and final accuracy of the homotopy are controlled by a residue-regulating vector and an expansion threshold. General procedures of implementing RRHM are demonstrated using the Duffing and Van der Pol-Duffing oscillators, where approximate solutions containing abundant frequency components are successfully obtained, yielding significantly better convergence rate and performance stability compared to the other conventional homotopy-based methods.
Keywords Residue-regulating Á Analytical asymptotic method Á Auxiliary residue functions Á Strongly nonlinear oscillator The coefficient of kth base function appeared in the mth-order residual computing error

Introduction
Solving the equations of motion of strongly nonlinear oscillators is an important subfield in nonlinear dynamics. Perturbation [1] is the most commonly used technique in this field, but it is applicable only when small or large parameters are included in the equation, which severely limits its universality. To overcome this problem, the homotopy analysis/ asymptotic method (HAM) and a class of its modifications have been developed to obtain approximate analytical solutions of nonlinear equations [2][3][4][5][6][7][8].
These methods use the concept of homotopy to deconstruct the nonlinear equation into a family of sub-equations that can be easily solved, thus lifting the restriction of the necessity for inherent or artificial small (or large) parameters. Through the HAM, an approximate series solution expressed by a set of preselected base functions can be sequentially solved from a family of so-called deformation equations associated with a homotopy of equations. And it mathematically includes a variety of commonly used non-perturbation methods, such as Adomian decomposition method [9], Lyapunov's artificial small parameter method [10], the d-expansion method [11], and the variational iteration method [12]. In brief, for a given nonlinear equation N /ðxÞ ½ ¼0, one can construct a homotopy of equations (zeroth-order deformation equation) ð1Þ where x denotes the independent variables, L is an auxiliary linear operator, q 2 0; 1 ½ is an embedding parameter, / x; q ð Þis the unknown function, u 0 x ð Þ is an initial guess solution, h is an auxiliary parameter, and HðxÞ is an auxiliary function. Expand the functions involved in (1) into Maclaurin series with respect to q, then differentiate the zeroth-order deformation equation m times with respect to q, divide it by m!, and set q ¼ 0, the mth-order deformation equation can be obtained as Through this procedure, a series of functions u m x ð Þ can be solved from the deformation equations order by order, and must be a solution of the original nonlinear equation if it converges as q ! 1. The HAM provides us a homotopy with much better practicality compared to the homotopy continuation method [13], because the convergence region and rate of the solved series solution can now be regulated by the auxiliary parameter h [14]. That is also the reason why h is usually regarded as the convergence-control parameter in HAM. Additionally, the convergence performance of the HAM can be further improved through a number of optimization techniques with a more generalized form of homotopy [2,[15][16][17] where h 2 , H 2 ðxÞ, and G / x; q ð Þ; q ½ are additional auxiliary parameter, function, and operator, respectively, A q ð Þ and B q ð Þ are embedding analytical functions. In theory, it is free to choose the abovementioned embedding functions, auxiliary parameters, functions, and operators to construct a good enough homotopy to obtain a convergent series solution with the HAM as long as it satisfies the basic rules of solution expression, coefficient ergodicity and solution existence [2,3]. Most importantly, it provides a convenient way to adjust the convergence of the homotopy. Therefore, the HAM exhibit reasonable applicability to both strong and weak nonlinearities and has been successfully applied in various fields of science and engineering involving nonlinearity [2,3,[18][19][20][21][22][23][24][25][26].
However, there can be infinite auxiliary linear operators that satisfy the three basic rules for a given nonlinear equation and a certain set of base functions when implementing above-mentioned homotopy methods, but yet there are no mathematical theorems which can guide us to choose a proper linear operator L that can guarantee a convergent and efficient homotopy. To be specific, when solving nonlinear oscillator problems with the HAM, it is straightforward to choose time-related trigonometric functions as base functions to express the steady-state solution [27,28] or exponential and trigonometric functions for transient state [29]. But various single-frequency or multi-frequency auxiliary linear operators with or without damping terms can all meet the basic requirements for constructing a homotopy [27][28][29][30][31][32], and there are no universally applicable rules on how to choose the effective ones. Currently, a proper choice on L requires not only sufficient mathematical knowledge, but also extensive experience in applying the HAM, especially for problems with strong and hybrid nonlinearities. What is worse, even for some well-investigated nonlinear problems with a commonly used choice of linear operator L but simply with slightly modified parameters (just like the examples later shown in this paper), the convergence will be very likely to fail and the proper choice on L becomes elusive. In short, given a kind of nonlinear problem, it is a known challenge to choose a universal auxiliary linear operator so that convergent solution can be obtained within a wide range of parameters, which greatly impairs the deployment of the HAM, particularly for complicated nonlinearity.
Different from the HAM, there is a class of methods that do not use auxiliary operators to construct a homotopy of equations but only apply the concept of homotopy to expand the solution into a power series with respect to q thus turning the nonlinear equation into a series of solvable sub-equations, such as the iterative homotopy harmonic balancing approach (IHHBA) [33,34], global residue harmonic balance method (GRHBM) [35,36], spreading residue harmonic balance method (SRHBM) [37], and forward residue harmonic balance method (FRHBM) [38].
Although these methods are easy to implement and do not involve any auxiliary operators, they all lack overall theoretical guidance and convergence-control techniques, thus not suitable for solving strongly nonlinear problems.
In this paper, we present a new asymptotic method based on the concept of homotopy for solving strongly nonlinear oscillators without the participation of any preselected additional auxiliary mathematical operators. This method is regarded below as the residueregulating homotopy method (RRHM), which initiates with a simple guess approximation formed by just a few base functions and can gradually reduce the residue with controllable convergence rate and accuracy. The RRHM is also capable of dealing with strong nonlinearities since it requires no participation of any small or large perturbation parameter and includes a convergence-control technique. It not only greatly facilitates the implementation of homotopy for solving strong nonlinearity by avoiding the difficulty in choosing a proper auxiliary linear operator but also demonstrates better convergence rate and performance stability than the traditional HAM and other above-mentioned methods. The basic ideas of the RRHM are introduced in Sect. 2. Then, examples of strongly nonlinear oscillators involving Duffing and Van der Pol-Duffing nonlinearities and the discussion are presented in Sect. 3 and Sect. 4 to demonstrate the advantages of our method and further illustrate the operational process in detail. Finally, the concluding remarks are presented in Sect. 5.
where S / x; q ð Þ; x; q ½ is an auxiliary residue-regulating operator, which satisfies and can be expanded into Maclaurin series with respect to the embedding parameter q as Then substitute (4) and (8) into (6), differentiate it m times with respect to q and divide it by m!, and finally set q ¼ 0, the mth-order deformation equation of the RRHM reads where R mþ1 u m ; x ð Þ are defined as (3). In order to distinguish it from the zeroth-order deformation equation when m ¼ 0, this equation will be referred to as the initial-order deformation equation hereafter. To solve the series solution from this family of deformation equations, we first choose a proper universal set of base functions B for the expression of the whole homotopy-series. Then, for each order of approximation, we choose a subset of B where N m is the number of base functions included in the subset, then correspondingly let where a ðu m Þ k are unknown parameters. Substitute (11) into the mth-order deformation equation, we have where e ð6 2U m Þ n x ð Þ are the base functions which belong to B but not to U m . By setting the coefficients of base functions belong to U m as zero, one can attain N m algebraic governing equations with N m unknown parameters from which all the unknown parameters in (11) can be solved. Sequentially, a homotopy-series / x; q ð Þ defined as (4) can be obtained and must satisfy the nonlinear equation considered if it converges when q ! 1. Note that the value of N 0 should not be too large if one wishes to obtain an analytical approximate solution, because the initial-order governing equations are nonlinear. But the higher order governing equations are all linear, so that the numbers of base functions which can be further taken into consideration are not limited. Similar to the traditional HAM, as the embedding parameter q increases from 0 to 1, / x; q ð Þ varies continuously from the initial approximation u 0 x ð Þ to at least one of the exact solutions of the considered nonlinear equation. However, the deformation equations in the RRHM do not symbolize the continuous deformation from a linear equation to the considered nonlinear equation, but a homotopy where S / x; q ð Þ; x; q ½ varies from the initial residual error to zero. In later content, we show that a better homotopy can be constructed without involving any auxiliary linear operator through the residue-regulating process. And the RRHM is also of great flexibility on account of the freedom of choosing auxiliary residue functions S m m ! 0 ð Þ and solution subsets U m m ! 0 ð Þ.

Construction of solution subsets
In the RRHM, it is totally free to select the base functions from B to determine the analytical expression and governing equations of the solution to each order deformation equation. However, in order to calculate higher order approximations with high efficiency, a rule that can be automatically executed for the construction of solution subsets U m m ! 1 ð Þ must be established. Here, we take the residual error of each order approximation as the standard to select base functions for U m m ! 1 ð Þ. Substitute the ðm À 1Þ thorder approximation / mÀ1 ¼ P mÀ1 n¼0 u n x ð Þ into the original nonlinear equation, we have the ðm À 1Þ thorder residual error expressed by the base functions Then put the base functions whose absolute values of the coefficients are sufficiently large together with those in U mÀ1 into U m : where w is an artificially chosen expansion threshold. By changing the value of w, one can modulate the numbers of base functions contained in the solution subsets and thus control the computing efficiency and accuracy of the RRHM.

Construction of the auxiliary residue functions
To construct a series of auxiliary residue functions S m m ¼ 0; 1; 2; 3; Á Á Á ð Þ that facilitate the convergence of the homotopy, one must understand primarily the role it plays in the process of order-by-order approximation and restrictions it should satisfy as a part of the homotopy. First of all, because u m can only nullify the coefficients of the base functions within U m in the mth-order deformation equation, the remaining terms expressed by other base functions must be nullified by S m to make (9) true. Secondly, as it is, these remaining terms can be considered as the mth-order residual computing error of the homotopy and might be obliterated in higher order deformation equations, because the higher order solution subsets will expand as the order increases until the residual error is small enough according to (14) and (15) and more base functions can be included into the residue-regulating process. Therefore, the residual computing error should be included in the higher order deformation equations, that is, in S n n [ m ð Þin some way. Finally, (7) must be true so that the convergent series is the solution of the considered equation. On these grounds, we put in a straightforward strategy to construct the auxiliary residue functions, that each S m m ! 0 ð Þ should consist of two parts, one of which contains the mth-order residual computing error (denoted as S ðqÞ m ) and the other is related to all lower order residual computing errors (denoted as S ðpÞ m ). Let and assume then it is straightforward to define S ðqÞ m as and Additionally, in order to reflect the idea of gradually eliminating the accumulated residual computing where v m m ! 1 ð Þare freely chosen residue-regulating parameters. We then have a recurrence formula Therefore, It can be seen that the expression of S is not prescribed in advance, but obtained in the order-by-order calculation. In order to make the series solution convergent, we must choose a proper residue-regulating vector Here, we introduce a plain form of the vector: where v c is called the uniform residue-regulating parameter used for adjusting the convergence of homotopy. And if the calculation yields  (16) and (20). Then S ðqÞ m is formed by all remaining terms as described in (17) and (18) after obtaining the solution u m . When solving nonlinear equations through this procedure, except the base functions and the initial solution subset, only two parameters need to be chosen artificially. Although the approach of constructing the auxiliary residue functions that satisfies the above conditions is not unique in principle, we will show that the convergent series solution of nonlinear equations considered can already be obtained effectively by such a plain approach with the RRHM due to the residue-regulating vector v, especially when combined with some convergence enhancement strategies. First, let us consider a nondimensional Duffing oscillator [39] with dual-frequency excitations: where c and k are the linear viscous damping coefficient and nonlinear stiffness coefficient, respectively, and the dot denotes differentiation with respect to t. F i cos X i t ð Þ i ¼ 1; 2 ð Þare two excitations with different frequencies and complex amplitudes. Substitute (26) into (9), then follow the standard procedures, we have the mth-order deformation equation where d is the Kronecker delta. To solve its periodic solution, we choose the periodic trigonometric functions as the base functions and denote them in a simplified manner for the convenience of expression: Note that the sine and cosine terms of the same frequency always appear together in the problems of nonlinear dynamics, thus this notation will not cause ambiguity. We simply choose as the initial subset of solution and accordingly form an initial approximation where a ðu 0 Þ 1 and b ðu 0 Þ 1 are unknown parameters. Substitute the initial approximation into the first-order deformation equation, we have 3 a are solved from the governing equations obtained by setting the coefficients of base functions belong to U 0 (cos X 1 t ð Þ and sin X 1 t ð Þ) as zero [31].
Meanwhile, the expression of S which is the opposite to the summation of all terms formed by base functions not involved in U 0 as described by (18). In addition, we have S according to (20) and (24). Substitute u 0 into the original nonlinear equation, we have the zeroth-order residual error expressed by the base functions: where x m;n all belong to B. The zeroth-order residual error equals to ÀS ðqÞ 0 , which is not necessarily true for higher order approximations. In this example, we use the frequency components with sufficiently large amplitudes in the residual error to form the solution subsets, that is where w is artificially chosen threshold. For instance, if we set w ¼ 0, then we have Y 0 :¼ X 2 ; 3X 1 f gand U 1 :¼ X 1 ; X 2 ; 3X 1 f g . By repeating this procedure, we can continue to calculate the higher order approximations and attain a convergent series solution of the nonlinear equation if the chosen base functions, initial solution subset, and values of w and v c are appropriate.
Consider the nonlinear equation (26) with c ¼ 0:2, k ¼ 1, F 1 ¼ F 2 ¼ 0:25, X 1 ¼ 1:0, and X 2 ¼ 1:5, both the traditional HAM and RRHM are employed to obtain its series solution expressed by the base functions. For the traditional HAM, we use a plain form homotopy as (1), where L f ½ ¼ o 2 f =ot 2 þ X 2 1 f , which is a common choice in the literature [2,27,31], and the initial guess approximation takes the form as (30). To improve the computational efficiency, we use the discrete squared residual error E m as a substitution for the analytical squared residual error to evaluate the accuracy of the approximations, which is where Ds ¼ T=N, T is the minimum period of the periodic homotopy-series, N is an integer and set to 2000 in this work. The valid region of h is determined by plotting E m versus h at different order of approximation m ¼ 5, 10, 15 and 20 as shown in Fig. 1a, which indicates that the optimal value of h is around -0.9 for the implementation of traditional HAM. Similarly, we determine the optimal value of v c as -0.8 by the E m À v c plot for the proposed RRHM as presented in Fig. 1b. Using the optimal values and the same initial guess solution u 0 t ð Þ ¼ 0.5589 cos t þ 0.3453 sin t [31], the first 50 order approximations are calculated through both methods. As shown in Fig. 2a, when using the traditional HAM, the residual error decays exponentially throughout the homotopy as the order of approximation increases; as for the homotopy constructed through the proposed RRHM procedure, the residual errors decay at a significantly higher speed and tend to be stable after reaching the limit values related to the threshold w. Then, we construct a new homotopy using the obtained 50th-order solution subset as the first-order solution subset and find that the predetermined first-order solution subset can further increase the convergence rate without changing the limit value of residual error. Moreover, the number of base functions involved in the approximations keeps rising unlimitedly when applying the traditional HAM, but for the RRHM, the expansion of the solution subsets is limited by artificial thresholds as shown in Fig. 2b. The correctness of these homotopy-series is proven by comparing with the numerical result obtained by high-accuracy Runge-Kutta (R-K) method as depicted in Fig. 3. The results show that the RRHM can achieve a higher convergence rate with a smaller number of base functions before reaching the limit value compared to the traditional HAM, and the final accuracy can be controlled conveniently by the threshold w. Note that since there are no more unknowns to be solved other than the coefficients of the base functions, it is needless to set initial conditions to provide additional governing equations when solving the examples presented in our paper with the HAM or RRHM [29,32]. The initial values of approximations will also converge to certain values if the series solution is convergent, and are used as initial conditions for obtaining numerical-forward-integration solutions via the R-K method. Then, we employ the two homotopy methods to solve the same nonlinear oscillator but forced by different excitations as the second example, that is The initial approximation for this example is u 0 t ð Þ ¼ 0.2293 cos t þ 0.3497 sin t. Accordingly, as shown in Fig. 4a, the valid region of h for the 20thorder approximation of the traditional HAM shrinks to a very narrow interval and the optimal value is very close to zero (0.01). On the contrary, the E m À v c plot obtained by the RRHM shown in Fig. 4b indicates that the valid region of v c only shrinks slightly and the optimal value of it is around -0.6. As a result, with the optimal values, through the RRHM, the homotopy still asymptotically approaches to the exact solution of the nonlinear equation at a high convergence rate with considerably fewer base functions involved, while the traditional HAM can hardly reduce the residual error as presented in Fig. 5. We only calculated the first 20 order approximations through the HAM, because the homotopy constructed does not tend to converge at all, and the number of base functions contained therein grow at a very high rate. As for the RRHM, the limit values of residual errors and the numbers of base functions involved are still effectively controlled by the threshold w. The 50th-order approximations obtained by the RRHM with three different values of the threshold w are all in good agreement with the numerical solution, while the approximation obtained by the traditional HAM is significantly different from the others due to its low accuracy as shown in Fig. 6. We attribute the inefficiency of the traditional HAM to the characteristics of the inverse operator of the chosen auxiliary linear operator. It tends to over-amplify the amplitudes of the frequency components close to X 1 , so h must be limited to a range very close to zero, otherwise it will cause divergence. In contrast, the proposed RRHM is significantly less sensitive to the change in excitations. Such an improved performance stability upon traditional HAM clearly demonstrates the advantages of our method in convergence control for complicated nonlinear problems.
Note that we do not exclude the possibility that a better homotopy can be constructed through the traditional HAM by choosing a more suitable auxiliary linear operator L, but it requires much more skills and experience than simply assigning proper values to w and v c in the RRHM. In addition, although one can also limit the base functions involved to a smaller range in advance when implementing the HAM, it is difficult to preselect the major frequency components before obtaining an accurate approximate solution, especially when there are a large number of frequency components involved. But with the proposed method, the major frequency components with relatively large amplitudes can be included during the approximation process according to the values of w as shown in the second example, thereby ensuring the effectiveness of our method in dealing with strong nonlinear problems. We then consider the Duffing oscillator with F 1 cos X 1 t ð Þ ¼ 0:25 cos t ð Þ and F 2 cos X 2 t ð Þ ¼ 0:25 cos 1:1t ð Þ. In this example, the first 200 order approximation is obtained through the RRHM with the v c ¼ À0:2, w ¼ 10 À3 and the same initial approximation as the first example as shown in Fig. 7 (a). It can be seen that the increasing complexity of solution leads to the increasing number of approximations required to achieve the same accuracy. We then use the 5th, 35th, and 65th order approximations from the first homotopy as the new initial approximations to construct the second homotopy with v c ¼ À1 and w ¼ 10 À4 . According to the results presented in Fig. 7b, c and Table 1, it can be seen that the performance of our method is very sensitive to the accuracy of the initial approximation and better initial approximation can considerably increase the efficiency of the homotopy. Therefore, even if the convergence rate can be relatively low due to the increasing complexity of the nonlinear equation or the poor accuracy of the initial approximation, it is still possible to use the iterative approach [40] that is to  6 Solution obtained by the R-K method, the 20th-order approximation obtained by the traditional HAM and 50th-order approximations obtained by the RRHM with w ¼ 10 À3 , w ¼ 10 À4 and 10 À5 (F 1 cos X 1 t ð Þ ¼ 0:1 cos t ð Þ, F 2 cos X 2 t ð Þ ¼ 0:1 cos 1:1t ð Þ) construct a new homotopy with a better initial approximation obtained by the RRHM, to finally achieve a high-accuracy approximation.
As shown in Figs. 3, 6 and 7c, the frequency spectra of three examples all consist of a series of equidistant spectral lines generated by nonlinear frequency mixing of the two harmonic excitations. The minimum spacing between adjacent spectral lines is DX ¼ X 2 À X 1 . When the minimum spacing is sufficiently small compared to the excitation frequencies, the spectrum divides into several clusters around the excitation frequency and their odd frequency multiplications due to the Duffing nonlinearity and the amplitudes of these spectral lines decrease  exponentially from the central to both sides within each cluster. In addition, as the minimum spacing decreases and the excitation amplitude increases, more frequency components need to be considered to achieve the same accuracy.

The Van der Pol-Duffing forced oscillator
In order to demonstrate that the RRHM can also find multiple solutions of a nonlinear equation, we then consider a forced oscillator with Van der Pol-Duffing nonlinearity [27,41]: where l ¼ 0:1 is the damping coefficient, F and X are the amplitude and frequency of the force, respectively. The dot denotes differentiation with respect to t. Using (9) as the zeroth-order deformation equation of the RRHM, we then have the mth-order deformation equation where d is the Kronecker delta. Similar to the examples of Duffing oscillator, we choose the base functions as The initial solution subset is chosen to be U 0 :¼ X f g and the initial guess solution is The artificial threshold w is set to be 10 -20 . First, setting F cos Xt ð Þ ¼ cos 2t ð Þ, three different initial approximations can be obtained by solving the initial-order nonlinear governing equations as shown in Table 2. We then calculate the first 10-order approximations with these initial guess solutions and optimal v c through the RRHM. As shown in Fig. 8a and b, three different convergent series solutions are obtained and are consistent with the results yielded from the traditional HAM as demonstrated in Ref. [27]. Among them, only solution 3 is a stable solution.
Then, the first 20-order approximations of the nonlinear equation (37) with the increased force amplitude of 5, 10, and 15 are calculated through both methods. The corresponding initial approximations are collected in Table 2 and the calculations are conducted with optimal h (for traditional HAM) or v c (for RRHM). As shown in Fig. 9a-c, the 20th-order approximations obtained by both the traditional HAM and RRHM are all consistently with the numerical solutions. For each example, the solution is unique and the amplitude increases with the excitation force, which indicates that there could be a bifurcation point between F ¼ 1 $ 5. Moreover, the RRHM achieves higher convergence rates than the traditional HAM for these nonlinear equations, and the performance of the RRHM is more stable since its convergence rates decay more slowly as the force amplitude increases as presented in Fig. 10.

Discussion
Through the above examples, we have shown the possibility of obtaining convergent homotopy-series solutions of strongly nonlinear oscillators by the proposed RRHM. The Duffing and Van der Pol nonlinearities considered in the examples are two basic and widely studied nonlinearities [39,42]. Oscillators with these two nonlinearities are not only rich in dynamic phenomena [39,[42][43][44], but also frequently used in engineering fields as simplified models to capture the main nonlinear properties of the actual situations [45,46]. Therefore, our method can be directly applied to deal with a large number of related problems, especially for some strongly nonlinear cases where a large number of frequency components are needed to obtain high-accuracy solutions, such as the problems of phononic frequency combs [47,48]. On the other hand, our method can also be potentially extended to solve analytical or semi-analytical solutions of other types of nonlinear equations [23][24][25][26][49][50][51][52][53][54][55][56][57][58] since there are no restrictions are imposed on the form of the considered equations when implementing the RRHM.
In the RRHM, an artificial threshold w has been introduced as a control parameter for the expansion of solution subsets, hence to regulate the analytical expressions of each order approximation, which can no longer be derived directly from the deformation equations due to the absence of auxiliary linear operator. In theory, it is impossible to choose an optimal value for w since it is not explicitly involved within the expression of series solutions. Fortunately, a convergent series solution can be obtained without imposing rigorous restriction on the value of w. In the practical operations of the above examples, as long as the solution subsets do not expand too fast to encumber the computational efficiency and in the meantime continuous reduction of the residual error occurs in the first several orders, we consider the value of w as valid, and so it is actually. In addition, the number of base functions included in the solution will also affect the final accuracy. Qualitatively speaking, in a convergent homotopy, the smaller w is, the more  accurate the convergent series solution can be. We note that it is in fact not necessary to keep w constant throughout the homotopy as shown in above examples, as one can continue to reduce its value after reaching the limit residual error to further increase the accuracy.
As for the residue-regulating vector v, it represents the intensity of the residual computing error expected to be obliterated in each order approximation to some extent.
The condition of 1 þ v m j j¼ 0 (1 [ 1 þ v m j j[ 0) means that the residual computing error accumulated up to the previous order is completely (partially) put into the next order deformation equation.
In above examples, we use the plain form of the vector and successfully regulate the convergence rate through the uniform residue-regulating parameter v c , which is, in effect, the convergence-control vector of the RRHM since it is somewhat similar in function to h of the traditional HAM [14]. We notice some aforementioned methods [33][34][35][36][37][38] exhibit some similarities with the RRHM proposed in this paper in the way of constructing the linear sub-equations. Nevertheless, these works do not include any residueregulating parameters that can control the convergence of the iteration or homotopy. In these methods, the accumulated residue in previous approximations is directly introduced to the equation of next order approximation, which is equivalent to our method with v c ¼ À1. And in this case, the convergence rate is lower than choosing the optimal v c , and sometimes it can even lead to divergence as shown in the first and second example. Additionally, these methods neither construct a general theoretical framework of the homotopy nor propose a general technique to select base functions for each order of approximation, thus resulting in obscurity on the effectiveness of solving complex problems and on the scope of applicability, too.

Conclusion
In this paper, we propose a new method for constructing an efficiently convergent homotopy for solving nonlinear equations, where a family of deformation equations can be constructed without the auxiliary linear operator as required in traditional HAM. The validity of our method has been demonstrated by solving a dual-frequency excited Duffing nonlinear oscillator and a forced Van der Pol-Duffing oscillator with the described approaches of constructing solution subsets and auxiliary residue functions. By properly choosing the uniform residue-regulating parameter v c and threshold for the expansion of solution subsets, w, series solutions with better convergence rates are obtained through our method. In addition, the convergence rate can be further improved by artificially selecting a better subset for the first order approximation or employing the iterative approach. Most importantly, we emphasize that our method is preferable in the sense of ease of use, since it only requires basic mathematical knowledge on the considered nonlinear equation when choosing the base functions and the initial set of solution subset, and then the series solution can be easily obtained following the steps described in this paper without manual intervention. Finally, the proposed RRHM can also be extended to other nonlinear problems due to our generalized formulation, which ensures its broad adaptability to various nonlinearities in the fields of science and engineering. Data availability The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request. Fig. 10 Comparison of the convergence rates in the first 20order approximations through the traditional HAM and the RRHM (F ¼ 5; 10; 15)