BLUES iteration applied to nonlinear ordinary differential equations for wave propagation and heat transfer
BBLUES iteration applied to nonlinear ordinary differentialequations for wave propagation and heat transfer
Jonas Berx and Joseph O. Indekeu Institute for Theoretical Physics, KU Leuven, B-3001 Leuven, Belgium (Dated: December 22, 2020)
Abstract
The iteration sequence based on the BLUES (Beyond Linear Use of Equation Superposition)function method for calculating analytic approximants to solutions of nonlinear ordinary differentialequations with sources is elaborated upon. Diverse problems in physics are studied and approximateanalytic solutions are found. We first treat a damped driven nonlinear oscillator and show thatthe method can correctly reproduce oscillatory behaviour. Next, a fractional differential equationdescribing heat transfer in a semi-infinite rod with Stefan-Boltzmann cooling is handled. In thiscase, a detailed comparison is made with the Adomian decomposition method, the outcome of whichis favourable for the BLUES method. As a final problem, the Fisher equation from populationbiology is dealt with. For all cases, it is shown that the solutions converge exponentially fast tothe numerically exact solution, either globally or, for the Fisher problem, locally. a r X i v : . [ n li n . PS ] D ec . INTRODUCTION Differential equations (DEs) are the cornerstones of modern physics and mathematics.While linear DEs are quite well understood, nonlinear DEs remain elusive objects. Notwith-standing this, most physical systems exhibit some form of nonlinearity. It is a challenge toobtain exact analytical solutions whenever possible. Some analytical or semi-analytical tech-niques such as the Adomian decomposition method (ADM), the homotopy analysis method(HAM) or perturbative techniques such as the soliton perturbation theory have been pro-posed and proved to be extremely useful [1–4]. When an external force is also present inthe 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 pro-vide. While linear DEs with sources can be treated by the theory of Green functions and theprinciple of superposition, nonlinear DEs violate the superposition principle and a solutioncannot be constructed in this way. In this paper, following up on our Letter [7], we investi-gate further how the practice of Green functions can be usefully extended to nonlinear DEswith an inhomogeneous source or sink, effectively using the superposition principle beyondthe 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 reaction-diffusion-convection DEs describing traveling wave fronts. These simple solutions are exactprovided a co-moving Dirac delta function source is added to the DE. In [9] it was observedthat for some problems one and the same simple exponential tail solution simultaneouslysolves the nonlinear DE with a Dirac delta source as well as a related linear DE with thesame Dirac delta source. This observation led to the idea to formulate an analytic methodthat uses the concept of Green function, and its convolution with an arbitrary source, be-yond the linear domain. This approach was named BLUES (Beyond Linear Use of EquationSuperposition) and it was suggested that the method can be useful as a perturbation expan-sion, the small parameter being the ratio of the width of the source to the decay length ofthe tail. Subsequently, in [7] quantitative calculations were performed, demonstrating thatthe method entails an analytic iteration procedure which converges exponentially rapidlyfor a variety of problems, and which is non-perturbative. There is no need for a small pa-rameter. Furthermore, the method was also used the other way around, starting from a2inear DE and freely adding a nonlinearity to it, instead of starting from a nonlinear DE andlooking for a related linear DE. Applications were presented to solitary wave solutions ofthe Camassa-Holm equation and traveling wavefront solutions of the Burgers equation. Ourcontribution now is to extend substantially the applicability of the method and demonstrateits accuracy for problems involving oscillatory waves, phenomena described by fractionaldifferential equations (FDEs), and systems with nonlinear growth. We present a detailedcomparison with an alternative method, the ADM. In the new applications we show that theconvergence of the method is sometimes local instead of the global convergence illustratedin [7]. The applications until now were for problems that can be reduced to an ordinary DEwith a co-moving source. This limitation is henceforth removed.This paper is organized as follows. In Section II we recall briefly the analytic iterationprocedure in order to make the paper sufficiently self-contained. In Section III, we study thedamped nonlinear oscillator and show that the method allows one to approximate correctlythe decaying oscillations. In Section IV, a model for heat transfer in a semi-infinite rodwith a fractional-order derivative is treated, and employed to show that the BLUES functionmethod can be extended to the arena of fractional differential equations. Also in this Sectionwe compare the performance of the method to that of the ADM. Finally, in Section V, awave solution of a nonlinear extension of the heat equation involving reaction, i.e., the Fisherequation, 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 ofthe efficient iteration procedure that was developed from it in [7]. One starts from a linearordinary DE which can be written as an operator L z acting on a function U ( z ) and assumesthat one knows the (piecewise analytic) Green function U ( z ) = G ( z ) which solves L z U ( z ) = δ ( z ) , (1)with suitable boundary conditions. The Dirac delta source compensates a possible discon-tinuity of the derivative of order n − z = 0 in the case of an n -th order DE. Next oneconsiders the linear DE with an arbitrary source ψ ( z ). The solution is then the convolutionproduct G ∗ ψ . 3ne may add a nonlinearity rather freely, but so that the same boundary conditions arerespected, and arrive at the nonlinear ordinary DE N z U ( z ) = ψ ( z ) , (2)A solution is then proposed in the form U ( z ) = B ∗ φ . The function B ( z ) is called BLUESfunction and it is taken to be the Green function of the linear DE, so B ( z ) = G ( z ). Thechallenge 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 φ ( z ) = ψ ( z ) + R z ( B ∗ φ )( z ) . (3)To obtain the solution to the nonlinear DE (2), equation (3) can be iterated in order tocalculate an approximation in the form of a sequence in powers of the residual R z . To zerothorder, the sources φ (0) ( z ) and ψ ( z ) are identical and the approximation is the convolutionproduct U (0) ψ ( z ) = ( B ∗ φ (0) )( z ) = ( B ∗ ψ )( z ) . (4)To n th order ( n ≥ B ( z ) [7] U ( n ) ψ = ( B ∗ φ ( n ) )( z ) = U (0) ψ ( z ) + (cid:16) B ∗ R z U ( n − ψ (cid:17) ( z ) . (5)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 Diracdelta source with reduced amplitude su tt − u xx + γu x + u = s δ ( x − ct ) , (6)in which the displacement u , time t , space x , amplitude γ and velocity c are also reducedso as to be dimensionless, and look for traveling wave solutions by transforming to thecoordinate z = x − ct and restricting U ( z ) ≡ u ( x, t ), αU zz + γU z + U = s δ ( z ) , (7)4here α = c −
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, B ( z ) ≡ , z < (cid:0) λz α (cid:1) e − γz α , z ≥ , (8)with λ = (cid:112) α − γ and source amplitude s = λ/
2. Now an arbitrary nonlinear term canbe added, which we choose to be the cubic-quintic function βU + ξU , where β and ξ aretuneable parameters. Altogether the nonlinear wave equation with an arbitrary source is αU zz + γU z + U + βU + ξU = sψ ( z ) , (9)where again the amplitude is s = λ/ z ) of the source ψ ( z ) is unity. This DE is a basic model for a myriad of physical systems. When β and ξ arechosen to be − /
3! and 1 /
5! respectively, the terms U + βU + ξU can be interpreted asthe first three nonzero terms in the sine Taylor series. Including higher-order terms in theseries, one can construct the nonlinear DE αU zz + γU z + sin U = sψ ( z ) , (10)which is an equation for the damped and driven Sine-Gordon model [10], often used to de-scribe the dynamics of Josephson junctions in superconductors [11, 12]. Another importantapplication of equation (9) is the cubic-quintic Duffing oscillator, which is used to describedamped harmonic motion in a nontrivial potential and has become a paradigm for the studyof chaos. For computational convenience we will only include the terms in the sine Taylorseries up to and including third order, so we will choose ξ = 0.Following the procedure outlined in Section II, we extract the residual operator andidentify its action on U , R z U ≡ − βs U (11)and use this to calculate higher-order approximants using B ( z ) in the iteration sequence (5).Note that R z B ( z ) (cid:54) = 0. For the source ψ we can choose, e.g., an exponential corner sourcewhich possesses a tunable dimensionless decay length K , ψ ( z ) = e −| z | /K K (12)5he choice of the source (12) is not limited to functions with non-zero 1-norm; we have alsoperformed calculations for sources with 1-norm zero and found that this has no significantimpact on the convergence or accuracy of the method.The zeroth-order approximant for the solution of (9) can be calculated by performing theconvolution integral (4) with BLUES function (8) and normalized exponential corner source(12) U (0) ψ ( z ) = Kλ C + e z/K , z < (cid:2) A sin (cid:0) λz α (cid:1) − Bλ cos (cid:0) λz α (cid:1)(cid:3) e − γz α + Kλ C − e − z/K , z ≥ , (13)where the constants A, B, C ± are introduced to simplify notation. They are given by com-binations of α , γ and K : A = 2 α − K γ + 2 αK α − K γ + 2 αK + K B = K γα − K γ + 2 αK + K C ± = α ± Kγ + K . (14)Higher orders can in principle be calculated using equation (5) but in practice this is notfeasible without the aid of mathematical software.In Fig. 1, a comparison between the numerically exact solution and the zeroth-, first-and second-order approximants is made. Note that the zeroth-order approximant almostcoincides with the BLUES function (also shown), for z >
0, while the second-order approx-imant is on top of the numerically exact solution. Next, we analyse the convergence of theiteration sequence to the exact solution. In Fig. 2 the values of the approximants at z = 4for orders n = 0 , n , which providesanother means of monitoring the convergence of the method. The residual in the BLUESfunction and the residual functions for n = 0 , U ( z ),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, where6 ( z ) U num U ψ ( ) U ψ ( ) U ψ ( ) - Figure 1: Traveling wave solution to the nonlinear wave equation (9) with an exponentialcorner source (12). The numerical solution U num (red, full line), the zeroth-order U (0) ψ (black, dashed line), the first-order approximant U (1) ψ (black, wider spaced dashed line) andthe second-order approximant U (2) ψ (black, dot-dashed line, calculated for z <
9) arecompared. The BLUES function (gray, solid line) is also shown. Parameter values are α = 3, γ = 1, β = 1, ξ = 0 and K = 1 / 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 calculatethe BLUES function (8). Because the BLUES function does not solve (9), also not in thislimit, the iteration sequence converges (exponentially fast) to a nontrivial function. Onecan calculate the first two terms by setting ψ ( z ) = δ ( z ) in equations (4) and (5). To zerothorder, 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 U (1) δ ( z ) = B ( z ) − βs ∞ (cid:90) −∞ B ( z − z ) B ( z ) d z . (15)Note that this does not provide the first-order correction in a series expansion in the param-eter β because higher-order iterations generate additional contributions linear in β .7 - U ψ ( n ) ( ) - - Log | Δ U ψ ( n , nu m ) ( ) | Figure 2: Wavefront values U ( n ) ψ ( z = 4) versus order n for the nonlinear wave equation (9).The numerically exact value (red, dashed line) is also shown. Inset:
A log semi-log plotof the increments | ∆ U ( n,num ) ψ (4) | ≡ | U ( n ) ψ (4) − U ( num ) (4) | of the approximants versus n , anda linear fit. Parameter values are α = 3, γ = 1, β = 1, ξ = 0 and K = 1 / 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 <α ≤ n ≥ U ( t = 0) = C , with C ≥ N t U = D αt U + U n = ψ ( t ) , (16)where D αt is the Riemann- Liouville fractional derivative defined as follows, for α > t > D αt f ( t ) ≡ m − α ) d m dt m (cid:90) t f ( τ )( t − τ ) α +1 − m d τ, m − < α < m ∈ N d m dt m f ( t ) , α = m ∈ N (17)where Γ( . ) is the gamma function. This equation has previously been studied in the contextof nonlinear heat transfer for the case that α = 1 / C = 0 and n = 4 (Stefan-Boltzmanncooling) [13, 14]. The calculations that follow are valid for all values of 0 < α ≤ n ≥ z B R z U num R z U ψ ( ) R z U ψ ( ) R z U ψ ( ) - - - R z U ψ Figure 3: Residual function R z U ( n ) ψ of the approximants ( n = 0 , ,
2) for the nonlinearwave equation (9) with exponential corner source (12). The residual operator applied tothe BLUES function and to the numerical solution (red, full line) are also shown.Parameter values are α = 3, γ = 1, β = 1, ξ = 0 and K = 1 / ψ ( 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 thedifferential 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 differentialoperator, and a fractional part 0 < β ≤ m + 1. In the remainder of this work, we will assume0 < α ≤ U n ( t ) term. We write the linearFDE in operator form L t U = D αt U = ψ ( t ) , (18)with arbitrary source ψ ( t ) and the initial condition chosen to be U (0) = 0. The Green9unction for (18) can now be calculated by considering a Dirac delta function source insteadof ψ ( t ) D αt G ( t, t (cid:48) ) = δ ( t − t (cid:48) ) , (19)where t − t (cid:48) ≥ G ( t, t (cid:48) ) = ( t − t (cid:48) ) α − Γ( α ) (20)Consequently the solution of the linear FDE (18) is the convolution integral of the Greenfunction (20) and the source ψ ( t ), which we will choose from now on to be the constantfunction ψ ( t ) = 1 for t ≥
0, as was done in references [13, 14].It is worth emphasizing that this application to heat transfer is fundamentally differentfrom the other ones considered in this paper as well as in [7] in that the source is not assumedto be originating from a disturbance that is “co-moving” with the solution. The variablehere is time and not a co-moving coordinate, and the source arises as a natural physicalingredient of the problem. Therefore this example constitutes a non-trivial extension of thedomain of applicability of the method not only in the type of DE (from DE to FDE) butalso in the character and interpretation of the source term in the DE.We obtain U (0) ( t ) = t (cid:90) G ( t, t (cid:48) ) ψ ( t (cid:48) )d t (cid:48) = 1Γ( α ) t (cid:90) ( t − t (cid:48) ) α − d t (cid:48) = t α α Γ( α ) (21)The residual operator is the difference between the operators of the linear FDE (18) and thenonlinear FDE (16) and is defined by the action on U ( t ), i.e., R t U = L t U − N t U = − U n (22)Now the p − th order approximant to equation (16) can be calculated by using the BLUESiteration sequence (5) U ( p ) ( t ) = U (0) ( t ) + t (cid:90) G ( t, t (cid:48) ) R t (cid:48) U ( p − ( t (cid:48) )d t (cid:48) (23)The first-order approximant to the nonlinear problem can easily be calculated using (20),1021) and the iteration sequence definition (23), with the choice n = 4, U (1) ( t ) = U (0) ( t ) − t (cid:90) ( t − t (cid:48) ) α − Γ( α ) t (cid:48) α α Γ ( α ) d t (cid:48) = 1 α Γ( α ) t α − Γ(1 + 4 α ) α Γ ( α )Γ(1 + 5 α ) t α (24)One can now iterate (23) to generate higher-order approximants to the solution of thenonlinear FDE (16). In Fig. 4 and Fig. 5, the approximants for different values of α arecompared with the numerically exact solution.For the choice α = 1 /
2, equation (16) is associated with the heat transfer equations fora 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 ≥ n ≤
2, all energy is eventuallyradiated away. For the remainder of this work, we will use nonlinear Stefan-Boltzmanncooling, n = 4. The zeroth-, first-, and second-order approximants for α = 1 / n = 4are, respectively, given by U (0) ( t ) = 2 (cid:18) tπ (cid:19) / U (1) ( t ) = 2 (cid:18) tπ (cid:19) / − (cid:18) tπ (cid:19) / U (2) ( t ) = 2 (cid:18) tπ (cid:19) / − (cid:18) tπ (cid:19) / + 20971524725 (cid:18) tπ (cid:19) / − (cid:18) tπ (cid:19) / + 8796093022208369208125 (cid:18) tπ (cid:19) / − (cid:18) tπ (cid:19) / (25)The approximants obtained by the BLUES function method can be compared to thoseobtained with the Adomian decomposition method (ADM) by making use of the followingrecursion relation [13] for the coefficients a n of the solution series of the ADM, U ADM ( t ) = ∞ (cid:88) n =0 a n t αn , (26)with a = 0 a = 1Γ( α + 1) a n +1 = − Γ( nα + 1)Γ( nα + α + 1) A n (27)11 ( ) U ( ) U ( ) U ( ) U ( ) ( t ) (a) U ( ) U ( ) U ( ) U ( ) U ( ) ( t ) (b) Figure 4: Solution to the nonlinear FDE (16) with constant source ψ ( t ) = 1 and fractionalorder (a) α = 1 / (b) α = 1 /
2. The numerical solution (red, full line) is comparedwith the approximants up to fourth order (black/gray lines).12 ( ) U ( ) U ( ) U ( ) U ( ) ( t ) (a) U ( ) U ( ) U ( ) U ( ) U ( ) t U ( t ) (b) Figure 5: Solution to the nonlinear FDE (16) with constant source ψ ( t ) = 1 and fractionalorder (a) α = 3 / (b) α = 1. The numerical solution (red, full line) is compared withthe approximants up to fourth order (black/gray lines).13here the A n , n ≥ A n = n (cid:88) l =0 l (cid:88) s =0 s (cid:88) 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 / a n = 0 for n (cid:54) = 4 k + 1, with k = 0 , , , ... . We will henceforth define the p -th order ADM approximantas the truncated series U ( p )ADM ( t ) ≡ p (cid:88) n =0 a n t αn , (29)A comparison is now made between BLUES and ADM. In Fig.6 the 21st-order approxi-mant for the ADM (containing five nontrivial exact terms) is shown, together with the 4thapproximant for the BLUES function method (containing many more terms, but also fivenontrivial exact ones – see further) and the numerically exact solution. ADM 21st orderBLUES 4th iteration t U ( t ) Figure 6: Comparison between the 21st-order approximant of the ADM (dotted line) andthe 4th iteration of the BLUES function method (dashed line) for α = 1 / n = 4. Thenumerically exact solution (red, full line) is also shown.The number of nonzero terms g ( i, n ) for the i -th approximant generated by the BLUES14unction method can be calculated for a general nonlinearity exponent n ∈ N , n ≥
1, i.e., g ( i, n ) = n − ( n − n i ) n ≥ i + 1 n = 1 (30)so for n = 4, the 4th approximant already contains g (4 ,
4) = 86 nonzero terms. Notethat the number of terms in BLUES increases exponentially with the order of iteration. Incomparison, the ADM generates a series with a number of terms which grows linearly withthe order of the approximation. Note that the ADM generates the exact coefficients in aseries expansion of the solution while the BLUES function method does not. In contrast,BLUES generates many more terms, of which only the lower-order ones are exact. Thehigher-order coefficients have not yet settled or converged to their exact value. We findempirically that in the BLUES function method in iteration p the first p + 1 coefficients areexact, which is a linear progression like in ADM. For example, in (25) the expressions for U (0) ( t ) and U (1) ( t ) contain the exact coefficients, while only the first three terms in U (2) ( t )are exact. Although the higher-order terms are not exact, it appears that their presence(in BLUES) leads to a more accurate approximant than their absence (in ADM), when thecomparison is made with equal numbers of exact terms.The BLUES function method generates, in each iteration, a (huge) number of scout termsthat probe the emerging series expansion and gradually gain precision. In this respect, theBLUES function method is reminiscent in spirit of a Pad´e approximation applied to a seriesexpansion. The coefficients a n in the BLUES function method saturate roughly linearlywith increasing order of iteration, as can be seen, e.g., when keeping track of the coefficient a of the t / term. This coefficient is first generated in the second approximant U (2) with a provisional value of a = 1 . U (3) the value of thiscoefficient increases to a = 19 . U (4) , attainingits exact value a = 30 . a of the t / term. This term is first generatedin the second approximant with a value of a = − . a = − . a = − . a = − .
387 in the fifth iteration, which is the exact result asfound by the ADM.The higher-order residual functions are shown in Fig.7 together with the residual operator1522) applied to the numerically exact solution (red, full line). Note that for this model theresiduals are not localized, in contrast to all previous examples. The approximants to thesolution of equation (16) do not converge to the correct numerical value at t → ∞ owing tothe divergence of the residual in every order of iteration. While the approximants divergefor larger values of t , a radius of convergence can still be identified. R t U ( ) R t U ( ) R t U ( ) R t U ( ) - - - - - R t U ( t ) t Figure 7: Residual R t U ( t ) = − U ( t ) for α = 1 / 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 describesthe propagation of a density u ( x, t ) u t − νu xx = 0 , (31)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, L z U ≡ − U z − kU zz = δ ( z ) , (32)16here 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 pointincluding z = 0) is the piecewise analytic exponential tail, B ( z ) = , z < − z/k , z ≥ c ( k ) = 1 /k . We now add a reaction-type nonlinearity (growthterm) and a source ψ ( z ) to obtain the forced Fisher equation [17] in co-moving coordinates,i.e., N z U = − U z − kU zz − kU (1 − U ) = ψ ( z ) (34)with boundary conditions U ( z → ∞ ) → U ( z → −∞ ) →
1. Equation (34) governs thedimensionless density of some bio-chemical substance or biological population experiencingdiffusion and growth. The limit at negative infinity signifies the saturation of the density atthe normalized value 1. The residual operator is now acting as follows R z U = kU (1 − U ) (35)If we now choose ψ ( z ) to be the exponential corner source (12), the zeroth-order approximantto the nonlinear DE (34) is the same as was calculated for the Burgers equation in [7], and isrepeated here in appendix A. Higher orders can easily be calculated by iteration but will notbe given here. We will, however, present the calculated first-order approximant for k (cid:54) = K in the appendix. Note that the first-order approximant approaches a constant which is notunity at negative infinity for z and consequently does not obey the boundary condition. Onecan calculate the non-trivial constant by considering the limit of the first-order approximantat negative infinite z , i.e., U c ≡ lim z →−∞ U (1) ( z ) = 1 + lim z →−∞ ∞ (cid:90) −∞ B ( z − z (cid:48) ) R z (cid:48) U (0) ψ ( z (cid:48) )d z (cid:48) = 1 + ∞ (cid:90) −∞ R z (cid:48) U (0) ψ ( z (cid:48) )d z (cid:48) = 1 + k (2 k + 4 k K + 6 kK + 3 K )4( k + K ) (36)In Fig.8a, the zeroth- and first-order approximants are shown together with the numeri-cally exact solution and the BLUES function (33). In Fig.8b, a zoomed-in representation of17he shoulder of the wavefront is shown. The numerical solution is compared with approxi-mants 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 → −∞ . Thisis a consequence of the lack of localization of the residual for higher orders. This is shownin Fig.9. The numerical residual function R z U num and the zeroth-order residual functionare localized, but higher-order residual functions are not anymore. This corresponds to adivergence of the approximants of higher orders.The local convergence of the approximants can nevertheless be assessed by studying thevalue of the approximants for a fixed value of z . In view of the divergence of the approximantsfor z → −∞ , one has to be careful in choosing the value of z . In this case, we have optedfor z = −
1, which can be seen to lie within a reasonable region of convergence. The resultsare shown in Fig.10. Note that within the region of convergence, the approximants convergeexponentially fast to the numerically exact solution.
VI. CONCLUSIONS AND OUTLOOK
In this paper we have extended a useful and accurate approach for solving nonlinearordinary DEs with sources, based on the BLUES function method introduced in [9] andfirst applied quantitatively in [7]. We have shown that the method can be applied to obtainuseful approximations for the oscillating wave solution of the damped nonlinear oscillatorand for traveling wavefront solutions of the Fisher equation. It can also be extended tothe arena of fractional DEs, as was shown for a heat transfer problem. In problems ofwave propagation, typically a partial DE is transformed to an ordinary DE by means of atraveling wave Ansatz and the source must be assumed to be co-moving. This limitationhas, however, been removed by considering a broader class of physical problems, as we haveillustrated 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 detailedcomparison of our iteration procedure with the Adomian decomposition method. It is foundthat both methods generate equal numbers of exact coefficients in a given order of the seriesexpansion (ADM) or of the iteration (BLUES). However, the BLUES method generates anexponentially large number of approximate coefficients of higher-order terms that are notpresent in the ADM. The presence of these higher-order terms appears to accelerate the18 ( z ) U num U ψ ( ) U ψ ( ) - - - - (a) U num U ψ ( ) U ψ ( ) U ψ ( ) U ψ ( ) U ψ ( ) - - - - (b) Figure 8: (a)
Traveling wavefront solution to the nonlinear Fisher DE (34) with anexponential corner source (12). The numerical solution U num (red, full line), thezeroth-order U (0) ψ (black, dashed line) and the first-order approximant U (1) ψ (black, widerspaced dashed line) are compared. The BLUES function (gray) is also shown. (b) Azoomed-in view around the shoulder of the wavefront. The approximants are shown up toand including 4th order. The latter (gray, full line) lies just below the numerical one (red,full line). Parameter values are k = 1 / K/k = 1 / z U num R z U ψ ( ) R z U ψ ( ) R z U ψ ( ) R z U ψ ( ) - - - - R z U ψ Figure 9: Residual function R z U ( n ) ψ of the approximants of order n = 0 , , , k = 1 / K/k = 1 / t and space x , which cannot be reduced to ordinary20 U ψ ( n ) ( - ) - - Log | Δ U ψ ( n , n - ) ( - ) | Figure 10: Wavefront values U ( n ) ψ ( −
1) at z = − n for the Fisher equation(34). The numerically exact value (red, dashed line) is also shown. Inset:
A log semi-logplot of the increments | ∆ U ( n,n − ψ ( − | of the approximants versus n , and a linear fit.These increments are defined as the difference of the n -th and ( n − k = 1 / K/k = 1 / Appendix A: First-order approximant for the Fisher equation
The zeroth-order approximant to the solution of the nonlinear Fisher equation (34) withexponential corner source (12) is U (0) ψ ( z ) ≡ ( B ∗ ψ )( z ) = 12 − KK + k e z/K , z < KK − k e − z/K − k K − k e − z/k , z ≥ , (A1)21hile for K = k the convolution product results in U (0) ψ ( z ) ≡ ( B ∗ ψ )( z ) = − e z/k , z < (cid:0) + z k (cid:1) e − z/k , z ≥ . (A2)To calculate the first-order approximation to the solution of the Fisher equation (34), theresidual operator (35) is applied to the zeroth-order approximant in accordance with theiteration equation (5) R z U (0) ψ ( z ) = − k Ke z/K ( Ke z/K − k + K ) ) ( k + K ) , z < , e − z ( k + 1 K )( K ( k + K ) e z/k − k e z/K )( ( k + K ) e z/k ( K − K − k ) e z/K ) − k e z/K ) ( k − K ) , z ≥ U (1 , ψ = U (1) ψ − U (0) ψ is presented here, which is the convolution product ∆ U (1 , ψ ( z ) = ( B ∗ R z U (0) ψ )( z )∆ U (1 , ψ ( z ) = k α (cid:16) K k + K e z/K − K e z/K + 2(2 k + 4 k K + 6 kK + 3 K ) (cid:17) , z < , k α β (cid:16) k e − z/k + k K (5 k + K ) γ e − z/k − k K e − αz/kK +4 K α e − z/K + K α k − K e − z/K (cid:17) , z ≥ , (A4)where α = k + K , β = k − K and γ = K − k . The cases for which K = k or K = 2 k mustbe treated separately but these (elementary) calculations will not be performed here. Notethat the first-order approximant U (1) ψ ( z ) approaches the constant value U c for z → −∞ ,which has already been presented in equation (36). [1] G. Adomian, Solving Frontier Problems of Physics: The Decomposition Method , 1st ed., Fun-damental Theories of Physics, Vol. 60 (Springer Netherlands, 1994).[2] S. Liao,
Homotopy Analysis Method in Nonlinear Differential Equations , 1st ed. (Springer-Verlag Berlin Heidelberg, 2012).[3] V. I. Karpman and E. M. Maslow, “Perturbation theory for solitons,” J. Exp. Theor. Phys. , 537–559 (1977).
4] J. P. Keener and D. W. McLaughlin, “Solitons under perturbations,” Phys. Rev. A , 777–790 (1977).[5] K. Maeda and T. Colonius, “A source term approach for generation of one-way acoustic wavesin the Euler and Navier-Stokes equations,” Wave Motion , 36 (2017).[6] A. Ermakov and Y. Stepanyants, “Soliton interaction with external forcing within theKorteweg-de Vries equation,” Chaos , 013117 (2019).[7] J. Berx and J. O. Indekeu, “Analytic iteration procedure for solitons and traveling wavefrontswith sources,” J. Phys. A-Math. Theor. , 38LT01 (2019).[8] J. O. Indekeu and R. Smets, “Traveling wavefront solutions to nonlinear reaction-diffusion-convection equations,” J. Phys. A-Math. Theor. , 315601 (2017).[9] J. O. Indekeu and K. K. M¨uller-Nedebock, “BLUES function method in computationalphysics,” J. Phys. A-Math. Theor. , 165201 (2018).[10] E. G. Ekomasov, A. M. Gumerov, R. V. Kudryavtsev, S. V. Dmitriev, and V. N. Nazarov,“Multisoliton dynamics in the Sine-Gordon model with two point impurities,” Braz. J. Phys. , 576 (2018).[11] Z. Gul, A. Ali, and A. Ullah, “Localized modes in parametrically driven long Josephsonjunctions with a double-well potential,” J. Phys. A-Math. Theor. , 015203 (2018).[12] I. O. Starodub and Y. Zolotaryuk, “Fluxon interaction with the finite-size dipole impurity,”Phys. Lett. A , 1419 – 1426 (2019).[13] J.-S. Duan, R. Rach, D. Baleanu, and A. Wazwaz, “A review of the Adomian decompositionmethod and its applications to fractional differential equations,” Commun. Frac. Calc. ,73–99 (2012).[14] A. Wazwaz and S. Khuri, “A reliable technique for solving the weakly singular second-kindVolterra-type integral equations,” Appl. Math. Comput. , 287 – 299 (1996).[15] K. Padmavally, “On a non-linear integral equation,” J. Math. Mech. , 533–555 (1958).[16] J. Keller and W. Olmstead, “Temperature of a nonlinearly radiating semi-infinite solid,” Q.Appl. Math. , 559–566 (1972).[17] R. Fisher, “The wave of advance of advantageous genes,” Ann. Eugen. , 355–369 (1937).[18] A.-M. Wazwaz, Partial Differential Equations and Solitary Waves Theory , Nonlinear PhysicalScience (Springer Berlin Heidelberg, Berlin, Heidelberg, 2009).
19] Y. Yang, Z. Yan, and D. Mihalache, “Controlling temporal solitary waves in the generalizedinhomogeneous coupled nonlinear Schr¨odinger equations with varying source terms,” J. Math.Phys. , 053508 (2015).[20] Z. Yan, X.-F. Zhang, and W. M. Liu, “Nonautonomous matter waves in a waveguide,” Phys.Rev. A , 023627 (2011).[21] T. J. Kippenberg, A. L. Gaeta, M. Lipson, and M. L. Gorodetsky, “Dissipative Kerr solitonsin optical microresonators,” Science (2018).[22] X. Ma, O. A. Egorov, and S. Schumacher, “Creation and manipulation of stable dark solitonsand vortices in microcavity polariton condensates,” Phys. Rev. Lett. , 157401 (2017).[23] M. V. Berry, “Minimal analytical model for undular tidal bore profile; quantum and Hawkingeffect analogies,” New J. Phys. , 053066 (2018).[24] M. V. Berry, “Minimal model for tidal bore revisited,” New J. Phys , 073021 (2019).[25] M. Kardar, G. Parisi, and Y.-C. Zhang, “Dynamic scaling of growing interfaces,” Phys. Rev.Lett. , 889–892 (1986).[26] K. A. Takeuchi, “An appetizer to modern developments on the Kardar-Parisi-Zhang univer-sality class,” Physica A , 77 (2018).[27] J. M. Horowitz and M. Kardar, “Bacterial range expansions on a growing front: Roughness,fixation, and directed percolation,” Phys. Rev. E , 042134 (2019)., 042134 (2019).