Introduction

The Boussinesq equation is a nonlinear partial differential equation of fourth order, as follows:

This equation was formulated as part of an analysis of long waves in shallow water. It was subsequently applied to problems in the percolation of water in porous subsurface strata. It also crops up in the analysis of many other physical processes. Many different closed-form, series approximation, and numerical solutions are known for particular sets of boundary and initial conditions. The α,β,γ elements are real-valued parameters that are physically meaningful. Top

In the following, we shall be applying commutative hypercomplex mathematics. This is a system of mathematics that obeys the axioms of the classical complex variables and behaves in all ways like the classical complex analysis, while treating a 4-D independent variable. In order to understand the following, the reader should first review the Hypercomplex Math page.

We shall present something new in the form of a "near-general" analytical solution for the Boussinesq equation (we shall explain) in terms of commutative hypercomplex mathematics. We shall show that, unlike for linear partial differential equations, there is not a variety of eigenfunctions from which a solution may be constructed in a variety of ways, but only one characteristic function that satisfies the equation when analyticity is required. The solution will be presented as an analytic function of one 4-D variable and as a complementary pair of functions of a classical complex variable. At the end, we shall reduce all expressions to familiar real functions that are verifiable solutions of the Boussinesq equation. Top

We emphasize analytical solution because the solution that we shall develop is analytic in the classical complex variable sense (i.e., is continuous and single-valued in the region of interest) and is "conservative" in physics terminology. Other solutions might be available that do not preserve, for example, total energy or mass in the full four-space. The solution is general in that we shall first transform the PDE into an ODE and then perform four integrations to obtain a solution with four arbitrary constants of integration. It is also general in the sense that it is a 4-D function that implicitly contains several different solution forms in 2-D and in one dimension and time, as we shall show.

One might recall that partial differential equations, especially those used to model the behavior of some material substance, typically describe the behavior of some variable, parameter, or physical effect about a point. They are typically derived from a force-energy-momentum-mass balance on an infinitesimal element at a point. As an aside, we note that they usually represent approximations only to first or second order of the effects that they describe, in order to keep the equations tractable. We might obtain a closed-form solution for the PDE, only to find that it does not accurately describe the extremes of the physical behavior that we are attempting to model. Top

If we obtain a characteristic function of the Boussinesq equation by whatever means, then the function will embody the same description of physical effects as the PDE and must be viewed and applied in the same way; i.e., as describing the variation of effects about a point, and not necessarily the macro behavior over all space. The Boussinesq equation was developed to describe the immediate, localized reaction of a tiny, incremental element of mass to given external forces and momentum and energy inputs. In physics terms, the integrated result is expressed in body-centered coordinates whose origin moves with the subject particle of mass. At any given instant of time and for given local conditions, the integrated result indicates how the particle of mass will move next within the body-centered frame. We continue to emphasize: The integrated result describes an immediate, localized reaction, and says nothing about the long-term motion of a given particle of mass. For that reason, we must do a finite-element-like numerical calculation in order to coordinate the motions and interactions of all the particles, thereby obtaining a view of the overall motion of the fluid. This view explains why there is not any analytic-function solution of the Boussinesq equation that models turbulent behavior in the large. Any "solution" is point-localized. Top

Hypercomplex Solution

As mentioned earlier, our basic approach to solution is to first convert the (nonlinear) Boussinesq equation to an ODE, then solve it by means of classical methods. To convert partial differentials to ordinary derivatives, we shall use the following consequences of the 4-D Cauchy-Riemann equations: Top

where Z=1x+iy +jz+kct. Refer to the Hypercomplex Math page for a derivation. These relations transform partial differentials into an ordinary derivative with respect to a 4-D independent variable, Z. Observe that the first two equalities on the left are identical to those for 2-D classical complex variables, and the two on the right are extensions for the 4-D case. In the rightmost relation, c is the characteristic speed of propagation in whatever medium we are working. These relations, being direct consequences of the 4-D Cauchy-Riemann conditions, are therefore continuity conditions. By applying them to the Boussinesq equation, we force continuity upon any solution(s). One might argue that we are simultaneously solving the Boussinesq equation with the Cauchy-Riemann equations. However, that would be somewhat misleading because the Cauchy-Riemann equations are satisfied by any analytic function. Top

Making the partial derivative conversions, we obtain:

where Z is a 4-D variable and c is the characteristic speed of propagation of effects. This equation is solvable by simple, direct methods. The first two integrations yield: Top

where A and B are arbitrary 4-D constants of integration. Now we have a problem in how to proceed, because of the AZ term. One approach would be to do a power series expansion of u(Z), put it into the differential equation, above, and solve for the coefficients. The resulting solution form would be useful from a practical, computational standpoint, because we are interested in the effects about a point, in this case the origin. However, if we could continue on and do two more integrations, then we would have the general solution, which would yield the most-complete theoretical insight. Lacking that, the next best thing is a "near-general" solution (three arbitrary constants vs. four for the general case), as follows. The constant A is arbitrary, so we are at liberty to set it to zero, leaving: Top

Now, we can multiply both sides of the equation by the first derivative of u and integrate, to obtain: Top

The element C is a third constant of integration (recall that we zeroed and discarded the A constant).

As an aside, we note that this result easily and explicitly shows that the KdV equation is a special case of the Boussinesq equation. In a previous work, (KdV page), we solved the KdV equation using the same methodology as here, and at this point in the analysis we obtained qualitatively the same first-order, nonlinear differential equation. The only substantive difference is that in the KdV equation the u² term has a constant coefficient, kc, which is just a special case of the parameterized form that we have here. In fact, it could be argued that the second-order, ODE form of the Boussinesq equation is just a KdV equation with forcing terms, the latter being the terms that were eliminated when we zeroed out the AZ term, above. Top

[The following arguments might appear convoluted, but bear with me. We shall eventually arrive at real-valued functions that are verifiable solutions of the Boussinesq equation.]

In the KdV analysis, we determined that an equation of the form that we have here can be integrated if we can factor the polynomial in u. W. F. Ames gives the solution in terms of the factors, at the same time establishing the functional form that the solution will have, here [W. F. Ames, Nonlinear ordinary differential equations in transport processes (Academic Press, New York, 1968), p. 55]. From simple theory of polynomials, there will be three factors and they will all be arbitrary, because they will all be functions of the arbitrary constants B,C, the fixed parameters α,β,γ, and simple constants. Let us denote the new, arbitrary factors by ρ,σ,ω. Here, the coefficients of the Boussinesq equation, hence the ODE polynomial, include parameters α,β,γ that have theoretical significance. We need to preserve them in the final solution, insofar as possible. Accordingly, we rearrange the first-order ODE as follows: Top

Clearly, we could solve for the ρ,σ,ω factors in terms of the coefficients of the corresponding polynomial. However, each would be a function of the arbitrary constants B,C, thus would themselves be arbitrary. We therefore simply adopt the ρ,σ,ω factors as surrogates for the arbitrary constants of integration. In this approach, the alpha parameter loses its separate identity and is merely folded into the three assumed arbitrary factors. That is not unreasonable, because the alpha parameter is the coefficient of the u_xx term, which is not the most important in giving the Boussinesq equation its characteristic behavior. Top

If we expand the factored form and compare corresponding terms,

then we see that the u² term in the above equation imposes an unavoidable condition

on the three factors, because all elements on the right are either constants or fixed parameters. This effectively reduces the number of independent factors to two. For convenience, we are free to use ρ,σ as the independent, arbitrary constants; then: Top

This represents no loss of generality, because at this point we have only two arbitrary constants of integration.

Indeed, there is no good reason to explicitly solve for the factors in terms of the (arbitrary) constants of integration. Instead, we can simply use the factored notation and the two independent factor symbols as surrogates for the (so far) two arbitrary constants of integration. The factors are arbitrary, because if we did factor the polynomial, they would be expressed in terms of the constants of integration, themselves arbitrary. With that insight, we can use Ames' solution for the fourth and final integration, yielding: Top

The fourth and final, 4-D constant of integration is D. The function sn[...] is the Jacobi elliptic function, in this case with a 4-D argument, Z. This solution for the ODE-form Boussinesq equation is easily verified in Mathematica, treating Z as the independent variable and treating all constants in the conventional way as regards differentiation. This result is a rather busy-looking formula whose form we can simplify and standardize by defining more conventional-looking parameters, to wit:

All elements on the right-hand sides of these definitions are either constant or arbitrary, so the newly-defined parameters are arbitrary. Note the difference in the new, italicized parameter k and the unitalicized, algebraic basis element k. Secondly, because the constants on the right are arbitrary and can be real or complex, the classical imaginary i can be subsumed into the arbitrary constants, so we do not further display it. With these observations, the solution is: Top

The parameters β,γ are from the original Boussinesq PDE, and the μ,λ,k constants are arbitrary. This is a "near-general" analytical solution for the Boussinesq equation (because it incorporates three arbitrary constants of integration, rather than four for the general case) in terms of the 4-D commutative hypercomplex variable Z. It is a solution of the ODE form of the Boussinesq equation. Though not the most general, it is a characteristic function of the PDE form of the Boussinesq equation that models all of the desired key behaviors, as we shall show. Top

Aside from the parameters β,γ, this is the same solution that we obtained for the KdV equation, which supports our contention that the Boussinesq equation can be viewed as a forced KdV equation. Moreover, we shall show that it models all of the common wave actions. Just from a calculus perspective, one can appreciate that the general solution would only add complication to our above expression, thereby adding complication to all of the basic wave motions that are modeled. We shall attempt, below, to explain the difference from the general solution.

All three constants μ,λ,k in the solution are arbitrary and are derived indirectly from the constants of integration, as indicated. They can be real, classical complex, or 4-D hypercomplex, and can be manipulated to satisfy boundary and/or initial conditions. In all cases, we have a means of interpretation, as we shall show. Top

In the u(Z) expression, the ± sign merely indicates that either positive or negative-going waves are allowed, so for convenience we shall discuss only the positive case from this point forward. Also, the i factor in front shall be assumed to be subsumed into the arbitrary constants, and will not be further shown. Top

We remind the reader that this is not the one and only, final solution for all problems in a Boussinesq medium. Just as does the Boussinesq PDE, it merely represents the behavior about a point. Given an external force, impulse, etc. on a small increment of material at the point, the characteristic function that we have found indicates how the disturbance will be propagated away from that point in the absence of any other disturbances, reflections, interferences, etc. in the medium. Of course, in a real problem all parts of the medium may be simultaneously receiving independent impulses, and the net external impulse on a specific increment of material at a given time is an integration over all the disturbances originating in all of the rest of the material increments, taking account of the fact that the speed of propagation may vary from point to point in the medium. Consequently, a solution for a given problem involving arbitrary boundary and initial conditions must be pieced together point by point within each small time increment, similar to what is done in a finite element method. It might seem that we have gained no advantage. However, the analytical solution can yield further theoretical insight and provide new avenues for numerical solution. Top

The above solution is in terms of one 4-D independent variable, Z. We can break this down into two separate solutions in terms of one classical complex variable, each, by writing it in the canonical form, as follows. Recall the definition of a commutative hypercomplex function or operator: Top

The eigenvalues are classical complex variables, and the above definition means that a differential operator that operates on the 4-D variable Z also operates separately on the two eigenvalues of Z, and that the 4-D solution u(Z) is the sum of two 2-D solutions u() and u():

i.e., the eigenfunctions of u(Z) in the canonical form are also solutions of the original PDE. Top

When we are given a 4-D expression such as u(Z) and we wish to expand it into canonical form, it is not so simple as writing two identical copies with just different independent variables. We must express every component of the u(Z) function, including all constants, in canonical form and then reduce the resulting expression into simplest canonical form, as follows. We first expand the arbitrary 4-D constants:

All of the subscripted components are real or classical complex numbers. The inequalities are required because of the way that the eigenvalues of any 4-D entity are constructed. [Review the definitions of and , above; they can be equal only at the point (0,0,0,0) or if they are of the form (x,y,0,0) or (x,0,0,0)]. Top

Recall that any 4-D commutative hypercomplex entity, such as an analytic function of a 4-D variable, can be represented in 4 X 4 real matrix form. If this were done for our u(Z) solution, then the eigenfunction components and their complex conjugates would be the eigenvalues of the matrix form of u(Z). Again, each of the eigenfunctions of u(Z) is a solution of a Boussinesq equation which is stated in terms of the associated 2-D independent variable or . Top

In the above terms, the eigenfunctions are:

We can carry this a step further and expand u(Z) into 4-D vector form:

This time, the individual vector components of u(Z) definitely are not solutions of a Boussinesq equation, because the Boussinesq equation is nonlinear. Top

Elementary Properties

We have obtained a solution u(Z) in terms of one variable having three space dimensions and time. One of the first questions that must be answered is, "Does it reduce to a solution of the original, one-dimensional Boussinesq PDE when the y,z components are set to zero?" This is easy enough to check. Setting y=z=0, we get: Top

Now, because we have established a classical analysis that is applied to 4-D functions (see the Hypercomplex Math page), we can treat the 4-D constants μ,λ,k and the algebraic basis element k as simple constants with respect to differentiation, much as we would the classical imaginary, i, taking account that k²=1 and k^-1=k. We leave it to the reader to put the above result into Mathematica and verify that it is a solution of the Boussinesq PDE. Top

Our reduced, one-dimensional solution u(x,ct), above, is not immediately illuminating about what we know of the physical situation that it is used to model. We can discern the true behavior when we examine the eigenfunctions u₁(), u₂() with the y,z coordinates set to zero: Top

As before, all of the subscripted constants are arbitrary real or complex quantities, save that complementary pairs cannot be equal unless both are of the form (0,0,0,0), (x,y,0,0), or (x,0,0,0). We can choose meaningful real values for the constants such that we have purely real expressions that are each solutions of the Boussinesq PDE. The sn(x,k) function (below) has the qualitative behavior of the ordinary sin(x) function, except that it has more-rounded peaks and period 4K(k), where K(k) is the complete elliptic integral of the first kind.

Consequently, the sn^-2(x,k) function represents a periodic wavetrain with the following qualitative behavior and period 2K(k):

The periodic spikes go to infinity, which is something that is never seen in real materials, such as in water. In the real eigenfunctions, above, the minimum amplitude is set by the square of the λ parameter. Moreover, the λ parameter modifies the x-scale, therefore the speed of propagation. This interacting behavior is observed, for example, in real water waves, where higher-amplitude waves propagate at greater speeds. Top

Actually, the sn(x,k) function is doubly periodic, and so our solution can model a field of wave spikes above a plane, such as the spiking wave swells seen under certain storm conditions on the ocean. Similarly for the 2-D array of quantum potential wells seen on the surface of a crystal lattice. Top

The infinite amplitude is clearly too idealized for real water waves. We could temper our solutions by doing orthogonal-function approximations of the sn^-2() terms with the same arguments, for example, λ₁ (x-ct), in such a way that the same general shape is maintained, but with finite amplitude and only near-discontinuous transition. Indeed, we could match the shape to what is observed in real wave experiments for specific fluids. With arguments as indicated, the approximation would then realistically exhibit the connection between overall wave amplitude and wave speed. We could then use the tailored solutions in numerical calculations for specific fluids. One such approximation is as follows. We note that sn²= 1-cn² and that 0≤cn²(x,k)≤1; therefore Top

In these terms, our 1-D, sn^-2 (x,k)-form solutions are converted to:

So far, we have made no approximations, but only an expansion. Note that the 1-D, sn^-2(x,k)-form eigenfunctions have periodic singularities, but no single term in the cn^2(n-1)(x,k)-form eigenfunctions exhibits any singularity. It is the sum of the terms that produces periodic singularities, as follows. Each cn^2(n-1)(x,k) term has amplitude unity and is superimposed and centered upon all of the previous terms. As one adds terms, one adds increments of amplitude by straight summation. Therefore, if we truncate the series, we are left with a finite-amplitude approximation of the original 1-D eigenfunctions. The truncated form would take on the appearance of an offset, pinched cn(x,k) function, thusly: Top

This is beginning to look like real water waves, and we can control the amplitude by the number of terms that we include in the truncation. Note, also, that throughout any truncations, we maintain the interaction between wave height and wave speed, such that higher-amplitude waves are propagated at higher speeds.

Next, we must account for the fact that the Boussinesq equation also models solitons (solitary waves). To do so, we can set the arbitrary parameter k to unity and use the identity sn(x,1)=tanh(x) to convert our 1-D eigenfunctions to:

The tanh(x) function is an s-shaped curve with amplitude unity and only one zero crossing, at x=0. Its square has amplitude unity except for one negative-going spike that drops to zero, shown below: Top

As a result, the tanh^-2(x) eigenfunctions represent a complementary pair of positive-going, solitary wave spikes moving in opposite directions. Because of the above, each solitary spike has the following shape:

The offset base for the real tanh^-2(x,t) terms in the eigenfunctions is set by the square of the λ parameter. The amplitude is infinite, again pointing out the limitations of the Boussinesq model for ordinary wave action. However, an infinite-amplitude solitary spike can be interpreted as a shock wave moving through the medium. An interesting experiment would be to incrementally manipulate the value of the k parameter in such a way as to show the transition between wave train action and solitary waves. Top

Similarly to what we did with the sn^-2(x,k)-form eigenfunctions, above, we can do an expansion and truncation of the tanh^-2(x,t) forms. We use the identity tanh²(x)=1-sech²(x) and the fact that 0≤sech²(x)≤1 to assert:

In these terms, our 1-D, tanh-form solutions are converted to: Top

Exactly as above, the tanh^-2(x,t) forms exhibit singularities, but no single term in the sech(x,t) expansion does so. As above, each sech^2(n-1)(x,t) term is superimposed and centered on all of the preceding ones, and each adds an increment of amplitude, eventually building to infinity. If we truncate the expansion after the sech⁴(x,t) term, for example, we can obtain a finite amplitude and a very realistic wave form that is often seen on the ocean: Top

Lastly, many published Boussinesq solutions for particular problem conditions are based around the sech²(x) function. This function is especially apropos for shallow-water, solitary wave behavior because it has a solitary hump shape and range 0≤sech²(x)≤1. Here, if we truncate after the first term, we are left with:

These forms have the same qualitative appearance as the sech⁴(x,t) truncation, above, but with a lower amplitude:

The relative height and width are adjustable via the λ parameter. For example: Top

These forms can represent low-hump waves (solitons) moving in opposing directions. The approximation is good because, although all the terms in the expansion have amplitude unity, we are raising 0≤sech²(x)≤1 fractional numbers to succeedingly higher powers, forcing them closer to zero. Therefore, each higher-order term represents a centered, superimposed, increasingly-narrower peak with amplitude unity. In a real water wave, their effect would be to increase the height of the wave and narrow its peak, trying to drive the amplitude to infinity. They even hint at the relative fluid motion within a wave as it is being formed. The shape approximation is best in the leading and trailing edges of the peak, because the higher-order terms have much narrower skirts than the sech²(x,t) term. It may well be that the higher-order terms are damped out in water, by processes and effects that are amenable to experiment. Top

We emphasize that all of the above idealized, propagating disturbances can occur only under well-prescribed circumstances. Basically, a single impulse or equally-spaced wavetrain impulses must propagate into a uniform, stationary medium. There must be no interferences, variable resistances, reflections, etc. If, on the other hand, there are multiple, random impulses at different points in the fluid, for example, then we must use a numerical, finite-element method of solution, using the characteristic functions given above in place of the original PDE. Top

As shown above, we can control wave height by truncating the series expansion. Further, we can even control the shape of the truncated form, as follows. We can run wave generation experiments measuring height, shape, and speed, then fit an approximating function to the data of the form:

In this approximation, the δ_i coefficients represent both a shaping and damping effect. For damping, they might decrease in magnitude as some function of the power of the succeeding sech(x,t) terms. For shaping, they might increase or decrease. They could be determined via a least-squares fit to experimental data. Once these are determined, one should then have a very accurate, specific model of the dynamics of a given fluid. These effects will vary from fluid to fluid, due to differing viscosity and other fluid properties. This type of experimental approximation can also be carried out on the
sn^-2(x,k) form of the eigenfunctions. Top

Implications of Forcing Term

We are now in a position to give some qualitative account of the effects of the forcing term AZ in the equation

_{Eqn. (A)}

What we say must be speculative because we do not have the general solution in hand. Recall that we did not take the AZ term into account when we obtained our "near-general" solution. However, as a special case, we have seen that we can have a u(x,ct) solution that is proportional to a tanh^-2[-λ(x-ct)] form moving to the right along the +x-axis. This is a solitary spike with a positive amplitude and, in the real world, it can have a rounded top. Let α,β,γ, A,B be ordinary, positive, real constants; then at any given time we must have

_{1-D Eqn. (A)}

All terms except the second derivative are specified positive, so the second derivative must be negative and of such magnitude as to balance all of the other terms. Consider the apex point of a rounded hump moving to the right. The farther right it goes, the larger the Ax forcing term becomes, hence the greater negative the second derivative term at the apex must be to balance it. The effect is to sharpen and narrow the wave. Given what we know about real fluids, conservation of energy and momentum dictates that as the wave narrows, its amplitude must go up. Therefore, the medium is resisting the motion of the disturbance and the wave "piles up" against the resistance. A moving shock discontinuity develops, regardless of the original shape of the disturbance. In a real fluid governed by the Boussinesq equation, a surface wave would narrow and increase in height with distance, eventually jetting up to an unsustainable height. A pressure wave within the fluid would eventually narrow into a discontinuous shock front. Top

The above scenario, while plausible, is not probable. Real waves tend to broaden and dissipate with distance. By arguments similar to the above, we can postulate a dissipative case. Let A be real and negative. With distance in the +x-direction, an increasingly-negative Ax quantity must be compensated for in the second derivative term, which itself must be negative at top in order to have a positive-going peak. Therefore, if Eqn.(A) is to balance, then the second derivative must become less negative with distance. The wave broadens in shape and, by conservation laws, loses amplitude. In either case that we have discussed, the effect is that wave height and half width vary with distance, which is an effect that is seen in the real world.

Conclusions

We have found a characteristic function for the Boussinesq equation that exhibits the key behaviors of wave action. We have integrated four times and have obtained three independent, arbitrary constants of integration. In one relatively-general function, the result implicitly models traveling wave trains, shock waves, and solitons, all depending upon the selection of values for the three arbitrary constants. In all cases, the solution reduces to a pair of complementary solutions with wave actions moving in opposing directions, either of which might be blocked by the initiating boundary. It shows that the Boussinesq equation overstates the possible and allowable real water wave height, but we have suggested a way to incorporate wave damping by fitting the characteristic function to real data.

We have discovered that the ODE form of the Boussinesq equation can be viewed as a forced KdV equation. We have solved the unforced part, and have obtained a KdV-like solution. The completely general solution will be more difficult to obtain in closed form, but should be readily amenable to power series expansion or numerical treatment, starting from a base of the twice-integrated (second-order) ODE form of the equation wherein the AZ forcing term first appears. We have argued that the forcing term can be used to account for dissipative behavior with distance. Top

Addendum (3/2/04): Obtaining a Boussinesq equation from the Korteweg-de Vries equation.

By use of analytic function theory, it is simple to derive a Boussinesq equation from the KdV equation. We start with the KdV equation in the following form:

Assuming that we would want an analytic solution (continuous and single-valued - why wouldn't we?), then we can use the Cauchy-Riemann conditions to assert:

Next, we take another partial derivative with respect to x to obtain:

This is looking Boussinesq-like, but we need a term. By use of the second-derivative relation from above, we can split the leftmost term into two parts by use of a fractional parameter 0 < < 1 :

Now, it is just a matter of dividing through by the coefficient of the first term and relabeling the remaining coefficients to obtain a Boussinesq equation, Top

The only thing out of the ordinary is that now two of the α,β,γ parameters implicitly include the basis element k (making them on longer strictly real quantities). However, this is no problem, because we have already seen that the solution of the Boussinesq equation includes 4-D constants.

Free hit counters

© 2003, 2006
Clyde M. Davenport
cmdaven@comcast.net