Clyde M. Davenport
cmdaven@comcast.net
© 2003, 2006, 2008
Updated 7/15/06
Introduction  Hypercomplex Solution  Elementary Properties  Forcing Terms 
Conclusions 
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 closedform, series approximation, and numerical solutions are known for particular sets of boundary and initial conditions. The α,β,γ elements are realvalued 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 4D 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 "neargeneral" 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 4D 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 singlevalued 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 fourspace. 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 4D function that implicitly contains several different solution forms in 2D 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 forceenergymomentummass 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 closedform 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 bodycentered 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 bodycentered frame. We continue to emphasize: The integrated result describes an immediate, localized reaction, and says nothing about the longterm motion of a given particle of mass. For that reason, we must do a finiteelementlike 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 analyticfunction solution of the Boussinesq equation that models turbulent behavior in the large. Any "solution" is pointlocalized. 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 4D CauchyRiemann equations: Top where Z=1x+iy+jz+kct. Refer to the Hypercomplex Math page for their origins. These relations transform partial differentials into an ordinary derivative with respect to a 4D independent variable, Z. Observe that the first two equalities on the left are identical to those for 2D classical complex variables, and the two on the right are extensions for the 4D 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 4D CauchyRiemann 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 CauchyRiemann equations. However, that would be somewhat misleading because the CauchyRiemann conditions are satisfied by any analytic function. Top Making the partial derivative conversions, we obtain: where Z is a 4D 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 4D 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 mostcomplete theoretical insight. Lacking that, the next best thing is a "neargeneral" 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 firstorder, nonlinear differential equation. The only substantive difference is that in the KdV equation the u^{2} 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 secondorder, 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 realvalued 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, 1968] gives the solution in terms of the factors, at the same time establishing the functional form that the solution will have, here. 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 firstorder 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 we know from commutative hypercomplex theory that u_{x x}=u_{t t} . Top If we expand the factored form and compare corresponding terms,
then we see that the u^{2} 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, 4D constant of integration is D. The function sn[...] is the Jacobi elliptic function, in this case with a 4D argument, Z. This solution for the ODEform 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 busylooking formula whose form we can simplify and standardize by defining more conventionallooking parameters, to wit: All elements on the righthand sides of these definitions are either constant or arbitrary, so the newlydefined parameters are arbitrary. Note the difference in the new, italicized parameter k and the unitalicized, algebraic basis element k. With these observations, the solution is: Top The parameters β,γ are from the original Boussinesq PDE, and the μ,λ,k constants are arbitrary. This is a "neargeneral" 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 4D 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 4D 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 negativegoing waves are allowed, so for convenience we shall discuss only the positive case from this point forward. 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 of 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 4D 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 4D variable Z also operates separately on the two eigenvalues of Z, and that the 4D solution u(Z) is the sum of two 2D 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 4D 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 4D constants: All of the subscripted components (excluding the ebasis elements) are arbitrary real or classical complex numbers, except that none should be zero. We wish to avoid singularities, later on. Top Recall that any 4D commutative hypercomplex entity, such as an analytic function of a 4D 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 2D independent variable or . Top In the above terms, the eigenfunctions are: We can carry this a step further and expand u(Z) into 4D 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, onedimensional 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 4D functions (see the Hypercomplex Math page), we can treat the 4D 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^{2}=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, onedimensional 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_{1}(), u_{2}() with the y,z coordinates set to zero: Top As before, all of the subscripted constants are arbitrary real or complex quantities. None should be zero, in order to avoid singularities later on. 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 morerounded 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 xscale, therefore the speed of propagation. This interacting behavior is observed, for example, in real water waves, where higheramplitude waves propagate at greater speeds. Top Something very significant has simply and naturally fallen out of the commutative hypercomplex mathematics and the canonical formulation of our result. We have a complementary pair of periodic, travelingwave solutions moving in opposite directions. Of course, either of the pair might be suppressed by the boundary conditions. Apparently, because of the way that the Boussinesq equation is written and in the case of the 1D eigenfunctions, above, we can have only positivegoing spikes, here. However, recall that the mostgeneral, 4D solution u(Z), above, had a ± sign in front of it, which for simple convenience we did not carry forward. Therefore, negativegoing spikes are possible. 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 2D 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 orthogonalfunction approximations of the sn^{2}() terms with the same arguments, for example, λ_{1}(xct), in such a way that the same general shape is maintained, but with finite amplitude and only neardiscontinuous 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^{2}=1cn^{2} and that 0 ≤ cn^{2}(x,k) ≤ 1; therefore we can use the simple geometric progression rule: Top In these terms, our 1D, sn^{2}(x,k)form solutions are converted to: So far, we have made no approximations, but only an expansion. Note that the 1D, sn^{2}(x,k)form eigenfunctions have periodic singularities, but no single term in the cn^{2(n1)}(x,k)form eigenfunctions exhibits any singularity. It is the sum of the terms that produces periodic singularities, as follows. Each cn^{2(n1)}(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 finiteamplitude approximation of the original 1D 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 higheramplitude 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 1D eigenfunctions to: The tanh(x) function is an sshaped curve with amplitude unity and only one zero crossing, at x=0. Its square has amplitude unity except for one negativegoing spike that drops to zero, shown below: Top As a result, the tanh^{2}(x) eigenfunctions represent a complementary pair of positivegoing, 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 infiniteamplitude 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^{2}(x)=1sech^{2}(x) and the fact that 0 ≤ sech^{2}(x) ≤ 1 to assert (by the geometric progression rule): In these terms, our 1D, tanhform 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(n1)}(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^{4}(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^{2}(x) function. This function is especially apropos for shallowwater, solitary wave behavior because it has a solitary hump shape and range 0 ≤ sech^{2}(x) ≤ 1. Here, if we truncate after the first term, we are left with: These forms have the same qualitative appearance as the sech^{4}(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 lowhump 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^{2}(x) ≤ 1 fractional numbers to succeedingly higher powers, forcing them closer to zero. Therefore, each higherorder term represents a centered, superimposed, increasinglynarrower 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 higherorder terms have much narrower skirts than the sech^{2}(x,t) term. It may well be that the higherorder 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 wellprescribed circumstances. Basically, a single impulse or equallyspaced 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, finiteelement 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 leastsquares 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 "neargeneral" 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}[λ(xct)] form moving to the right along the +xaxis. 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 _{1D 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 (i.e., the flow velocity must slow with distance). 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 a perfect 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 +xdirection, an increasinglynegative Ax quantity must be compensated for in the second derivative term, which itself must be negative at top in order to have a positivegoing peak. Therefore, if Eqn.(A) is to balance, then the second derivative must become less negative with distance (i.e., the flow velocity must increase 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 relativelygeneral 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 KdVlike 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 twiceintegrated (secondorder) 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 Kortewegde 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 singlevalued  why wouldn't we?), then we can use the CauchyRiemann conditions to assert: Next, we take another partial derivative with respect to x to obtain: This is looking Boussinesqlike, but we need a term. By use of the secondderivative 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 4D constants. © 2003, 2006, 2008
