New Research In
Physical Sciences
Social Sciences
Featured Portals
Articles by Topic
Biological Sciences
Featured Portals
Articles by Topic
 Agricultural Sciences
 Anthropology
 Applied Biological Sciences
 Biochemistry
 Biophysics and Computational Biology
 Cell Biology
 Developmental Biology
 Ecology
 Environmental Sciences
 Evolution
 Genetics
 Immunology and Inflammation
 Medical Sciences
 Microbiology
 Neuroscience
 Pharmacology
 Physiology
 Plant Biology
 Population Biology
 Psychological and Cognitive Sciences
 Sustainability Science
 Systems Biology
Averaging and renormalization for the Korteveg–deVries–Burgers equation

Contributed by Alexandre J. Chorin, July 2, 2003
Abstract
We consider traveling wave solutions of the Korteveg–deVries–Burgers equation and set up an analogy between the spatial averaging of these traveling waves and realspace renormalization for Hamiltonian systems. The result is an effective equation that reproduces means of the unaveraged, highly oscillatory, solution. The averaging enhances the apparent diffusion, creating an “eddy” (or renormalized) diffusion coefficient; the relation between the eddy diffusion coefficient and the original diffusion coefficient is found numerically to be one of incomplete similarity, setting up an instance of Barenblatt's renormalization group. The results suggest a relation between selfsimilar solutions of differential equations on one hand and renormalization groups and optimal prediction algorithms on the other. An analogy with hydrodynamics is pointed out.
There are many equations of interest in science whose solution is too complicated to be computed with sufficient accuracy, and one may be interested in finding an effective equation, easier to work with, whose solution is an average of the real solution. Recently, there has been a major effort to find rational ways to derive such equations for nonlinear problems (see, e.g., refs. 1 and 2).
Optimal prediction methods (3) have been developed to find effective equations where the average is over a probability density for unresolved degrees of freedom. In the case of systems in thermal equilibrium the derivation of the reduced equations simplifies considerably and turns out to be closely related to renormalization group (RNG) methods (4).
It is of great interest to extend these simpler optimal prediction methods to nonHamiltonian systems, and we attempt to do so here by analogy; a more ambitious extension of RNG ideas to systems far from equilibrium, including diffusive systems, has been presented in refs. 5–7. We work here with the specific example of the Korteveg–deVries–Burgers (KdVB) equation, with boundary conditions that give rise to traveling waves. The problem has a single dimensionless parameter that we call R (to highlight an analogy with the Reynolds number of fluid mechanics). For large values of R the traveling waves are highly oscillatory; we consider a spatial (not ensemble) average of these traveling waves and look for an equation whose solution approximates this average (see Fig. 1). This problem was chosen because of its apparent simplicity, because previous work by Barenblatt et al. (8, 9) has suggested interesting conjectures about its solution and also because of interesting connections with fluid mechanics.
In the following sections I describe the problem, review the equilibrium averaging theory for Hamiltonian systems, set up the analogy for the KdVB equation, and provide a numerical analysis of its validity. I also compare the numerical results to earlier analytical work and point out interesting scaling relations that emerge from the calculations. Note that the averaging/renormalization presented here is applied to stationary solutions; the extension to timedependent problems will be presented elsewhere. A key feature of the work is that the methods are applied to partial differential equations and thus require a careful consideration of scaling.
The KdVB Equation, Its Numerical Solution, and an Asymptotic Analysis
Consider the equation [1] with boundary conditions [2] where the subscripts denote differentiation, x is the spatial variable, t is time, and α ≥ 0, ε^{2}, U ≥ 0 are given constants. The boundary conditions create a traveling wave solution moving to the right (toward +∞) with velocity U/2. We pick as length scale , as time scale T = ε/U^{3/2}, set x = x′/L, t = t′/T, u′ = u/(L/T), and drop the primes; the result is the equation [3] with u_{x}(–∞) = 0, u(+∞) = 0, u(–∞) = 1; is the Reynolds number. For R ≤ 1 the traveling wave has a monotonic profile, whereas for R > 1 the profile is oscillatory, with oscillations whose wave length is of order 1 (10). At zero diffusion (R = ∞) the wave train extends to infinity on the left with a singularity near the location where there would be a shock in the absence of dispersion (11). For finite R the wave train is damped as is the singularity, and the solution of Eq. 3 tends to 1 as x decreases even in the steady state.
The numerical solution of Eq. 3 requires some care. As usual we proceed by fractional steps corresponding to the several terms in Eq. 3. The fractional steps that correspond to the linear terms are performed via a Fourier expansion, with each coefficient advanced in time exactly (the earliest reference I could find for this procedure in this context is ref. 12). For the nonlinear term we use the fourthorder LaxWendrofflike method of Zalesak (13); its accuracy is second order in time and fourth order in space, and so is the accuracy of the overall fractional step scheme, although in practice the error in time is negligible.
Eq. 3 is solved numerically in a domain –X + t/2 ≤ x ≤ X + t/2, in which the traveling wave is stationary. X is picked large enough so that there are no significant oscillations of the solution outside this domain. The boundary conditions are imposed at the edges of this domain rather than at infinity.
The solution of Eq. 3, shifted by –t/2, tends to a steady state. For R ∼20 the convergence is smooth and rapid (see Fig. 2, where the residual variation is caused by the translation); for greater values of R oscillations can persist for a long time, and reaching a steady state at an acceptable cost may require a little averaging in time. At the steady state we average the solution at each point x over the region (x – ℓ/2, x + ℓ/2) and call the result ū, the averaged solution at x; ℓ is a dimensionless multiple of the unit of length L specified above. We are looking for an effective equation g(v, v_{x}, v_{xx},...) = 0 whose solution v approximates this average; this solution can be expected to be smoother than the solution of Eq. 3 and thus require fewer mesh points for an accurate numerical solution.
In refs. 8 and 9, Barenblatt and colleagues have given an asymptotic derivation of an effective equation for this problem. They assumed that R ≫ 1, that the front was averaged over a length much larger than the wave length of the dispersive waves but much smaller than the width of the wave train (i.e., of the region where u is neither the constant 1 nor the constant 0). They also assumed that the dispersion in the effective equation was negligible and thus the effective equation had the form [4] with an effective viscosity (not necessarily constant) ν_{eff} < 1. If one writes u = ū + u′, and the spatial average of u′ is zero, substitution into Eq. 4 yields: [5] An asymptotic analysis then suggested that ν_{eff} was approximately a constant, at least far from the edges of the oscillation region, and a numerical calculation supported this conclusion to some extent. We shall call Eq. 5 the hydrodynamic approximation and the value of ν_{eff} it yields ν_{hydro}.
Renormalization and Averaging for Hamiltonian Systems
To explain what will be done below for the KdVB equation, we quickly review optimal prediction at equilibrium (3), which for a Hamiltonian system is a special case of realspace renormalization (4). Consider a system of ordinary differential equation [6] where φ, x, and R are ndimensional vectors with components φ = (φ_{1},..., φ_{n}), etc., and t is the time. Assume that the system is Hamiltonian: there exists a function H = H(φ) such that for i odd and for i even. The system (Eq. 6) then preserves the canonical probability density Z^{–1}exp(–H/T), where Z is the partition function.
Suppose we want to calculate the average value of a function A = A(φ) with respect to the canonical density. One can do so by Markov chain Monte Carlo based on the Hamiltonian H, or one can solve Eq. 6 in time, average the values of A, and hope for ergodicity.
Divide the component of φ into groups: , , so that , and similarly x = (x̂, x̃), R = (R̂, R̃). Note that R̂ = R̂(φ), i.e., a subset of components of R is in general a function of all of the components of φ. Suppose the function A depends only on the values of the reduced set of m variables , then there is no need to calculate the full statistics of φ. Indeed, one can approximate the system (Eq. 6) by [7] where is the conditional expectation of R̂(φ) given ] is a function of only (3) and is the best approximation of R̂(φ) by a function of in the mean square sense. Hald's theorem (3) states that the reduced system (Eq. 7) is Hamiltonian as well, with a new Hamiltonian , where ; we have folded the temperature T into the Hamiltonian and assume the reader can insert the appropriate factors of T into the equations of motion. The partition function of the reduced system equals that of the original system, , and if the initial data for the reduced system are drawn from the new canonical density Z^{–1}e^{–Ĥ}, then the joint probability density of the variables equals their marginal density in the full system. The average of the function can now be computed by Monte Carlo based on the Hamiltonian Ĥ or by solving Eq. 7 and relying on ergodicity. However, in general, the reduced Hamiltonian system does not produce accurate time evolutions for the components (see ref. 3); it only reproduces the correct equilibrium statistics.
This construction is closely related to renormalization, in the present context, to realspace renormalization (14). Associate the variables φ_{i}, i = 1,..., with points on a regular lattice. Divide the points into groups, each containing the same number q of variables, each group having the same shape and the same number of lattice points (for example, one could divide a 1D lattice into groups of two points). Pick as components of one of the variables in each group (for example, in the case of a 1D lattice divided into groups of two, take the leftmost of each pair). Pick as components of the variables that are not in . Then write and relabel the components of φ^{(1)} so that they are attached to the points in the original lattice. This transformation: H^{(0)} → H^{(1)}, φ^{(0)} → φ^{(1)}, followed by relabeling, is a RNG transformation and a spatial averaging. Note that the Hamiltonian H^{(1)} was obtained from H^{(0)} by averaging the right sides of the equations of motion (see Eq. 7).
With a suitable representation for the Hamiltonians H^{(0)}, H^{(1)}, this operation can be repeated, and one obtains a sequence of Hamiltonians H^{(0)}, H^{(1)}, H^{(2)},.... The fixed points of the transformation H^{(}^{n}^{)} → H^{(}^{n}^{+}^{1)} are the critical points of the system if the system is infinite. The sequence of Hamiltonians can be expected to converge to one of the stable fixed points. If the Hamiltonians have a representation of the form , where the ψ_{i}s are a suitable basis and the expansion has a finite number of terms, then the matrix A with elements is readily computed (14), the eigenvalues of this matrix for H^{(}^{n}^{)} near critical points determine the critical exponents of the system (see any textbook on critical phenomena, e.g., refs. 14 and 15). The corresponding eigenvectors are the “scaling fields” at the critical points, i.e., fields that scale simply near those points.
An Analogy Between the RNG and the Averaging of the KdVB Equation
We draw an analogy between the conditional expectations that define the renormalized variables in the equilibrium Hamiltonian case and an averaging in space that defines “renormalized” variables for solutions of the KdVB equations that are stationary in a moving frame. Averaging over an increasing length scale corresponds either to more renormalization steps or, equivalently, to renormalization with a greater number of variables grouped together. We pick a class of equations in which to seek the “effective” equation, the one whose solutions best approximate the averages of the true solution in the meansquare sense. The choice of meansquare approximation in the KdVB case corresponds to the use of L_{2} norms implied by the use of conditional expectations in the Hamiltonian case, and the choice of a class of equations in which to look for the effective equation is analogous to the choice of a basis for the representation of the Hamiltonian. The calculation of the best coefficients in the chosen class of effective equations corresponds to the evaluation of the coefficients in the series for the renormalized Hamiltonians. Note that in the Hamiltonian case we average the right sides of the equations and in the analogous KdVB case we attempt to average the solutions; this must be so because in the KdVB case we do not have the Hald theorems that guarantee that averaging the right sides produces the correct statistics for the solutions. One expects that the effective equation would have smoother solutions than the original equation and would require fewer mesh points to be properly approximated; in this sense the number of variables is decreased, although we choose not to change the scale of the averaged solution and thus have no analogue of the relabeling of the variables in the Hamiltonian case.
We initially looked for an effective equation in the class of equations of the form [8] where ν ≥ 0, α ≥ 0, β ≥ 0 were constants. This form was suggested by the work of Barenblatt et al. (8, 9). The average solution is where u is a solution of Eq. 2, and ℓ is dimensionless. The problem is to find the values of the parameters in the effective equation that minimize The shift z is needed because the problem is translation invariant, and the numerical procedures can produce a shift that has no intrinsic significance. We found that the last term in Eq. 8 had little effect on the minimum of I and could be safely omitted. The effective equation thus has exactly the form of the original equation but a different value of the coefficient ν (or alternately, a different Reynolds number). The dispersive term can also be omitted (as was done by Barenblatt et al., ref. 8) (see also below), but we do not do so mostly for esthetic reasons.
Dimensional Analysis and Similarity
I briefly remind the reader of the fundamentals of similarity theory (5, 6). Suppose a variable a is a function of variables a_{1}, a_{2},..., a_{m}, b_{1}, b_{2},..., b_{k}, where a_{1},..., a_{m} have independent units whereas the units of b_{1},..., b_{k}, can be formed from the units of a_{1}, a_{2},..., a_{m}. Then there exist dimensionless variables where the α_{i}, α_{ik} are powers, such that Π is a function of the Π_{i}: [9] This is just a consequence of the requirement that a physical relationship be independent of the units of measurement. Suppose the variables Π_{i} are small and assume that the function Φ (about which we know nothing at this stage) has a nonzero finite limit as its arguments tend to zero; then Π ∼ constant, and one finds a power monomial relation between a and the a_{i}. The resulting relation is invariant under a group of scaling transformations generated by changes in the units. A similar argument works if the Π_{i}s are very large. If the function Φ does not have the assumed limit, it may happen that for Π_{1} small, , where the dots denote lowerorder terms, q is a constant, the other arguments of Φ have been omitted and Φ′ has a finite nonzero limit. One can then obtain a power monomial expression for a in terms of the a_{i} and b_{1}, with undetermined powers that must be found by other means. The resulting power relation is an incomplete similarity relation, and the corresponding group of transformations under which the relationship is invariant is Barenblatt's RNG. A relation between Barenblatt's RNG and more standard definitions of the RNG is discussed in ref. 7.
A logarithmic change of variables can transform the problem of determining exponents in incomplete similarity into a wave propagation problem where the wave velocity is unknown (see refs. 5 and 6, and applications in ref. 2). In our current problem the wave velocity is determined in advance.
Numerical Results
In summary, we consider the equation [10] approximated numerically as described above. As before, we define consider an effective equation of the form where v = v(ν, x) satisfies the same boundary conditions as u, and look for ν = ν_{eff}, which minimizes The general situation is shown in Fig. 1, where we have plotted u, ū, for R = 60, and v for ν = ν_{eff} = 4.42; z = –0.4.
In Fig. 3 we plot I(ν) as a function of ν for R = 60, ℓ = 12. I(ν) first decreases fast, then slowly, and then increases. For smaller values of R and ℓ, there is a local minimum at the point where the rapid decrease becomes a slow decrease; since this local minimum disappears for larger values of ℓ we pay no attention to it here. The value ν_{eff} is ν at the minimum of I(ν); the flatness of the curve at that point, especially for large values of R, can make the minimization of I rather painful, and rather than use sophisticated minimization routines we found the minima by tabulation. Note the resemblance of the procedure to system identification methods.
In Figs. 4 and 5 we plot ν_{eff} as a function of R with ℓ = 20, in Fig. 4 in regular coordinates and in Fig. 5 in loglog coordinates. We see that ν_{eff}, the renormalized or eddy viscosity, increases as the “bare” viscosity 1/R decreases. This is plausible because as 1/R decreases the fluctuations in u increase and the portion of the x axis in which there are oscillations increases. From Fig. 5 we conjecture that ν_{eff} is proportional to R^{α} with α = 3/4 (a leastsquare fit gives α = 0.7506). If indeed this is the case, one can obtain an effective equation from this scaling relation even when R is too large for an affordable solution of the original equation.
As we have shown in the preceding section, one has where Φ is an unknown dimensionless function, and we are interested in the case where both R and ℓ are large. The result here is that Φ(R, ℓ) = R^{3/4}Φ(ℓ), an incomplete similarity relation.
Note that in fluid mechanics one also expects ν_{eff}, the eddy viscosity, to increase with the Reynolds number R defined by the bare (i.e., real) viscosity. For example, one can use the scaling laws for the intermediate region of wallbounded flow (16) to derive an eddy viscosity of the form v_{eff} = C_{1} – C_{2}(1/ln R), with C_{1} ≥ 0, C_{2} ≥ 0. However, the increase is much slower. This qualitative difference is caused by the geometric facts already mentioned and also by the lack of any stochastic structure in the oscillatory wave train (for the importance of this last element, see e.g., ref. 17). In either case, the effect of the bare viscosity propagates from one scale of description to the next and is never forgotten. Note that in the fluid dynamics problem incomplete similarity also plays a key role, as must be the case if the effect of the bare viscosity propagates to averaged descriptions of the phenomena.
In Fig. 6 we display the dependence of ν_{eff} on ℓ for R = 60. What is surprising at first sight is that v_{eff} is not approximately constant as ℓ changes, but, of course, it cannot be because of the boundary conditions. This is the kind of phenomenon that one encounters in wave propagation analyses of timedependent scaling, and the relationship between the present problem and those other studies, if any, remains to be elucidated. The relationship between ν_{eff} and ℓ thus displays an incomplete similarity as well.
In Fig. 7 we show the meansquare difference I at ν = v_{eff} between the mean solution ū and the solution v of the effective equation, for R = 60 and various values of ℓ. As expected, this difference goes down monotonically. If the dispersive term in our effective solution is dropped the curve is changed very little and we do not display the result. Barenblatt et al.'s conjecture that the dispersive term is not significant after averaging is thus justified.
Note that when I = I(ν_{eff}) decreases, the mean solution is increasingly well approximated by a selfsimilar solution (the solutions of vv_{x} = νv_{xx} with our boundary conditions, for different values of ν, can be mapped on each other by a change of the length scale). A selfsimilar solution is approached as ℓ increases and one approaches a fixed point of the RNG; indeed, we conjecture that selfsimilar solutions are the scaling fields of RNG fixed points (i.e., solutions that scale simply near those fixed points), and thus Barenblatt et al.'s RNG is a restriction of more general RNGs to the eigenspace of the relevant eigenvalues in the neighborhood of the fixed points. I shall discuss this issue in detail in a separate article.
Finally, we consider the validity of the asymptotic analysis of Barenblatt et al. (see Eq. 5). Barenblatt et al. looked for an equation satisfied exactly by a mean solution, whereas we are content to look for an equation that produces an approximation to that mean in an L_{2} sense; as we have just seen, the difference decreases when ℓ increases. Clearly, Eq. 5 does not hold in the major part of the wave system where v_{x} ∼ 0 when v is identified with ū for large ℓ, and thus ν_{hydro} is very large. However, one can check whether it holds in the transition region where ū ∼ v varies rapidly. The width of the transition region can be defined as x_{2} – x_{1}, where x_{1},x_{2} satisfy v(x_{1}) = 0.75 and v(x_{2}) = 0.25; ν_{hydro} is then The ratio ν_{hydro}/v_{eff} is plotted in Fig. 8 for ℓ = 20 and several values of R. The error in the ratio is substantial because x_{2} – x_{1} is not large and the number of mesh points in the transition zone is rather small (of order 10). Note that one of the assumptions in Barenblatt's theory is not satisfied; with our parameter values the width of the transition zone is comparable to the wave length of the dispersive waves. Under these conditions, the closeness of the ratio to one is a surprising confirmation of Barenblatt's theory.
These calculations were run in a computational domain of width 2X between 200 and 400 (depending on the length of the wave train), generally with 1,024 mesh points, and with a time step k equal to the smallest of: 0.3 h (to keep the advection term stable), 2h^{2}/R, 2h^{3} (to keep overall accuracy), where h is the mesh size. The calculations were run up to a time T = 200.
Conclusions
We have set up an analogy between realspace renormalization (and thus optimal prediction algorithms) and the spatial averaging of the stationary KdVB equation. The analogy leads to incomplete similarity relations between bare and renormalized or eddy diffusion, and allows us to check and extend earlier averaging calculations; it also leads to interesting conjectures regarding the relation between renormalization and selfsimilarity. A procedure similar in principle can be constructed for timedependent problems, but requires substantial further technicalities, as will be explained in subsequent work. Note that the most interesting parts of our analysis result from the fact that we are dealing with partial differential equations rather than with systems of ordinary differential equations.
The specific relation between viscosity on different scales is quantitatively different here from what is seen in hydrodynamics, as it should be in this very special model, but the model, like the Navier–Stokes equations, does exhibit a persistent effect of the bare viscosity as the scale of the model increases, expressed through an incomplete similarity. One of the reasons for picking the KdVB model for study is the fact, explained in ref. 18, that vortex stretching in the Navier–Stokes equations induces dispersive behavior and is the main source of smallscale oscillations, like the dispersion here.
Acknowledgments
I thank Professors P. Colella and M. Wright for helpful numerical advice and Professor G. I. Barenblatt for suggesting the problem and for many helpful discussions. This work was supported in part by National Science Foundation Grant DMS 9732710, and in part by the Office of Science, Office of Advanced Scientific Computing Research, Mathematical, Information, and Computational Sciences Division, Applied Mathematical Sciences Subprogram, of the U.S. Department of Energy, under Contract DEAC0376SF00098.
Footnotes

↵* Email: chorin{at}math.berkeley.edu.

Abbreviations: KdVB, Korteveg–deVries–Burgers; RNG, renormalization group.
 Copyright © 2003, The National Academy of Sciences
References
 ↵
E, W. & Engquist, B. (2003) Commun. Math. Sci. 1, 87–132.
 ↵
Rowley, C., Kevrikides, I., Marsden, J. & Lust, K. (2002) Nonlinearity 16, 1257–1282.
 ↵
 ↵
 ↵
Barenblatt, G. I. (1996) Scaling, SelfSimilarity, and Intermediate Asymptotics (Cambridge Univ. Press, Cambridge, U.K.).
 ↵
Barenblatt, G. I. (2003) Scaling (Cambridge Univ. Press, Cambridge, U.K.).
 ↵
Goldenfeld, N. (1992) Lectures on Phase Transitions and the Renormalization Group (Perseus, Reading, MA).
 ↵
Barenblatt, G. I. & Shapiro, G. I. (1984) Izv. Atmospheric Oceanic Phys. 20, 210–214.
 ↵
Barenblatt, G. I., Ivanov, M. & Shapiro, G. I. (1985) Arch. Ration. Mech. Anal. 87, 293–303.
 ↵
Bona, J. & Schonbek, M. (1985) Proc. R. Soc. Edinburgh A 101, 207–226.
 ↵
Gurevich, A. & Pitaevskii, L. (1974) Sov. Phys.JETP 38, 291–297.
 ↵
Fornberg, R. & Whitham, G. (1978) Proc. R. Soc. London A 289, 373–404.
 ↵
Zalesak, S. (1984) in Advances in Computer Methods for Partial Differential Equations V, eds. Vichnevetsky, R. & Stepleman, R. (IMACS, New Brunswick, NJ), p. 491.
 ↵
Burkhardt, T. & van Leeuwen, J. (1982) RealSpace Renormalization (Springer, Berlin).
 ↵
Kadanoff, L. (2000) Statistical Physics: Statics, Dynamics, and Renormalization (World Scientific, Singapore).
 ↵
 ↵
Lochak, P. & Meunier, C. (1988) Multiphase Averaging for Classical Systems (Springer, New York).
 ↵