1 Introduction

“If the proof is correct, then no other recognition is needed”

by Grigori Y. Perelman

On the homotopy multiple-variable method

and its applications in the interactions

of nonlinear gravity waves

Shi-Jun LIAO

State Key Laboratory of Ocean Engineering

School of Naval Architecture, Ocean and Civil Engineering

Shanghai Jiao Tong University, Shanghai 200240, China

( Email address: sjliao@sjtu.edu.cn )

Abstract The basic ideas of a homotopy-based multiple-variable method is proposed and applied to investigate the nonlinear interactions of periodic traveling waves. Mathematically, this method does not depend upon any small physical parameters at all and thus is more general than the traditional multiple-scale perturbation techniques. Physically, it is found that, for a fully developed wave system, the amplitudes of all wave components are finite even if the wave resonance condition given by Phillips (1960) is exactly satisfied. Besides, it is revealed that there exist multiple resonant waves, and that the amplitudes of resonant wave may be much smaller than those of primary waves so that the resonant waves sometimes contain rather small part of wave energy. Furthermore, a wave resonance condition for arbitrary numbers of traveling waves with large wave amplitudes is given, which logically contains Phillips’ four-wave resonance condition but opens a way to investigate the strongly nonlinear interaction of more than four traveling waves with large amplitudes. This work also illustrates that the homotopy multiple-variable method is helpful to gain solutions with important physical meanings of nonlinear problems, if the multiple variables are properly defined with clear physical meanings.

Key Words nonlinearity, wave resonance, multiple-variable, HAM

1 Introduction

In his pioneering work about nonlinear interaction of four gravity waves in deep water, Phillips [1] gave the criterion condition for wave resonance


where with () is the wave angular frequency, and denotes the wave number. Note that, Phillips resonance condition (1) works only for weakly nonlinear waves with small amplitudes, because is the angular frequency of linear theory for a single gravity wave with small amplitude.

Especially, in case of , Phillips [1] showed that, if and , then a steady-state solution for the triad did not exist, and the amplitude of the third wave, if initially zero, would grow linearly in time. This conclusion was confirmed by Longuet-Higgins [2] via perturbation theory and supported by some experiments [3, 4]. Besides, Benney [5] solved the equations governing the time dependence of the resonant modes and studied the energy-sharing mechanism involved.

Although half century passed since Phillips’ pioneering work [1], there exist still some open questions about nonlinear interaction between gravity waves. First, Bretherton [6] pointed out that the perturbation scheme used by Phillips breaks down for large time, and suggested by investigating a one-dimensional dispersive wave model that the amplitude of each wave component should be bounded. Are the amplitudes of each component waves bounded when the resonance condition is exactly satisfied? Besides, based on perturbation method, Phillips [1] gave the condition (1) of wave resonance for only four waves. What is the resonance condition for more waves with large amplitude? It seems difficult to apply the perturbation scheme used by Phillips [1] and Longuet-Higgins [2] to answer these questions, because the related algebra was daunting and “extraordinarily tedious”, as mentioned by Phillips [7].

Perturbation techniques are powerful analytic tools for nonlinear equations, especially when there indeed exists a small physical parameter, i.e. perturbation quantity: the solution is often expressed in a series of the small physical parameter, and the original nonlinear equation is transformed into a sequence of (mostly linear) sub-problems. Frankly speaking, perturbation techniques greatly enrich our knowledge and deepen our understandings about many nonlinear problems. It has a golden times before the times of computer when people calculated mostly by hand and wrote the “analytic” results on one page (or a few pages) of paper. However, some restrictions of perturbation techniques are well-known. First of all, it depends too strongly on the small physical parameters, but unfortunately many nonlinear problems have no such kind of small parameters. Besides, perturbation approximations often break down when the so-called perturbation quantity increases. All of these disadvantages greatly restrict the applications of perturbation techniques.

In history, the calculating tools have a deep impression on the mathematical methods. One example is the development of numerical methods such as the finite element method, the finite difference method, and so on. Fortunately, we are now in the times of computer: a computer can do more than fundamental operations and save a great lot of data in diskette within a seconds! Besides, there are some powerful symbolic computation software such as Maple, Mathematica, MathLab and so on. Today, using these software on a laptop, one can deduce lengthy formulas, calculate a rather complicated expression, and save analytic results on diskette in a seconds which might be hundreds of papers long if printed out. So, when the pencil and paper are replaced by the keyboard, diskette and CPU of a personal computer, the revolution of analytic approximation methods comes, although it is much later than the revolution of numerical approximation methods.

What kind of analytic approximation methods should we have in the times of computer? Can we overcome the restrictions of the traditional perturbation methods? The answer is rather positive and stimulant. Note that, the traditional concept of “analytic” solution came into being in the times of pen and paper, and it was traditionally believed that a “analytic” result must be short and should be written on less than one paper. Today, it needs only a few seconds to save a rather complicated formula in diskette which might be even more than hundreds of papers. Besides, one needs only a few seconds to calculate such a lengthy expression by a personal computer. In fact, one spends much less time to save and calculate such lengthy formula by computer than to calculate and write a simple formula with half-page length by hand. Therefore, in the times of computer, the concept of “analytic” result must be modified, and an “analytic” formula can be very lengthy. Secondly, due to the high performance of symbolic computation of a computer, it should be easy to obtain high-order approximations by means of the new analytic methods. And more importantly, some efficient approaches should be provided to ensure the high accuracy of analytic approximations even for strongly nonlinear problems. In a short, the new analytic method should belong to the times of computer, and besides should overcome the restrictions of traditional methods mentioned above.

Such a kind of analytic method has been developed in 1990s when the early symbolic computation software just appeared: based on the homotopy, a fundamental concept in topology, the homotopy analysis method (HAM) was first proposed by Liao [8, 9, 10, 11, 12, 13, 14, 15, 16]. Different from perturbation methods, the HAM has nothing to do with any physical parameters. More importantly, different from other previous analytic methods, the HAM provides us a simple way to ensure the convergence of solution series. Besides, it has been proved that the HAM logically contains the previous non-perturbation methods such as Lyapunov artificial small parameter method, Adomian decomposition method, the -expansion method and so on. Therefore, the HAM is valid for strongly nonlinear problems and thus is more general. The HAM has been successfully applied to solve different types of nonlinear ODEs and PDEs in science and finance.

The idea of multiple-scales of perturbation methods has clear physical meanings. In this article, based on the homotopy analysis method (HAM), we propose a multiple-variable technique for general nonlinear differential equations, which keeps the clear physical meaning of multiple-scales but abandons completely the small physical parameters of perturbation techniques. In §2, the nonlinear interaction of primary periodic traveling waves in deep water is used as an example to describe the basic ideas of this approach. In §3, we give convergent series solution of a fully developed system of two primary waves even when Phillips’ resonance condition is exactly satisfied. Besides, some interesting results related to multiple resonant waves are reported. Especially, we reveal that Phillips resonance wave condition is mathematically equivalent to the zero eigenvalue of the eigenfunction related to the resonant wave. In §4.1, the expression of the eigenvalue of a fully developed wave system for arbitrary number of traveling waves is derived, and then a resonance condition for arbitrary number of small-amplitude waves is given. In §4.2, a more general resonance condition is further obtained from the physical points of view, which logically contains Phillips’ resonance condition (1) but works for arbitrary number of traveling waves with large amplitudes. In §5, the concluding remarks, open questions and discussions are given. The detailed mathematical derivation is given in Appendix A, and the proof of a convergence theorem is given briefly in Appendix B.

2 Mathematical description

Let denote the vertical co-ordinate, the horizontal co-ordinates, the time, the free surface, respectively. The three axises of are perpendicular to each other, with the unit vector , respectively, i.e. , where is the multiplication dot. Assume that the vorticity is negligible and there exists a potential for the velocity that and


in an incompressible flow, where

On the free surface , the pressure is constant, which gives from Bernoulli’s equation the dynamic boundary condition


where is the acceleration due to gravity. Besides, vanishes following a particle, which gives the kinematic boundary condition


Combining the above two equations gives the boundary condition


where and . On the bottom, it holds


For details, please refer to Phillips [1] and Longuet-Higgins [2] .

Without loss of generality, let us first consider a wave system basically composed of two trains of primary traveling gravity waves in deep water with wave numbers and the corresponding angular frequencies , respectively, where (i.e. the two traveling waves are not collinear). Due to nonlinear interaction, this wave system contains an infinite number of wave components with the corresponding wave number , where are integers. Assume that the wave system is in equilibrium so that each wave amplitude is constant, i.e. independent of the time. Let denote the angles between the positive -axis and the wave number vectors and , respectively, where , i.e. the -axis is perpendicular to the wave numbers . Then,


where and .

Write . According to linear gravity wave theory, the two trains of primary waves traveling with the wave numbers and are given by respectively, where


So, the above two variables have very clear physical meaning: for the considered problem, the wave profile must be a periodic function of and . Mathematically, using these two variables, the time should not appear explicitly for a fully develop wave system. In other words, one can express the potential function and the wave surface for two trains of primary traveling waves. In this way, the mathematical expressions of the unknown potential function and wave surface are clear, with clear physical meanings.

Obviously, it holds

and similarly







where is used, and

In general, it holds


for arbitrary functions and .

Similarly, we have


Then, the governing equation reads


which has the general solution


where are integers and are integral constants.

For the sake of simplicity, define


Using the new variables and , the dynamic boundary condition (3) becomes


On the free surface , the kinematic boundary condition (5) reads




On the bottom, it holds


Given two angular frequencies and , our aim is to find out the corresponding unknown potential function and the unknown free surface , which are governed by the linear partial differential equation (16) subject to two nonlinear boundary conditions (19) and (20) on the unknown free surface , and one linear boundary condition (21) on the bottom. Here, it should be emphasized that, by means of the two independent variables and , the time does not appear explicitly in the unknown potential function and wave surface. More importantly, these new variables have clear physical meanings. This greatly simplifies the problem solving, as shown later in this article.

3 Homotopy-based approach

As mentioned before, the two variables and have clear physical meanings and the solutions of considered problem should be periodic functions of and . From physical points of view, it is clear that the wave surface should be in the form


where is the amplitude of wave component . Note that (17) is the general solution of the governing equation (16). So, the corresponding potential function should be in the form




and is unknown coefficient independent of . Note that the potential function defined by (23) automatically satisfies the governing equation (16). The above expressions are called the solution expressions of and , respectively, which have important role in the frame of the homotopy analysis method.

For simplicity, define a nonlinear operator


where the angular frequencies are given. This nonlinear operator is based on the nonlinear boundary condition (20). Note that it contains a linear operator


Let denote an auxiliary linear differential operator with the property . As mentioned in many articles, one of advantages of the homotopy analysis method is the freedom on the choice of the auxiliary linear operator. Based on the results of linear wave theory, i.e.


we choose such an auxiliary linear operator


Then, let denote the embedding parameter, an auxiliary parameter (called convergence-control parameter), an initial approximation of the potential function with , respectively. We construct such a parameterized family of equations (called the zeroth-order deformation equation) in the embedding parameter :


subject to the two boundary conditions on :





Besides, at the bottom, it holds


Thus, using the property of the auxiliary linear operator (28), we have when the following relationships




which provide us the initial approximations of the potential function and the free surface . Since , when , Eqs. (29) to (32) are equivalent to the original equations (16),(19), (20) and (21), respectively, provided


So, as the embedding parameter increases from 0 to 1, deforms (or varies) continuously from the initial approximation to the unknown potential function , so does from 0 to the unknown wave profile , respectively. Mathematically speaking, Eqs. (29) to (32) define two homotopies:

This is exactly the reason why Eqs. (29) to (32) are called the zeroth-order deformation equations.

Note that we have freedom to choose the value of the convergence-control parameter . Assuming that is so properly chosen that the Taylor series


exist and converge at , then we have due to (35) the homotopy-series solution



are called homotopy-derivatives. In the above formulas, the relationships (33) and (34) are used. At the th-order of approximations, we have

The equations for and can be derived directly from the zeroth-order deformation equations (29) to (32). Substituting the series (36) into the governing equation (29) and the boundary condition (32) at bottom, and equating the like-power of the embedding parameter , we have


subject to the boundary condition at bottom


where . It should be emphasized that the two boundary conditions (30) and (31) are satisfied on the unknown boundary , which itself is now dependent upon the embedding parameter , too. So, it is relatively more complicated to deduce the corresponding equations. Briefly speaking, substituting the series (36) and (37) into the boundary condition (30) and (31) with , then equating the like-power of , we have two linear boundary conditions on :








The detailed derivation of the above equations and the definitions of , , , , , , are given in Appendix A. Note that, the sub-problems for and are not only linear but also decoupled: given and , it is straightforward to get directly, and then is obtained by solving the linear Laplace equation (40) with two linear boundary conditions (41) and (42). Thus, the high-order deformation equations can be easily solved by means of the symbolic computation software.

Liao [8] proved in general that the homotopy-series solutions satisfy the original nonlinear equations as long as they are convergent. Similarly, we have such a theorem:

Convergence Theorem The homotopy-series solution (38) and (39) satisfy the original governing equation (16) and the boundary conditions (19), (20) and (21), provided that


where are defined by (126) and (44), respectively.

A mathematical proof of the above convergence theorem is given briefly in Appendix B. Because the terms and are by-products in solving high-order deformation equations, the above theorem provides us a convenient way to check the convergence and accuracy of the homotopy-series solution. For this reason, we define the residual error squares


for the th-order approximations of and . According to this convergence theorem, the homotopy-series solution (38) and (39) satisfy the original equation and all boundary conditions if and as . Besides, the values of and indicate the accuracy of the th-order approximation of and , respectively.

Note that the auxiliary linear operator (28) has the property


where is defined by (24). Thus, mathematically speaking, has an infinite number of eigenfunctions with the corresponding eigenvalue


In short,

Besides, defined by (24) automatically satisfies the governing equation (16) and the boundary condition (21) at the bottom. Therefore, the inverse operator is defined by


Note that the inverse operator has definition only for non-zero eigenvalue . When , we have


corresponding to the criterion of the so-called “wave resonance” mentioned by Phillips and Longuet-Higgins. Thus, the number of eigenfunctions with zero eigenvalue, denoted by , is the key for solving the problem.

When , it holds

which equals to zero only when (note that corresponds to and thus is not considered here). Similarly, when , the eigenvalue equals to zero only when . So, there exist at least two eigenfunctions and whose eigenvalues are zero, i.e.


for any constants and independent of and . Therefore, it holds in case of two primary waves. As mentioned by Phillips [1] and Longuet-Higgins [2], the criterion (52) of wave resonance can be satisfied for some special wave numbers and angular frequencies. Thus, when (52) is satisfied in case of and , where and are integers with , there exist three eigenfunctions , and

whose eigenvalues are zero, i.e.


for any constants and independent of and . Without loss of generality, Longuet-Higgins [2] discussed a special case and , corresponding to the eigenfunction

In §3 of this article, it implies and when in case of two primary waves, if not explicitly mentioned.

According to the definitions (28) and (45), it holds


which gives the definition of the linear inverse operator


Using this inverse operator, it is easy to solve the linear Laplace equation (40) with two linear boundary conditions (41) and (42), as illustrated below. Here, we emphasize that the above inverse operator has definition only for non-zero eigenvalue