BLUES iteration applied to nonlinear ordinary differential equations for wave propagation and heat transfer

The iteration sequence based on the BLUES (Beyond Linear Use of Equation Superposition) function method for calculating analytic approximants to solutions of nonlinear ordinary differential equations with sources is elaborated upon. Diverse problems in physics are studied and approximate analytic solutions are found. We first treat a damped driven nonlinear oscillator and show that the method can correctly reproduce oscillatory behaviour. Next, a fractional differential equation describing heat transfer in a semi-infinite rod with Stefan-Boltzmann cooling is handled. In this case, a detailed comparison is made with the Adomian decomposition method, the outcome of which is favourable for the BLUES method. As a final problem, the Fisher equation from population biology is dealt with. For all cases, it is shown that the solutions converge exponentially fast to the numerically exact solution, either globally or, for the Fisher problem, locally.


I. INTRODUCTION
Differential equations (DEs) are the cornerstones of modern physics and mathematics.
While linear DEs are quite well understood, nonlinear DEs remain elusive objects. Notwithstanding this, most physical systems exhibit some form of nonlinearity. It is a challenge to obtain exact analytical solutions whenever possible. Some analytical or semi-analytical techniques such as the Adomian decomposition method (ADM), the homotopy analysis method (HAM) or perturbative techniques such as the soliton perturbation theory have been proposed and proved to be extremely useful [1][2][3][4]. When an external force is also present in the form of an inhomogeneous source or sink [5,6], one often turns to numerical solutions.
These solutions do not offer the richness of information that analytical solutions can provide. While linear DEs with sources can be treated by the theory of Green functions and the principle of superposition, nonlinear DEs violate the superposition principle and a solution cannot be constructed in this way. In this paper, following up on our Letter [7], we investigate further how the practice of Green functions can be usefully extended to nonlinear DEs with an inhomogeneous source or sink, effectively using the superposition principle beyond the linear domain.
Relative to previous works on this method our contribution is situated as follows. In [8] the usefulness was demonstrated of simple exponential tail solutions of nonlinear reactiondiffusion-convection DEs describing traveling wave fronts. These simple solutions are exact provided a co-moving Dirac delta function source is added to the DE. In [9] it was observed that for some problems one and the same simple exponential tail solution simultaneously solves the nonlinear DE with a Dirac delta source as well as a related linear DE with the same Dirac delta source. This observation led to the idea to formulate an analytic method that uses the concept of Green function, and its convolution with an arbitrary source, beyond the linear domain. This approach was named BLUES (Beyond Linear Use of Equation Superposition) and it was suggested that the method can be useful as a perturbation expansion, the small parameter being the ratio of the width of the source to the decay length of the tail. Subsequently, in [7] quantitative calculations were performed, demonstrating that the method entails an analytic iteration procedure which converges exponentially rapidly for a variety of problems, and which is non-perturbative. There is no need for a small parameter. Furthermore, the method was also used the other way around, starting from a linear DE and freely adding a nonlinearity to it, instead of starting from a nonlinear DE and looking for a related linear DE. Applications were presented to solitary wave solutions of the Camassa-Holm equation and traveling wavefront solutions of the Burgers equation. Our contribution now is to extend substantially the applicability of the method and demonstrate its accuracy for problems involving oscillatory waves, phenomena described by fractional differential equations (FDEs), and systems with nonlinear growth. We present a detailed comparison with an alternative method, the ADM. In the new applications we show that the convergence of the method is sometimes local instead of the global convergence illustrated in [7]. The applications until now were for problems that can be reduced to an ordinary DE with a co-moving source. This limitation is henceforth removed. This paper is organized as follows. In Section II we recall briefly the analytic iteration procedure in order to make the paper sufficiently self-contained. In Section III, we study the damped nonlinear oscillator and show that the method allows one to approximate correctly the decaying oscillations. In Section IV, a model for heat transfer in a semi-infinite rod with a fractional-order derivative is treated, and employed to show that the BLUES function method can be extended to the arena of fractional differential equations. Also in this Section we compare the performance of the method to that of the ADM. Finally, in Section V, a wave solution of a nonlinear extension of the heat equation involving reaction, i.e., the Fisher equation, is treated. We close the paper with our conclusions and outlook.

II. THE BLUES ITERATION METHOD
Here we recapitulate concisely the BLUES function concept [9] from the viewpoint of the efficient iteration procedure that was developed from it in [7]. One starts from a linear ordinary DE which can be written as an operator L z acting on a function U (z) and assumes that one knows the (piecewise analytic) Green function U (z) = G(z) which solves with suitable boundary conditions. The Dirac delta source compensates a possible discontinuity of the derivative of order n − 1 at z = 0 in the case of an n-th order DE. Next one considers the linear DE with an arbitrary source ψ(z). The solution is then the convolution product G * ψ.
One may add a nonlinearity rather freely, but so that the same boundary conditions are respected, and arrive at the nonlinear ordinary DE A solution is then proposed in the form U (z) = B * φ. The function B(z) is called BLUES function and it is taken to be the Green function of the linear DE, so B(z) = G(z). The challenge is to calculate the associated source φ(z) for the given source ψ(z), knowing that B * ψ solves the linear DE with source ψ(z). For achieving this one defines a residual operator R z ≡ L z − N z and makes use of the implicit identity To obtain the solution to the nonlinear DE (2), equation (3) can be iterated in order to calculate an approximation in the form of a sequence in powers of the residual R z . To zeroth order, the sources φ (0) (z) and ψ(z) are identical and the approximation is the convolution To nth order (n ≥ 1), the approximate solution can be found by iterating (3) and taking the convolution product with B(z) [7] U (n) We now proceed to applications beyond what was presented in [7].

III. THE DAMPED NONLINEAR OSCILLATOR
We start from a general linear wave equation in one dimension with a co-moving Dirac delta source with reduced amplitude s in which the displacement u, time t, space x, amplitude γ and velocity c are also reduced so as to be dimensionless, and look for traveling wave solutions by transforming to the coordinate z = x − ct and restricting U (z) ≡ u(x, t), where α = c 2 − 1. Derivatives (with respect to t, x or z) have been denoted by subscripts.
This DE is solved by the Green function, to be used as BLUES function, with λ = 4α − γ 2 and source amplitude s = λ/2. Now an arbitrary nonlinear term can be added, which we choose to be the cubic-quintic function βU 3 + ξU 5 , where β and ξ are tuneable parameters. Altogether the nonlinear wave equation with an arbitrary source is where again the amplitude is s = λ/2 and the 1-norm (i.e., the integral over z) of the source ψ(z) is unity. This DE is a basic model for a myriad of physical systems. When β and ξ are chosen to be −1/3! and 1/5! respectively, the terms U + βU 3 + ξU 5 can be interpreted as the first three nonzero terms in the sine Taylor series. Including higher-order terms in the series, one can construct the nonlinear DE which is an equation for the damped and driven Sine-Gordon model [10], often used to describe the dynamics of Josephson junctions in superconductors [11,12]. Another important application of equation (9) is the cubic-quintic Duffing oscillator, which is used to describe damped harmonic motion in a nontrivial potential and has become a paradigm for the study of chaos. For computational convenience we will only include the terms in the sine Taylor series up to and including third order, so we will choose ξ = 0.
Following the procedure outlined in Section II, we extract the residual operator and identify its action on U , and use this to calculate higher-order approximants using B(z) in the iteration sequence (5).
Note that R z B(z) = 0. For the source ψ we can choose, e.g., an exponential corner source which possesses a tunable dimensionless decay length K, The choice of the source (12) is not limited to functions with non-zero 1-norm; we have also performed calculations for sources with 1-norm zero and found that this has no significant impact on the convergence or accuracy of the method.
The zeroth-order approximant for the solution of (9) can be calculated by performing the convolution integral (4) with BLUES function (8) and normalized exponential corner source where the constants A, B, C ± are introduced to simplify notation. They are given by combinations of α, γ and K: Higher orders can in principle be calculated using equation (5) but in practice this is not feasible without the aid of mathematical software.
In Fig. 1, a comparison between the numerically exact solution and the zeroth-, firstand second-order approximants is made. Note that the zeroth-order approximant almost coincides with the BLUES function (also shown), for z > 0, while the second-order approximant is on top of the numerically exact solution. Next, we analyse the convergence of the iteration sequence to the exact solution. In Fig. 2  where U > 0, that they are zero (with a cubic dependence on z) where the curves of Fig.   1 change sign, and that they are positive in the domain around their first minimum, where Figure 1: Traveling wave solution to the nonlinear wave equation (9) with an exponential corner source (12). The numerical solution U num (red, full line), the zeroth-order U  U < 0. This is obvious in view of the simple form (11).
In the limit K → 0, the chosen source converges to the Dirac delta source used to calculate the BLUES function (8). Because the BLUES function does not solve (9), also not in this limit, the iteration sequence converges (exponentially fast) to a nontrivial function. One can calculate the first two terms by setting ψ(z) = δ(z) in equations (4) and (5). To zeroth order, the convolution is identical to the BLUES function. To first order in the iteration, the solution with a Dirac delta function source is given by Note that this does not provide the first-order correction in a series expansion in the parameter β because higher-order iterations generate additional contributions linear in β.

IV. FRACTIONAL HEAT TRANSFER EQUATION
In this section we consider the following nonlinear ordinary fractional differential equation (FDE) for U (t) defined on the semi-infinite real line t ∈ [0, ∞) with differential order 0 < α ≤ 1 and exponent n ≥ 1 and with initial condition U (t = 0) = C 0 , with C 0 ≥ 0.
where D α t is the Riemann-Liouville fractional derivative defined as follows, for α > 0 and t > 0, where Γ(.) is the gamma function. This equation has previously been studied in the context of nonlinear heat transfer for the case that α = 1/2, C 0 = 0 and n = 4 (Stefan-Boltzmann cooling) [13,14]. The calculations that follow are valid for all values of 0 < α ≤ 1 and n ≥ 1.
ψ of the approximants (n = 0, 1, 2) for the nonlinear wave equation (9) with exponential corner source (12). The residual operator applied to the BLUES function and to the numerical solution (red, full line) are also shown.
It has been shown that if ψ(t) is a piecewise continuous bounded function, equation (16) is guaranteed to have a unique solution. If ψ(t) is nondecreasing in an interval 0 < t < s, s ∈ (0, ∞) then the solution is also nondecreasing in that interval [15,16]. Note that the differential order can in principle be higher than α = 1. One can then separate the order α = m + β in an integer part m ∈ N corresponding to a regular integer-order differential operator, and a fractional part 0 < β ≤ 1 which again corresponds to a fractional differential operator. The FDE (16) should consequently be supplemented with additional boundary (or initial) conditions up to a number m + 1. In the remainder of this work, we will assume 0 < α ≤ 1.
Following the steps in the BLUES procedure outlined in Section II, we can simply start from the linear DE obtained by dropping the nonlinear U n (t) term. We write the linear FDE in operator form with arbitrary source ψ(t) and the initial condition chosen to be U (0) = 0. The Green function for (18) can now be calculated by considering a Dirac delta function source instead of ψ(t) where t − t ≥ 0 because the problem is formulated on the semi-infinite real line. This Green function is readily calculated to be Consequently the solution of the linear FDE (18) is the convolution integral of the Green function (20) and the source ψ(t), which we will choose from now on to be the constant function ψ(t) = 1 for t ≥ 0, as was done in references [13,14].
It is worth emphasizing that this application to heat transfer is fundamentally different from the other ones considered in this paper as well as in [7] in that the source is not assumed to be originating from a disturbance that is "co-moving" with the solution. The variable here is time and not a co-moving coordinate, and the source arises as a natural physical ingredient of the problem. Therefore this example constitutes a non-trivial extension of the domain of applicability of the method not only in the type of DE (from DE to FDE) but also in the character and interpretation of the source term in the DE.
We obtain The residual operator is the difference between the operators of the linear FDE (18) and the nonlinear FDE (16) and is defined by the action on U (t), i.e., Now the p−th order approximant to equation (16) can be calculated by using the BLUES iteration sequence (5) The first-order approximant to the nonlinear problem can easily be calculated using (20), (21) and the iteration sequence definition (23), with the choice n = 4, One can now iterate (23) to generate higher-order approximants to the solution of the nonlinear FDE (16). In Fig. 4 and Fig. 5, the approximants for different values of α are compared with the numerically exact solution.
For the choice α = 1/2, equation (16) is associated with the heat transfer equations for a semi-infinite solid [16] with external heating ψ(t) and either linear Newton (for n = 1) or nonlinear Stefan-Boltzmann cooling (for n = 4). In [16], it was shown that for n ≥ 3, some of the energy entering the solid will remain, while for n ≤ 2, all energy is eventually radiated away. For the remainder of this work, we will use nonlinear Stefan-Boltzmann cooling, n = 4. The zeroth-, first-, and second-order approximants for α = 1/2 and n = 4 are, respectively, given by The approximants obtained by the BLUES function method can be compared to those obtained with the Adomian decomposition method (ADM) by making use of the following recursion relation [13] for the coefficients a n of the solution series of the ADM, with a 0 = 0  where the A n , n ≥ 1 are defined as follows A n = n l=0 l s=0 s k=0 a n−l a l−s a s−k a k (28) From equations (27) and (28) it can be deduced that for α = 1/2 one has a n = 0 for n = 4k + 1, with k = 0, 1, 2, .... We will henceforth define the p-th order ADM approximant as the truncated series A comparison is now made between BLUES and ADM. In Fig.6 the 21st-order approximant for the ADM (containing five nontrivial exact terms) is shown, together with the 4th approximant for the BLUES function method (containing many more terms, but also five nontrivial exact ones -see further) and the numerically exact solution.  The higher-order residual functions are shown in Fig.7 together with the residual operator for larger values of t, a radius of convergence can still be identified. Figure 7: Residual R t U (t) = −U 4 (t) for α = 1/2 for different orders of iteration. The numerically exact residual R t U num is also shown (red, full line).

V. FISHER EQUATION
As a starting point for our final example, consider the diffusion equation which describes the propagation of a density u(x, t) with ν the (dimensionless) diffusion coefficient. Adopting a traveling-wave Ansatz z = x−ct, the diffusion equation becomes a linear ordinary DE. We add a co-moving Dirac delta source, where k is a dimensionless constant. We consider the wavefront boundary conditions U z (z → −∞) = 0 (and U (z → −∞) > 0) and U (z → ∞) = 0. The exact solution (in every point including z = 0) is the piecewise analytic exponential tail, and the wavefront velocity is c(k) = 1/k. We now add a reaction-type nonlinearity (growth term) and a source ψ(z) to obtain the forced Fisher equation [17] in co-moving coordinates, i.e., with boundary conditions U (z → ∞) → 0 and U (z → −∞) → 1. Equation (34) governs the dimensionless density of some bio-chemical substance or biological population experiencing diffusion and growth. The limit at negative infinity signifies the saturation of the density at the normalized value 1. The residual operator is now acting as follows If we now choose ψ(z) to be the exponential corner source (12), the zeroth-order approximant to the nonlinear DE (34) is the same as was calculated for the Burgers equation in [7], and is repeated here in appendix A. Higher orders can easily be calculated by iteration but will not be given here. We will, however, present the calculated first-order approximant for k = K in the appendix. Note that the first-order approximant approaches a constant which is not unity at negative infinity for z and consequently does not obey the boundary condition. One can calculate the non-trivial constant by considering the limit of the first-order approximant at negative infinite z, i.e., In Fig.8a, the zeroth-and first-order approximants are shown together with the numerically exact solution and the BLUES function (33). In Fig.8b, a zoomed-in representation of the shoulder of the wavefront is shown. The numerical solution is compared with approximants up to fourth order. Note that while all approximants obey the boundary condition U (z → ∞) → 0, only the zeroth-order approximant approaches unity for z → −∞. This is a consequence of the lack of localization of the residual for higher orders. This is shown in Fig.9. The numerical residual function R z U num and the zeroth-order residual function are localized, but higher-order residual functions are not anymore. This corresponds to a divergence of the approximants of higher orders.
The local convergence of the approximants can nevertheless be assessed by studying the value of the approximants for a fixed value of z. In view of the divergence of the approximants for z → −∞, one has to be careful in choosing the value of z. In this case, we have opted for z = −1, which can be seen to lie within a reasonable region of convergence. The results are shown in Fig.10. Note that within the region of convergence, the approximants converge exponentially fast to the numerically exact solution.

VI. CONCLUSIONS AND OUTLOOK
In this paper we have extended a useful and accurate approach for solving nonlinear ordinary DEs with sources, based on the BLUES function method introduced in [9] and first applied quantitatively in [7]. We have shown that the method can be applied to obtain useful approximations for the oscillating wave solution of the damped nonlinear oscillator and for traveling wavefront solutions of the Fisher equation. It can also be extended to the arena of fractional DEs, as was shown for a heat transfer problem. In problems of wave propagation, typically a partial DE is transformed to an ordinary DE by means of a traveling wave Ansatz and the source must be assumed to be co-moving. This limitation has, however, been removed by considering a broader class of physical problems, as we have illustrated in the case of the nonlinear ordinary heat transfer DE.
In the application, in Section IV, to the heat transfer problem we have made a detailed comparison of our iteration procedure with the Adomian decomposition method. It is found that both methods generate equal numbers of exact coefficients in a given order of the series expansion (ADM) or of the iteration (BLUES). However, the BLUES method generates an exponentially large number of approximate coefficients of higher-order terms that are not present in the ADM. The presence of these higher-order terms appears to accelerate the   ψ of the approximants of order n = 0, 1, 2, 3 for the nonlinear Fisher DE with exponential corner source (12). The residual operator applied to the numerical solution (red, full line) is also shown. Parameter values are k = 1/3 and convergence of the approximation.
One can think of further applications of the method. For an overview of pertinent types of DEs for which this method could be tried, see for example [18]. We have in mind, for example, systems described by the nonlinear Schrödinger (or Gross-Pitaevskii) equation.
This DE has been studied intensively in the past decades and is of general interest in the research on Bose-Einstein condensates, nonlinear optics and semiconductor physics [19][20][21][22].
In the field of fluid mechanics, when studying solitary waves on shallow water, it would be interesting to see whether the BLUES function method can be used to calculate physically relevant profiles of tidal bores in the context of the recent minimal analytic model for this problem [23,24]. Other opportunities arise when one considers different kinds of DEs, e.g., nonlinear partial DEs, either deterministic or stochastic [25,26] and systems of coupled DEs [27].
Finally, we announce that the method is also applicable to nonlinear partial differential equations (PDEs), for example in time t and space x, which cannot be reduced to ordinary DEs. For PDEs the initial condition can be employed as the source term in the method. This is the subject of a substantial further development on which we will report in a future communication. The increment between the zeroth-order and first-order approximants ∆U is presented here, which is the convolution product ∆U 2k+K e 2z/K − 4K 3 e z/K + 2(2k 3 + 4k 2 K + 6kK 2 + 3K 3 ) , z < 0, k 8α 2 β 2 4k 5 e −2z/k + 8k 2 K 3 (5k 2 +K 2 ) γ e −z/k − 8k 2 K 3 e −αz/kK +4K 3 α 2 e −z/K + K 4 α 2 2k−K e −2z/K , z ≥ 0, where α = k + K, β = k − K and γ = K 2 − 4k 2 . The cases for which K = k or K = 2k must be treated separately but these (elementary) calculations will not be performed here. Note that the first-order approximant U (1) ψ (z) approaches the constant value U c for z → −∞, which has already been presented in equation (36).