The BLUES function method applied to partial differential equations and analytic approximants for interface growth under shear

An iteration sequence based on the BLUES (beyond linear use of equation superposition) function method is presented for calculating analytic approximants to solutions of nonlinear partial differential equations. This extends previous work using this method for nonlinear ordinary differential equations with an external source term. Now, the initial condition plays the role of the source. The method is tested on three examples: a reaction-diffusion-convection equation, the porous medium equation with growth or decay, and the nonlinear Black-Scholes equation. A comparison is made with three other methods: the Adomian decomposition method (ADM), the variational iteration method (VIM), and the variational iteration method with Green function (GVIM). As a physical application, a deterministic differential equation is proposed for interface growth under shear, combining Burgers and Kardar- Parisi-Zhang nonlinearities. Thermal noise is neglected. This model is studied with Gaussian and space-periodic initial conditions. A detailed Fourier analysis is performed and the analytic coefficients are compared with those of ADM, VIM, GVIM, and standard perturbation theory. The BLUES method turns out to be a worthwhile alternative to the other methods. The advantages that it offers ensue from the freedom of choosing judiciously the linear part, with associated Green function, and the residual containing the nonlinear part of the differential operator at hand.


I. INTRODUCTION
It is a challenge, in exact sciences and theoretical physics in particular, to obtain useful analytical approximations to solutions of nonlinear differential equations (DEs). In this context the Adomian decomposition method (ADM), the homotopy analysis method (HAM) or perturbative techniques such as the soliton perturbation theory have proven useful [1][2][3][4].
In two recent papers [5,6], we demonstrated how the practice of Green functions can be usefully extended to nonlinear ordinary differential equations (ODEs) that are inhomogeneous, featuring a source or sink, by effectively using the superposition principle beyond the linear domain. In the present paper, we extend the approach to nonlinear partial differential equations (PDEs) and present an application to the physics of interface growth in a soft condensed matter system under shear flow.
To situate this development, we briefly recall the history of the BLUES function method.
In [7] exponential tail solutions of nonlinear reaction-diffusion-convection ODEs describing traveling wave fronts with co-moving sources were studied. In [8] it was noted that an exponential tail solution may simultaneously solve the nonlinear ODE and a related linear ODE, both with a co-moving Dirac delta source. This led to an analytic method that uses the Green function beyond the linear domain, named BLUES ("Beyond Linear Use of Equation Superposition"). Next, in [5,6] it was shown how to develop the method into a non-perturbative and rapidly converging analytic iteration procedure. One may start from a linear DE and freely add a nonlinearity. Applications were given to solitary waves, oscillatory waves, nonlinear growth and transport of heat, and the method was extended to fractional ODEs and to sources that need not be co-moving. Now, we extend the approach to nonlinear PDEs, e.g., in time t and one space coordinate x, which cannot be reduced to ODEs. For PDEs the initial condition serves as the source and no external source must be added. We will compare the BLUES iteration with four other methods: the Adomian decomposition method (ADM) [1,9], the variational iteration method (VIM) [10], the VIM with Green function (GVIM) [11], and straightforward perturbation theory (PT).
The setup of this work is as follows. In Section II we extend the BLUES function method to the arena of PDEs in two variables, one of which is time. We restrict our attention in this paper to operators with a first derivative in time. In Section III we illustrate the method for three simple exactly solvable PDEs and compare the different methods. In Section IV we set the stage for a physical problem by applying the method to a general power-law convective nonlinearity. Next, in Section V we introduce and study a simple model for the time evolution of a growing fluid interface under shear. In Section VI we conclude and present an outlook.

II. THE BLUES FUNCTION METHOD FOR A NONLINEAR PDE
Here we extend the BLUES iteration method originally developed for ODEs [5,6,8]

to
PDEs in time and one space variable. The crucial role of the extrinsic source (or sink) term in the context of the ODE will now be taken over, simply, by the intrinsic initial condition of the solution of the PDE. Consequently, the extension of the method to PDEs entails a conceptual simplification rather than complication, and allows one to increase substantially the range of physics problems that can be tackled.
Let us start from a linear PDE which can be written as an operator L t,x acting on a function u(x, t), say a density subject to diffusion, and let us attempt to solve L t,x u(x, t) = 0, for t > 0, (1) with initial condition Since the problem is linear the solution u(x, t) can be written as the convolution G * f of the initial condition f (x) with the Green function G(x, t), which satisfies L t,x G(x, t) = 0, for t > 0, with Dirac-delta initial condition The solution to the linear problem is the (single-variable) convolution, which reads For simplicity we restrict our attention to PDEs that involve only the first derivative w.r.t. to time, specifically L t,x u = u t +L x u, with u t ≡ ∂u/∂t andL x a time-independent linear operator. For our purposes, it is convenient to rewrite the PDE by invoking the initial condition f (x) through the action of a Dirac-delta source in time. The following time and space integral, which is a two-variable convolution u(x, t) = G * f δ, solves the rearranged inhomogeneous linear PDE, which is equivalent to the original linear PDE, This identity holds by virtue of the fact that L t,x contains only a first derivative w.r.t.
time t. This derivative generates two terms. The boundary term (the value of the integrand at t = t) exactly produces the right-hand-side of (6), in view of (4). The second term is contained in the action of L t,x , when it is moved inside the integral over t . That contribution, however, vanishes as one can verify by careful inspection. We conclude that G * f δ solves the PDE for all t > 0.
The initial condition is retrieved by examining the limit t → 0. Firstly, the solution u(x, t) as given by the time and space integral G * f δ obviously vanishes for t < 0 − by definition, so u(x, t < 0) = 0. However, this solution "jumps" to the initial condition function f (x) at t = 0 + through the action of δ(t ) and by the fact that the Green function becomes a spatial Dirac-delta in view of (4). The space integral then produces f (x). For t > 0 the solution evolves, in a continuous manner, from this initial condition.
Using this representation of the PDE, which naturally features an intrinsic source term expressing the initial condition, we can now generalize the BLUES iteration procedure from nonlinear ODEs to nonlinear PDEs. One may add a nonlinearity rather freely to the PDE, while preserving the simple form of the time-dependent part, withÑ x a time-independent nonlinear operator, and arrive at the nonlinear PDE with intitial condition, as before, The BLUES function method now proposes to construct a solution u(x, t) to the equiv- Clearly, this PDE coincides with the original nonlinear PDE (8) for t > 0 and we will shortly examine its behavior at t = 0. The function B(x, t) is called BLUES function and it is taken to be the Green function of an arbitrary but conveniently chosen linear operator L t,x related to N t,x . The challenge is to calculate the new associated source φ(x, t) knowing that B * f δ solves the linear PDE (6) with initial condition f (x) and source term f δ. Note that φ(x, t) need not be separable and in general it is not.
The initial condition is generated correctly, since, by definition, u(x, t < 0) = 0 and subsequently u(x, t = 0 + ) = f (x), provided three conditions are fulfilled. The first is that N t,x u = 0, for u = 0. The second condition is that the associated source φ decomposes as follows into a separable singular term and a (non-separable) smooth term ζ, which is to be calculated analytically: The third condition is that for all finite x the functionÑ x f (x) be finite. For nonlinear operators these are not obvious and must be checked.
For this calculation one defines a (time-independent) residual operator R x ≡ L t,x − N t,x and makes use of the implicit identity which follows directly from the Green function property of B w.r.t. the chosen linear PDE.
To obtain the solution to the nonlinear PDE (8) with initial condition (9), equation (11) can be rewritten and iterated, in order to calculate an approximation in the form of a sequence in powers of the residual R x . In zeroth iteration, and in nth iteration (n ≥ 1), Consequently, the nth analytical approximant to the solution of the nonlinear PDE is found through the two-variable convolution

III. TEST CASES FOR THE METHOD
A. Reaction-diffusion-convection equation Let us start with a simple example, in which the convolutions are all of single variable type. Unless otherwise stated the functions, variables and parameters are reduced (dimensionless). Consider the nonlinear reaction-diffusion-convection PDE [12] which can be used to describe, e.g., the propagation of a chemical of density u through the combined mechanisms of diffusion, nonlinear convection and reaction, This unbounded initial condition is rather unphysical but will serve as an ideal testbed for the comparison of the different approximation methods, as in this case a simple exact solution of (16) can be found. We will now consider the methods mentioned in Section I and compare their results. The ADM and VIM both produce the following sequence of approximants, . . .
which converges slowly to the exact solution Note that the sequence (18) . .
which converges to the exact solution (19) for n → ∞ as well.
We now turn to the BLUES function method, and follow the scheme outlined in Section II. First, the PDE (16) with initial condition f (x) is rewritten as follows defined on (x, t) ∈ R × [0, ∞) and the initial condition u(x, 0) = f (x) = e −x has been converted to a source term by multiplication with a Dirac-delta function in the temporal coordinate. Choosing the linear operator simple and without spatial derivatives, one can define the associated linear PDE with source ψ(x, t) ≡ f (x)δ(t) as follows, which is solved by u(x, t) = f (x)G(t), with G(t) the Green function for L t . Note that we omitted the linear term u xx from the linear part L t,x of the operator N t,x . This judicious choice, which is a distinct feature of the BLUES strategy, not only simplifies the calculations but also considerably improves the convergence.
We obtain a step function with exponential tail, and the solution U (t) for the linear problem with arbitrary source ψ(t), for t > 0, is since G(τ < 0) = 0 and s > 0.
We next define the residual operator R x as the difference between the linear and the nonlinear operator, i.e., R x = L t − N t,x , so and set up the iteration sequence based on (14) and (15) for the solution to (21), where the BLUES function B(τ ) is the Green function G(τ ) of (23) for the chosen linear operator L t , whose action is given in (22). The zeroth approximant is the convolution of the BLUES function with the source ψ(x, t), Iterating through the procedure (26), one finds the following sequence of approximants . . .
which converges to the exact solution (19) for n → ∞. Note that each approximant is bounded and useful for all t by virtue of the overall factor e −2t .
We can now compare the results of the three different methods. Since all three methods converge to the known exact solution (19), one can define an error function E (n) (x, t) as the absolute value of the difference between the nth approximant and the exact solution u ex , which for n = 4 and x = 1 results in (16e) −1 . Note that the errors for both the ADM and VIM and for the GVIM are monotonically increasing in time and hence the approximations decrease in accuracy for large values of t. In contrast, for the BLUES function method the error vanishes in the limit t → ∞ and this method provides the fastest convergence for all t > 0. The reason for this improved performance is that the choice of the linear operator part in the BLUES function method is free and can be tailored so as to render all the approximants well bounded for all times.

B. Porous medium equation with growth or decay
The second example is in the realm of fluid mechanics: the nonlinear porous medium equation [13] with linear growth or decay, with m > 1 and β ∈ R. We consider a density w(x, t) in one space dimension with initial condition w(x, 0) = f (x) = x. Unless otherwise stated the functions, variables and parameters are reduced (dimensionless). We will only consider a quadratic nonlinearity, m = 2, which allows us to write (31) as follows The components of the solution generated by the ADM are and in the limit n → ∞ this converges to the exact solution where the sign of β indicates whether there is growth or decay. Note that the ADM generates term by term the exact coefficients of the powers of t in the Taylor expansion in time of the solution.
The VIM produces the following sequence of approximants to the solution of (32), which also converges to the exact solution (35). Note that VIM and ADM produce different results. The VIM does not immediately give the exact coefficients but recursively adjusts them until they saturate at the exact value.
Next, the GVIM produces the sequence . . .
For n ≥ 2, the approximants (37) are invariable. The GVIM in this case produces the exact solution (35) already in the second iteration and contributions from higher iterations are zero.
We now turn to the BLUES function method. The PDE (31) with initial condition w(x, 0) = f (x) can be rewritten as a nonlinear PDE with a source ψ( Choosing the linear operator to be of the same form as the successful one used in the previous section, one can define the associated linear PDE with the same source term, and we recall the Green function for this linear operator, Note that in this case the linear operator is chosen by simply dropping (only) the nonlinear term in N t,x . We now obtain the residual operator R x , which acts as follows on the function w, and set up the iteration sequence for the solution to (38) where the BLUES function B(τ ) is the Green function G(τ ) of (40) for the chosen linear operator L t , whose action is given in (39). The zeroth approximant is the convolution of the BLUES function and the source ψ(x, t), i.e., Iterating further according to the procedure (42), one finds the following sequence of approximants for m = 2 . . .
which, remarkably, produces the exact solution (35) to (31) already in the first iteration.
Higher iterations remain at this "fixed point". In Fig.2 we compare the results from each of the above methods and also compare their errors, at the level of this first iteration.

C. Nonlinear Black-Scholes equation
For the following example, let us look at the field of economics. Unless otherwise stated the functions, variables and parameters are reduced (dimensionless). The Black-Scholes equation describes the value V (S, τ ) of an option for some underlying asset price S ∈ [0, ∞) over a period τ ∈ [0, T ], with T the time of maturity, that is, the last moment on which an option can be exercised. After expiration or maturity, the option contract will cease to exist and the buyer cannot exercise their right to buy or sell. The underlying asset price S is a stochastic variable and follows a geometric Brownian motion. In [14], the authors consider a nonlinear Black-Scholes PDE for V (S, τ ), which assumes that the market is incomplete through the combined feedback effects of illiquid markets and large trader effects. In this PDE S is treated as a continuous variable, which we name s, and s and τ are treated as independent variables. This PDE is the following, with t the time until expiry, t = T − τ , u the value function, u(s, t) ≡ V (S, τ ), σ the volatility, r the risk-free interest rate. The constant ρ is a measure of the liquidity of the market. In order to ensure that feedback effects from hedging generate so-called volatility smiles, one has to choose this liquidity parameter to be negative [15,16]. We consider the initial condition u(s, is the starting price of the asset. In [17], the authors study the solution of (45) by means of the ADM. This gives the following sequence of component functions of the solution, The solution is the sum of all the component functions u i (s, t), This claim can easily be verified by noticing that the component functions u i (s, t) are the coefficients of the Taylor series of the exact solution [14], The VIM produces the following sequence of approximants to the solution of (45), which converges slowly to the exact solution (48).
Next, the GVIM produces the following iterates Finally, we study the BLUES method. As usual, we first rewrite equation (45) with the inclusion of a source ψ(s, t) = f (s)δ(t), i.e., and consider the associated linear operator we have used in the previous examples together with the source ψ(s, t), with Green function, Note that in this example, the linear operator is chosen judiciously by not only dropping the nonlinear term but some linear terms as well. Hence, the residual, whose action is defined still contains two linear terms. The zeroth approximant is the convolution of the BLUES function (53) and the source ψ(s, t), The BLUES function method generates the following sequence of approximants In Fig. 3 we compare the results from each of the above methods and also compare their errors, at the level of the 3rd approximant or 3rd order (n = 3). Note that we have not If one were to fix T > 0 at a finite value, it is obvious that the accuracy of the approximate solutions for all of the above procedures decreases for t → T , i.e., for increasing remaining time until end of contract.

IV. DIFFUSION EQUATION WITH GENERAL NONLINEARITY
We now set the stage for the analysis of a nonlinear PDE associated with a simple physical model for the growth of an interface between two fluids that are subject to shear flow, by first considering a more general nonlinear PDE from a technical viewpoint. The heat equation with diffusion constant D > 0 and general nonlinearity u m u n x , where m, n ≥ 0 is given by the PDE, with Gaussian initial condition u( and boundary conditions u(|x| → ∞, t) = 0. As before, we adopt the notation N t,x u to denote the nonlinear operator acting on u(x, t). The associated linear PDE of our choice is the one-dimensional heat equation describing normal diffusion, with the same initial condition and the same boundary conditions. This linear PDE has Green function In the small time limit t → 0, the Green function (60) approaches a Dirac-delta distribution δ(x). The solution to the diffusion equation with the Gaussian initial condition f (x) can be calculated by convoluting f (x)δ(t) with the kernel G(x, t), Integrating over time and space gives which is itself a decaying Gaussian with mean zero and with variance Σ 2 (t) ≡ σ 2 + 2Dt.
This solution u (0) serves as the zeroth iteration in the BLUES scheme. One now considers the residual operator R x = L t,x − N t,x which can be applied to the zeroth approximant (62), Convoluting the previous expression with the Green function (60) results in where S 2 (t, s) ≡ 2D(t − s) + Σ 2 (s)/(m + n), which can be interpreted as a variance. Further, c(t, s) ≡ (Σ 2 (s)/S 2 (t, s))/(m + n) and α(t, s) ≡ (m + n)(S 2 (t, s)/Σ 2 (s))/(4D(t − s)). The spatial integral can be calculated exactly where Γ(n) is the gamma function and 1 F 1 (a, b, z) is the confluent hypergeometric function of the first kind [18]. This spatial integral can equivalently be expressed in terms of the Hermite polynomials H n (z) in the following way, Ξ(x, t, s, m, n) ≡ −i 2 n π α(t, s) n+1 H n i α(t, s)c(t, s)x .
We list here the following useful properties for the hypergeometric functions and for the Hermite polynomials: The first correction to the zeroth approximant (62) now becomes For some choices of (m, n) this can be simplified greatly. In the next section we discuss a physical system which features two such cases combined, (1, 1) and (0, 2).

V. INTERFACE GROWTH UNDER SHEAR
We propose a minimalistic model for the growth of an interface between two fluids near two-phase coexistence and subject to an externally imposed shear flow. On the one hand, we exploit the finding that the growing interface between a stable and an unstable domain in a kinetic Ising model at low temperature can be described by including in the effective growth equation a Kardar-Parisi-Zhang (KPZ) nonlinearity which allows for lateral growth [19][20][21][22].
On the other hand, we make use of the growth equation proposed for studying interface fluctuations under shear flow, including a Burgers type of nonlinearity [23] which allows for a background linear shear flow imposed on the phase-separated fluid [24,25]. We combine the two growth equations but limit ourselves to the minimal setting of two-dimensional systems (i.e., a one-dimensional interface) and the deterministic version of the equation. We ignore thermal noise and postpone an application to the stochastic DE until later work.
Our starting point is, as usual, the Edwards-Wilkinson equation for interface growth [26], which, in its deterministic version, reads where h(x, t) is the height of an interface that fluctuates, measured relative to a (horizontal) straight reference line (along x). This reference line is co-moving with the growing interface and therefore a velocity term v is omitted in (72). D is a diffusion coefficient (proportional to the interfacial tension whose action is to smoothen the interface).
A cartoon of the physical setting is shown in Fig. 4. Following Bray et al. [24,25] we include an externally imposed shear flow. The motivation, in part, for this was that there is an interesting subtle competition between the smoothing of an interface under shear and the roughnening of an interface under thermal noise. Later studies elucidated interface confinement under shear using Monte Carlo simulation [27,28]. Incorporating a (horizontal) shear velocity profile v x (y) amounts to invoking the total time derivative, Ah + B and we can choose a reference frame co-moving at the mean velocity, so B = 0. We thus add a Burgers convective nonlinearity to the PDE.
Next, following Devillard and Spohn [20] we recognize that the interface growth, ignoring the lattice anisotropies of the model, is in the direction normal to the local tangent. This growth, in which a stable domain overtakes an unstable one, is driven by a pressure difference, or chemical potential difference, with respect to two-phase coexistence (i.e., a non-zero external magnetic field in the Ising model). Incorporating this lateral growth amounts to invoking the KPZ geometric correction, where v is the velocity of the growing interface. Since the term v is already absorbed in (72) we need to add only the gradient-squared term to the PDE. Altogether we obtain the nonlinear PDE where A is the shear rate.
This PDE combines the Burgers and KPZ nonlinearities but, we recall, ignores thermal noise. When taken separately, each of these two nonlinearities amount to exactly solvable PDEs, but to our knowledge not when combined. This makes it worthwhile to derive a useful analytical approximant to the solution of the combined equation. Note that in our physical context extra terms proportional to h or h 2 are not present in (75) because in the absence of shear flow we require translational invariance of the growth equation along the y-direction. In addition, we require translational invariance along x. Also note that in terms of the scaling properties of interface growth the Burgers term is the dominant perturbation [24,25] and the KPZ term is subsidiary. We do not discuss these properties here.
There is an alternative route to the PDE (75) which is worth pointing out. One may start from the stochastic KPZ equation for interface growth and couple it to the stochastic Navier-Stokes (NS) equation for the velocity field v, by replacing the time derivative in KPZ by the total time derivative, as in (73), and invoking the NS equation for v. This system of coupled DEs was proposed and studied in [29]. If, in that system, one ignores the random force in the stochastic NS equation and imposes a (deterministic) shear flow velocity profile, and if one also ignores thermal noise in the KPZ equation, one arrives again at (75).
We now proceed to the calculations and adapt the notation slightly in order to be conform with that of previous sections. We define the nonlinear operator, acting on the function with α and β real parameters. For the linear operator L t,x we choose the entire linear part of N t,x , which is the linear diffusion operator. The residual operator R x (cf. Section IV), is then defined through By doing so, the nonlinear problem would be suited to be tackled by perturbation theory (PT), if the terms that feature the parameters α and β can be considered to be small compared to the terms of the linear part. This brings us in position to compare the BLUES iteration, which is non-perturbative, to a direct perturbation expansion, keeping in mind that the former makes no assumptions on the magnitude of the nonlinear terms. What we find is akin to our observations in the treatment of ODEs [6]. The BLUES iteration generates a sequence that is in general different from summing up the terms a series expansion, except possibly in the first iteration in which the BLUES result may coincide with that of 1st-order PT.
We consider two different initial conditions, corresponding to distinct physical situations.
The first is a single (Gaussian) interface protrusion or "bump", for which we will illustrate the method at the level of the zeroth and first iteration only, and show its close similarity to 1st-order PT. The second one is a (sinusoidal) periodic interface front, for which we will study the time evolution to higher level in the iteration scheme. For that case, we will perform a detailed comparison of the results from ADM, VIM, GVIM, BLUES and PT.

A. Gaussian initial condition
First, we will consider the situation of a solitary interface bump that can be modeled by a Gaussian initial condition u(x, 0) = f (x), given in equation (58). We assume the boundary conditions u(|x| → ∞) = 0. The associated linear PDE is the heat equation (59).
The zeroth approximant is now the decaying Gaussian solution (62) of the linear equation.
Using equation (71) twice, once for the convective nonlinearity (Burgers) and once for the nonlinear lateral growth (KPZ), the first approximant can be calculated analytically. We report here the result (a detailed calculation can be found in Appendix A), Note that the effects introduced by the convective nonlinearity contain only odd functions of x, and the effects introduced by the nonlinear growth contain only even functions of x.
In the first iteration the effect of nonlinearity is a simple superposition of the individual nonlinear effects, i.e., nonlinear convection and nonlinear growth. Only in higher iterations does the interplay (mixing) between these different effects take place.
At this level of approximation, the BLUES approximant u (1) coincides with the result of straightforward PT to first order in α and β. This is not surprising in view of the fact that the chosen residual operator coincides with the nonlinear part of the differential operator, which is precisely the "perturbation" when α and β are considered small. We have also performed the ADM and VIM calculations for this case. These methods are, however, not suitable here because they produce large oscillations that grow uncontrollably both in time and in higher orders of approximation. We will return to these methods when we consider a periodic interface undulation.
In the first iteration of the nonlinear problem we obtain, This is a non-decreasing function of time for β < 0, hence the bump grows as a consequence of the lateral growth correction, even when there is no overall (vertical) growth along y in the co-moving frame. Note that the parameter α does not enter the equation. The shear flow only moves particles along x and does not influence the bump size but only its shape.
In Fig. 5  These results coincide with, respectively, those of standard zeroth and 1st-order PT in the parameters α and β. In the course of time the bump distorts to the right (for α > 0) and grows somewhat (for β < 0) until its size saturates.

B. Space-periodic initial condition
For convenience and simplicity, in this example we will work with dimensionless variables x and t, as well as dimensionless u, D, α and β. To study a space-periodic interface contour, we can choose the following trigonometric initial condition f (x) One can now apply the residual operator (77) to (81). After simplifying the result by using trigonometric power reduction identities, the residual is The first approximant to the solution of equation (76) can be calculated by convoluting the residual (82) with the Gaussian Green function, making use of the following identities Hence, the first approximant is Higher approximants can be calculated with moderate effort. In Fig. 7 we show the first three BLUES approximants together with the numerically exact solution for a fixed time t = 1/3.
Next, in Fig.8 we compare the numerical solution and the fourth BLUES approximant with the 4th-order VIM and ADM results at t = 1/3.  Parameters are D = α = 1 and β = −1.
so we will consider the following perturbation expansion u P for the solution u P T (x, t) = u 0,0 (x, t) + αu 1,0 (x, t) + α 2 u 2,0 (x, t) + βu 0,1 (x, t) and we assume, within PT, to avoid ambiguity, that α and β are of the same order of magnitude. Performing the expansion and solving the resulting linear PDEs yields the expressions given in Appendix B for the perturbative solution u P T up to, and including, second order in α and β.
Note that PT generates terms of second order in α and β, i.e., α 2 , β 2 and αβ, and Fourier modes up to and including the third harmonic (with respect to the period of the initial condition). In contrast, in the second iteration the BLUES function method does not yet provide the exact coefficients of the 2nd-order terms. Furthermore, this method also generates terms of higher order in α and β, e.g., α 3 , β 3 , α 2 β, etc., and Fourier modes of the fourth harmonic are also already present in the second approximant. We provide also the full expressions of the 2nd approximant in Appendix B and compare them quantitatively with PT.
In Fig.9, we compare the second BLUES approximant with 2nd-order PT at t = 2/3. Finally, in Fig. 10 we show the various n = 2 approximations (ADM, VIM and BLUES) for a fixed spatial coordinate x = π. We remark that the 2nd-order approximations for the ADM and VIM coincide exactly for x = π.  are given by with c p = (a p − ib p )/2. In Fig.11  It is conspicuous that BLUES iteration progresses differently from PT. There is even a qualitative difference. For long times the asymptotic behavior of the BLUES approximants agrees with the numerically exact solution in that all the harmonics decay to zero. This is not always the case in the PT (e.g., for |c 1 (t)| and |c 2 (t)| in Fig. 12).
An interesting quantity is the asymptotic "size excess" ∆ of the solution as a consequence of the lateral growth correction of the interface. This is given by the long-time limit of c 0 (t), The numerically obtained precise value for the size excess is ∆ num = 0.2356, while the nth BLUES approximants give ∆

VI. CONCLUSIONS
The extension of the BLUES function method to PDEs (in time and one other variable) presented here represents a significant broadening of the scope of the method. In previously reported applications to ODEs, a (co-moving) source term had to be added to the differential equation, corresponding to a physical input external to the problem and inevitably somewhat ad hoc. In contrast, in the present application to PDEs, the source term is a natural intrinsic ingredient, being the initial condition of the problem.
In its formulation for a PDE, the BLUES iteration can be compared with various other approaches, and we have made this comparison for 4 methods: the ADM, VIM, GVIM and PT. We have observed that the BLUES iteration often provides better convergence towards the (numerically exact) solution, and offers a qualitative advantage in attaining correct asymptotic behavior for long times. This favourable position appears to result from the freedom in the method to tailor the linear operator part of the problem, so that an "optimal" Green function becomes available (which forces correct asymptotics), and to add the remainder of the DE as a residual that contains the nonlinear operator part but also whatever remaining linear part that was not chosen to be incorporated in the linear operator.
This freedom is instrumental, and requires physical insight on the side of the user before unleashing the calculation.
A physical application has been supplied, which deals with the motion and growth of an interface in a phase-separated fluid subject to shear flow. In a minimalistic model, we have combined, at the deterministic level (low-temperature approximation) the features of KPZlike lateral growth and Burgers-like convection due to shear. A detailed study was made of a space-periodic version, with comprehensive Fourier analysis of the evolving contour. In this example, a scrutiny was made of the similarities, and differences, between (non-perturbative) BLUES and PT. This comparison has turned out, once again, to favour the former.
For the future, we envision an extension of the method to stochastic DEs for which the noise can play the role of an external source. We also consider an application to coupled DEs for which the Green function is a matrix exponential. In closing we note that our restriction, in this paper, to a first-order time derivative, is not necessary. We announce that the method can be readily adapted to study second-order time derivatives, by converting the problem to a system of coupled first-order PDEs. The initial conditions for the function and its derivative can then both be included, after suitable multiplication with a delta function.
Finally, combining equations (A3) and (A6), the first correction to the zeroth approximant becomes which can easily be solved and subsequently simplified by noticing that 2Σ 2 (t) − Σ 2 (2t) = σ 2 to give the following expression for the correction in first iteration to the zeroth approximant of (76) This can now be rearranged to yield equation (78).

Appendix B: Fourier coefficients for the space-periodic interface contour
In this Appendix we discuss in detail the Fourier coefficients of various harmonics that are generated by the BLUES iteration at the level of the second approximant (n = 2) for the problem of the time evolution of the periodic interface contour and compare them with their counterparts in 2nd-order PT. We first present, for p ∈ {0, 1, 2, 3}, the real pth coefficients calculated by both methods and then discuss them with the aid of two figures, 13 and 14.
• b 1 BLUES: b 1 (t) = e −Dt − (α 2 + 4β 2 )(e −Dt − 2e −3Dt + e −5Dt ) 32D 2 (B10) Note that in the second approximant for a 0 terms of order α 2 β and β 3 are generated, which are absent in 2nd-order PT. Als note that a 1 (first harmonic) and a 3 (third harmonic) are both proportional to αβ, as in PT. Importantly, in the BLUES function method a 1 (t) and a 3 (t) tend to zero for long times, in agreement with the numerical solution, whereas the 2nd-order PT expressions tend to non-zero constants (see also Fig. 13). In this respect the BLUES iteration is qualitatively superior. The coefficient a 2 (t) (second harmonic) has a first order in β contribution which is the same in both methods, and an additional α 2 β contribution in the second BLUES approximant. In both methods the result is very close to the numerical solution (see Fig. 13). Note that a 4 (t) (fourth harmonic) is generated in 2nd-iteration BLUES but is absent in 2nd-order PT. This is a consequence of the fact that the BLUES function method is non-perturbative and already generates higher harmonics in a lower iteration than the perturbation series.
As for the b n (t), the coefficient b 1 (t) (first harmonic reflecting the initial condition) contains the zeroth approximant, which is (of course) the same in both methods. Moreover, the entire expressions for b 1 (t) coincide in 2nd-iteration BLUES and 2nd-order PT (see also Fig. 14). The coefficient b 2 (t) (second harmonic) has a first order in α contribution which is the same in both methods, and an additional αβ 2 contribution in the 2nd BLUES approximant.
In both methods the result is nearly the same but both are somewhat off of the numerical solution (see Fig. 14). Importantly, in the BLUES function method b 3 (t) (third harmonic) tends to zero for long times, in agreement with the numerical solution, whereas the 2ndorder PT expression tends to a non-zero constants (see also Fig. 14). In this respect the BLUES iteration is again qualitatively superior. Finally, b 4 (t) (fourth harmonic) is present