BLUES function method in computational physics
Joseph O. Indekeu Kristian K. Müller-Nedebock




Institute for  Theoretical Physics, KU Leuven, B-3001 Leuven, Belgium


Department of  Physics, Stellenbosch University, Stellenbosch, South Africa




We introduce a computational method in physics that goes “beyond linear use of equation superposition" (BLUES). A BLUES function is defined as a solution of a nonlinear differential equation (DE) with a delta source that is at the same time a Green's function for a related linear DE. For an arbitrary source, the BLUES function can be used to construct an exact solution to the nonlinear DE with a different, but related source. Alternatively, the BLUES function can be used to construct an approximate piecewise analytical solution to the nonlinear DE with an arbitrary source. For this alternative use the related linear DE need not be known. The method is illustrated in a few examples using analytical calculations and numerical computations. Areas for further applications are suggested.


BLUES function method in computational physics
    Kristian K. Müller-Nedebock
    December 30, 2023
==============================================





§ SETTING THE STAGE: AN EXAMPLE FROM POPULATION DYNAMICS

Our aim in this contribution is to provide a method that can assist in constructing a solution to a nonlinear differential equation (DE) with a source term using superposition. In the recent past some efforts in this direction were made, notably by Frasca who developed an extension of the Green's function method in the domain of nonlinear differential equations or nonlinear systems <cit.>. Here, we propose a different approach. In order to set a concrete stage for illustrating the method we give an example in the field of population dynamics <cit.> and the physics of active matter and recall the Fisher equation <cit.>, 

    ∂ u/∂ t = α u(1-u) + D ∂^2u/∂ x^2,

which describes the growth of a (dimensionless) concentration u(x,t)  up to a certain level (here normalized to 1) and its diffusion. The growth rate is α and the diffusion constant is D. After the transformation: t ←α t, x ←√(α/D) x, we arrive at a dimensionless representation,

    ∂ u/∂ t = u(1-u) + ∂^2u/∂ x^2


The following generalization of the Fisher equation allows for convective motion <cit.>, 

    ∂ u/∂ t + ∂ h(u)/∂ x = u(1-u) + ∂^2u/∂ x^2


For the function h(u) we consider the form h(u) = a_0 + a_1 u + a_2 u^2. The first-order term represents a constant convective velocity, which only shifts the propagation velocity. We therefore retain the quadratic term and write h(u) =  k u^2/2 <cit.>, where k is a positive dimensionless constant. 

It was recently observed that when a delta-source term in the co-moving frame is added,

    ∂ u/∂ t + k u ∂ u/∂ x - u(1-u) - ∂^2 u/∂ x^2 =  1/kδ(x-c(k) t),

where the function c(k) is chosen to match the propagation velocity self-consistently, 
the DE possesses the following traveling wave solution with an exponential tail <cit.>,

    u(x,t) ≡ U(z) =
    
    1,     z<0 
    
    e^-z/k,     z ≥ 0,

where z = x-ct, with c the velocity of the traveling wave, and k is required to be positive to provide a physical solution. The wave front head is at z=0. It is readily verified that (<ref>) solves (<ref>) with the dispersion relation

    c(k) = k + 1/k

The amplitude of the delta-function source is determined by integrating the DE, expressed in the co-moving frame, from z= 0^- to z=0^+ and requiring that the result be an identity. For further convenience, we rewrite the partial DE as an ordinary DE in the variable z,

    k(k U -c(k))∂ U/∂ z - kU(1-U) - k∂^2 U/∂ z^2 =  δ(z),

and we specify once and for all the boundary conditions that we are interested in, U(z→ -∞) = 1 and U(z→∞) = 0.

Our paper is organized as follows. Section 2 deals with the derivation of a suitable linear DE with a delta source, that is solved by the same piecewise analytic function as the nonlinear DE, which is therefore a Green's function. Section 3 defines the BLUES function and explains how it can be used using superposition of solutions. Section 4 identifies the source which must be combined with the nonlinear DE in order that the piecewise analytic superposition be an exact solution of the nonlinear DE. Section 5 presents an example of a  DE with a gradient-squared nonlinearity for which the same method can be applied. Section 6 gives a perturbation theory viewpoint in which our approach corresponds to the leading (zeroth-order) contribution. This zeroth-order approximation can be used in quite general circumstances. Section 7 illustrates the potential of the method in a few test cases by comparing the analytic zeroth-order solution with the numerical solution. A brief conclusion and outlook on further applications is given in Section 8. 



§ DERIVING A RELATED LINEAR DIFFERENTIAL EQUATION

We now attempt to linearize the DE in U about U=1. We write U= 1 +ϵ and obtain the following linear DE in zeroth order (i.e., neglecting ϵ but not its derivatives),

    k(k-c(k))∂ U/∂ z -  k∂^2 U/∂ z^2=  δ(z)

Remarkably, this DE is again solved by (<ref>) with the same dispersion relation (<ref>). 

Making use of (<ref>) we can now define a Green's function G with the help of the linear operator

    L_z ≡ -∂/∂ z - k ∂^2 /∂ z^2,

so that 

    L_z G(z-z_0) =  δ(z-z_0)

Clearly, the Green's function G(z) coincides with (<ref>).
The solution U^(ℓ)_ϕ of the linear DE L_z U^(ℓ)_ϕ (z) = ϕ(z) with an arbitrary source ϕ(z) is then given by

    U^(ℓ)_ϕ(z) =  ∫_-∞^∞ G(z-z_0) ϕ (z_0) dz_0 ≡ (G∗ϕ)(z)

For example, for a source that is an exponentially decaying “corner" with decay constant K and amplitude A,

    ϕ (z) = A e^-|z|/K,

the solution to the linear DE is

    U^(ℓ)_ϕ(z) =
    
    AK (2-K/K+k e^z/K ),     z<0 
    AK/K-k (K e^-z/K -  2k^2/K+k e^-z/k ),     z ≥ 0,

This solution, as well as its first and second derivatives, are continous in z=0. The third derivative jumps from -AK^-1(K+k)^-1 at z=0^- to A(2K+k)(Kk)^-1(K+k)^-1 at z=0^+. (In general the singularity in the solution is two orders higher than that in the source.) Note that setting A=(2K)^-1 in (<ref>) allows one to recover the Dirac δ-function source of (<ref>) in the limit K → 0. The solution (<ref>) then tends to the solution (<ref>).



§ BEYOND LINEAR USE OF EQUATION SUPERPOSITION

Since (<ref>)  simultaneously solves the linear and the nonlinear DE with the same delta source, we attempt to use its convolution with an arbitrary source beyond the domain of the linear DE and therefore term it a BLUES (Beyond Linear Use of Equation Superposition) function B(z).
We now test to what extent the superposition

    U^(ℓ)_ϕ(z) = ∫_-∞^∞ B(z-z_0) ϕ (z_0) dz_0 ≡ (B∗ϕ)(z)

can be a useful approximation to the (unknown) solution of the nonlinear DE with the same source function ϕ (z),

    N_z U ≡ -c(k)k∂ U/∂ z + k^2 U ∂ U/∂ z- k U(1-U) - k∂^2 U/∂ z^2=  ϕ

Subtracting the nonlinear DE with operator N_z from the linear one allows one to define a residual function R_z U(z), with R_z the residual operator whose action is defined through

    R_z U ≡ L_z U -  N_z U = k(1-U) (k∂ U/∂ z + U)

Note that for an arbitrary constant trial concentration U =  O(1) of order unity the residual is of order k. 
As long as R_z U(z) is everywhere small (compared to k) the error involved in replacing the true solution by the approximation U^(ℓ)_ϕ(z) can considered to be small. Note that the residual is zero in the BLUES function, R_z B=0. 

Clearly, the residual must vanish at the boundaries in order for the solution to satisfy the boundary conditions. This implies that, by taking z → - ∞ in (<ref>), that the source cannot be arbitrary, but must be normalized as follows, 

    ∫_-∞^∞ ϕ (z) dz = 1,

in order for the superposition (<ref>) to be of any use. It is conspicious from (<ref>) that for the proper choice A=(2K)^-1 in (<ref>) the residual is small for small K. Therefore our method allows one to obtain a good approximation to the solution of the nonlinear DE with a strongly localized, but not necessarily infinitely sharply localized, source. Let us quantify this further.

The residual function for the normalized exponential corner source (with A=(2K)^-1 )
is found to be

    R_z U^(ℓ)_ϕ(z) =
    Kk/2(K+k) e^z/K ( 1-1/2e^z/K ),     z<0 
    k/2 (1- 1/K-k (K/2 e^-z/K -  k^2/K+k e^-z/k ) ) e^-z/K,     z ≥ 0,

Inspection of this function for small K shows that its magnitude is of order K, so that the approximation is good for K ≪ k.

A measure of the significance of the error can be obtained by integrating the residual over space. This leads to the residual functional R̂[U]. Since the integral of the exponential tail in (<ref>) is k, we can consider R̂[U] to be small when its value is small compared to k^2. For our example we obtain the value

    R̂[U^(ℓ)_ϕ]= K^2k(3K+4k)/4(K+k)^2

For small K (strongly localized source) this is of order K^2, so that the error is of order (K/k)^2 compared to unity. 



§ CONSTRUCTING THE EXACT SOURCE

We now move on to a different perspective, which is a non-perturbative one. Inspection shows that the residual function can be used to redefine the source  so that the solution U^(ℓ)_ϕ(z) of the linear DE with source ϕ(z) solves the nonlinear DE with source ψ (z) ≡ϕ (z) - R_z U^(ℓ)_ϕ(z). This correspondence can be written in the form

    U^(ℓ)_ϕ(z) = U^(nonℓ)_ψ(z),

with U^(nonℓ) the solution to the nonlinear DE.
We conclude that solutions of the linear DE with an arbitrary source, constructed using the Green's function method, also solve the nonlinear DE exactly for a different, but related, source. The newly constructed source makes up for the drawback that a nonlinear DE does not allow one to apply superposition of solutions. To distinguish the two sources, we denote the function ϕ henceforth by resource and the function ψ by source. The resource is used to construct both the solution (by superposition) and the source (by using the residual operator) of the nonlinear DE.

In this approach it is neither necessary that the source be strongly localized, nor that the residual function be small. We illustrate this with the following example. We let the decay constant K of the source function (<ref>) tend to the decay constant k of the traveling wavefront, and choose A=(2k)^-1, so that the resource function is

    ϕ (z)=
    1/2k e^z/k,     z<0 
    1/2k e^-z/k,     z ≥ 0

This simplifies our system to a one-parameter problem with a resource that is localized about z=0 and spread out to the same extent as the wavefront. Using (<ref>) and taking the limit K → k we obtain the following exact solution 

    U^(ℓ)_ϕ(z) = U^(nonℓ)_ψ(z)=
    
    1-1/4 e^z/k,     z<0 
     (3/4 + z/2k ) e^-z/k,     z ≥ 0,

and the exact source (with respect to the nonlinear DE) in the limit K → k satisfies

    ψ (z) =ϕ (z)-  R_z U^(ℓ)_ϕ(z) =
     (1/2k- k/4 )  e^z/k + k/8e^2z/k,     z<0 
     (1/2k- k/2 ) e^-z/k +  (3k/8 + z/4 ) e^-2z/k,     z ≥ 0,

which is well localized about z=0. Note that some restrictions on k are in order to ensure positivity of the source. 



§ AN EXAMPLE WITH A GRADIENT-SQUARED NONLINEARITY

In order to diversify somewhat our cases we present an example with a gradient-squared nonlinearity, very common in physics problems. The nonlinear partial DE with delta source,

    ∂ u/∂ t  - u(1-u) - k^2/2(∂ u/∂ x)^2 - (1 + k^2/2 u) ∂^2 u/∂ x^2 =  (k/2 +1/k)  δ(x-c(k) t),

contains a gradient-squared term and also possesses an exact wavefront solution corresponding to the simple sharp tail (<ref>) and again with the same velocity dispersion (<ref>).
We recast this partial DE into the ordinary form

    -(k+1/k) ∂ U/∂ z  - U(1-U) - k^2/2(∂ U/∂ z)^2 - (1 + k^2/2 U) ∂^2 U/∂ z^2 =  (k/2 +1/k)  δ(z),

and point out that the related linear DE, obtained by replacing (∂ U/∂ z)^2 by U(∂^2 U/∂ z^2) and then approximating U by unity but keeping all the derivatives intact, 

    -(k+1/k) ∂ U/∂ z  - (1 + k^2) ∂^2 U/∂ z^2 =  (k/2 +1/k)  δ(z),

is solved by the same function (<ref>), which is therefore a BLUES function. In this case the residual operator is

    R_z U ≡ L_z U -  N_z U =  2k/k^2+2 (U(1-U)  - k^2 (1-U/2) ∂^2 U/∂ z^2   + k^2/2(∂ U/∂ z)^2 )

The further calculational development of this example is a simple exercise along the lines of the previous sections.



§ PERTURBATION THEORY VIEWPOINT

In this section we return to situations in which the source is strongly localized and the residual function is small, so that the solution does not differ much from the solution with a delta source, except in a small neighbourhood of the center of the source. We speculate on a formal perturbation expansion and identify our method as the leading (zeroth-order) implementation of it. 
Assume that the function B solves the nonlinear DE with operator N_z and delta source,

    N_z B(z-z_0) = δ (z-z_0),

and that B is at the same time a Green's function of the linear DE with operator L_z,

    L_z B(z-z_0) = δ (z-z_0),

then we call B a BLUES function. It has the property that its residual vanishes,

    R_z B(z) =  L_z B(z) -  N_z B(z) = 0


Recalling the solution to the linear DE with an arbitrary source ϕ,

    U^(ℓ)_ϕ(z) =  (B∗ϕ)(z),

we know from the previous sections that we can view ϕ as a resource with respect to the nonlinear DE and that the superposition U^(ℓ)_ϕ(z) solves the nonlinear DE with the source

    N_z  (B∗ϕ)(z) = ϕ (z) -  R_z U^(ℓ)_ϕ(z)  ≡ϕ -  R (B∗ϕ) ≡ψ

Now we turn the problem around and ask how, for an arbitrary source ψ (z), we can obtain a useful approximation to the function U^(nonℓ)_ψ(z) that solves the nonlinear DE.

Since B∗ϕ solves the nonlinear DE, we attempt to express this function in terms of ψ. To this end we write

    ϕ  = ψ  +  R (B∗ϕ),

and obtain

    U^(nonℓ)_ψ(z) = (B∗ϕ)(z) = ∫_-∞^∞  dz_0 B(z-z_0)  ( ψ (z_0)   +  R_z_0 (B∗ϕ)(z_0)   )

If we were to continue using (<ref>) iteratively, we would arrive at the following conjectural expansion in powers of the residual, 

    (B∗ϕ)(z) = ∫_-∞^∞  dz_0 B(z-z_0)  ψ (z_0)   + ∫_-∞^∞  dz_0 B(z-z_0)  R_z_0∫_-∞^∞  dz_1 B(z_0-z_1)ψ (z_1)  
       +   ∫_-∞^∞  dz_0 B(z-z_0)  R_z_0∫_-∞^∞  dz_1 B(z_0-z_1) R_z_1∫_-∞^∞  dz_2 B(z_1-z_2)ψ (z_2) +   ...

which in compact notation reads,

    U^(nonℓ)_ψ = B∗ϕ  = ∑_n=0^∞ B∗ [  R B ]^n ψ

The first term in this series is the zeroth order term in the residual and is identical to the superposition we have proposed based on the BLUES function. It is useful provided the source is properly normalized. The second term, of first order in the residual, is finite provided we use the properly normalized source. However, the addition of this term is likely to prevent the solution from satisfying the required boundary conditions. Favorable situations arise whenever the BLUES function is itself localized and decays to zero for |z| →∞. In such case (assuming the source is localized, too) all the terms in the expansion decay to zero for |z| →∞. For more general boundary conditions (such as ours), the perturbation theory viewpoint is still incomplete, but the leading term is well defined. In the next section we test its usefulness.



§ COMPARISON WITH THE NUMERICAL SOLUTION

 
In this section we compare the approximate (zeroth-order) solution B∗ψ with the numerical  solution of the nonlinear DE with some source ψ. We consider the nonlinear DE (<ref>), which we rewrite as

    k(k U -k-1/k) ∂ U/∂ z-  kU(1-U) - k∂^2 U/∂ z^2=  ψ

with BLUES function (<ref>) and we first choose the properly normalized exponential corner source

    ψ (z) =  1/2Ke^-|z|/K

The zeroth-order solution is already given in (<ref>) for arbitrary amplitude A, so we obtain

    (B∗ψ)(z)=
    
    1-K/2(K+k) e^z/K,     z<0 
    1/2(K-k) (K e^-z/K -  2k^2/K+k e^-z/k ),     z ≥ 0

In Figs. 1 and 2 we compare this piecewise analytic approximation (red dashed line) to the numerically exact solution (black solid line) of the nonlinear DE. We also show the BLUES function B(z) (blue dotted line). For K/k ≪ 1 the approximation accurately follows the numerical solution. Even for K/k=1/4 the agreement is still very good, as Fig.1 shows. For larger values of K/k the deviations become apparent, as Fig.2 shows for K/k = 1/2. Broadening the source by increasing K/k further causes the numerical solution to change sign at some value of z, to become non-monotonic and approach zero from below for large z. A source with this broadness is unphysical in our context because the concentration variable U must be positive. On the other hand, the approximation B∗ψ given in (<ref>) is by construction always positive. It appears to be a useful approximation in the entire domain of physically acceptable sources. 

Note that for constructing the approximation it is not necessary to know the related linear DE for which B(z) is the Green's function. The quality of the approximation depends on the sharpness of the source (quantified in our example by the number K/k). However, if the related linear DE is available, it becomes possible to compare the chosen source to the “exact" source, which allows one to give an additional visual appreciation of the quality of the approximation.




In Figs. 3 and 4 we include the relevant sources together with the piecewise analytic approximation  B∗ψ. We recall that  B∗ψ exactly solves the related linear DE 

    - ∂ U/∂ z - k∂^2 U/∂ z^2=  ψ,

while the same function B∗ψ exactly solves the nonlinear DE (<ref>) with the different, but related source χ, which follows from subtracting (<ref>) from (<ref>),

    χ (z)= ψ(z) -  R_z (B∗ϕ)(z) =
    1/2Ke^z/K-Kk/2(K+k) e^z/K ( 1-1/2e^z/K ),     z<0 
    1/2Ke^-z/K-k/2 (1- 1/K-k (K/2 e^-z/K -  k^2/K+k e^-z/k ) ) e^-z/K,     z ≥ 0

The sources ψ given in (<ref>) and χ given in (<ref>) are illustrated in Figs. 3 and 4. Clearly, for small K/k (sharp source) the two sources nearly coincide and it is therefore no surprise that the approximation we propose is accurate. Increasing K/k causes the two sources to differ more and this has consequences for the quality of the approximate solution to the nonlinear DE.



Next, we demonstrate the method for the same DEs but with a step (piecewise constant) source, properly normalized in accord with (<ref>),

    ψ (z) =  
    1/σ,     -σ/2≤ z≤σ/2
    
    0,     elsewhere

This choice of source has interesting advantages, as we will illustrate.
The zeroth-order solution is readily calculated,

    (B∗ψ)(z) =
    
    1,    z < -σ/2
    k/σ (1 -  e^-z+σ/2/k )+ 1/σ(σ/2-z),     -σ/2≤ z ≤σ/2
    2k/σsinhσ/2k  e^-z/k,     z > σ/2

This approximation is particularly useful because outside the support of the source, for |z| > σ/2, it coincides with the functional form of the exact solution to the nonlinear DE without source, and respects the boundary conditions. Let us clarify this. For z < -σ/2, (B∗ψ)(z) equals unity, which trivially solves the nonlinear DE with boundary condition U(z→ -∞) → 1. For z > σ/2, (B∗ψ)(z) is of the form C  exp(-z/k), which solves the nonlinear DE exactly in this domain. The constant C>0 represents the translational invariance of the tail and embodies a shift s through C =  exp(s/k). The value of C is fixed by the requirement of continuity of the solution at the matching point z=σ/2 where the nontrivial inner solution in z ∈ [-σ/2,σ/2] joins the tail. Clearly, the asymptotic boundary condition U(z→∞) → 0 is also satisfied.

In Fig. 5 we compare this piecewise analytic approximation (red dashed line) to the numerical solution (black solid line) of the nonlinear DE. We also show the BLUES function B(z) (blue dotted line). Note that for z> σ/2, (B∗ψ)(z) is simply a shifted exponential tail with the same decay constant k as the BLUES function. Also the numerical solution is of this simple form in this domain, but with a slightly different shift. For σ/k ≪ 1 the approximation  accurately follows the numerical solution. Even for σ/k=1/2 the agreement is still very good, as Fig.5 shows. For larger values of σ/k the deviations become larger.  For this choice of step source the approximation B∗ψ given in (<ref>) as well as the numerical solution are always positive. 

In order to quantify the difference of the numerical solution and the piecewise analytic one, we can expand both in powers of ϵ≡ (z+σ/2)/k about the left edge of the step source, z=-σ/2, and for ϵ >0. Expanding (<ref>) to fourth order we obtain

    (B∗ψ)(z) =
    1-k/σ (ϵ^2/2 - ϵ^3/3! + ϵ^4/4! +  O(ϵ^5) ),        -σ/2≤ z ≤σ/2,

whereas the readily calculated expansion of the actual solution to the nonlinear DE with step source (<ref>) reads

    U^(nonℓ)_ψ (z) =
    1-k/σ (ϵ^2/2 - ϵ^3/3! + (1+k^2)ϵ^4/4! +  O(ϵ^5) ),        -σ/2≤ z ≤σ/2

Clearly, the difference starts to appear only from the fourth-order term and is such that the actual solution lies slightly below the piecewise analytic approximation. 

Likewise, we can expand in powers of μ≡ (z-σ/2)/k about the right edge of the step source, at z=σ/2, and for μ >0. We again carry the expansions through to four relevant terms. Expanding (<ref>) to third order we obtain

    (B∗ψ)(z) =
    k/σ ( (1-e^-σ/k) +  (e^-σ/k-1) μ -  e^-σ/kμ^2/2 +  e^-σ/kμ^3/3! +   O(μ^4) ),        -σ/2≤ z ≤σ/2,

whereas the readily calculated expansion of the actual solution to the nonlinear DE with step source (<ref>) reads, with the notation U^+≡ U^(nonℓ)_ψ (σ/2),

    U^(nonℓ)_ψ (z) =
    U^+ -U^+ μ + (U^+ - k/σ) μ^2/2 +  (k/σ(1+k^2) -(1+k^3/σ ) U^+  )  μ^3/3! +  O(μ^4),     -σ/2≤ z ≤σ/2

The value of U^+ can be calculated perturbatively by requiring that, in each order of the expansion, U^(nonℓ)_ψ (-σ/2) =1. This leads to the following sequence of approximations,

    U^+ = (1+σ/k)^-1,     ,


    U^+ = (1+σ/2k) (1+σ/k+ σ^2/2k^2)^-1,     ,

and

    U^+ =  (1+σ/2k-(1+k^2)σ^2/6k^2)  (1+σ/k+ σ^2/2k^2-(1+k^3/σ)σ^3/6k^3)^-1,

Incidentally, the fourth-order result is less practical because it involves a quadratic equation in U^+. Inspection of the expansions (<ref>) and (<ref>) shows that i) the value of U^+ is in general different from the approximation k/σ (1-e^-σ/k), and ii) the expansions are formally the same up to and including the term of second order in μ and would be identical to that order if we would work with equal values for the leading term U^+. 




We now return to the properties of the approximate solution (<ref>). The residual of this approximate solution  is readily calculated and is everywhere positive (which easily follows using the inequality exp(1-x) ≥ 1-x).
Next, we include the relevant sources  in the figures and recall that the piecewise analytic function B∗ψ exactly solves the related linear DE  (<ref>)
while the same function B∗ψ exactly solves the nonlinear DE (<ref>) with a different, but related source χ, which is readily calculated, 

    χ (z)= ψ (z)  - R_z (B∗ψ)(z) =
    1/σ-k/σ (σ/2 -z)   ( 1/σ (σ/2 +z)- 
    1/σ [ 1 - e^-z + σ/2/k ]   ),     -σ/2≤ z ≤σ/2
    
    0,     elsewhere

The sources ψ given in (<ref>) and χ given in (<ref>) are illustrated in Figs. 6 and 7 for two different values of σ and for k=1. Clearly, for small σ/k (sharp source) the two sources nearly coincide and our approximation is accurate. Increasing σ/k causes the two sources to differ more.






§ CONCLUSION AND OUTLOOK

We conclude that an interesting computational opportunity arises if one knows the solution B of a nonlinear DE with a delta-function source and if it coincides with the Green's function of a related linear DE. In this case the Green's function superposition that solves the linear DE for an arbitrary source also solves the  nonlinear DE for a different, but related, source. Alternatively, an approximate analytic solution to the nonlinear DE can be constructed to zeroth order in perturbation theory.  In this case the source is arbitrary apart from its normalization, which must be imposed in accord with the boundary conditions. 
This alternative use is quite general, as the related linear DE need not be known. The sharpness of the source is the key factor determining the accuracy of the approximation. Computations to higher order in the perturbation theory appear feasible for cases in which B itself is localized, with sufficiently rapid asymptotic decay.

The usefulness of the method has been illustrated in several examples using step sources or exponential corner sources in the context of nonlinear reaction-diffusion-convection equations. Since we have used superposition of equations involving Green's functions but for nonlinear operators that do not commute with the integration, we have used superposition beyond the linear theory, whence the acronym BLUES (“beyond-linear-use-of-equation-superposition") for the function B.

Nonlinear differential equations for which this approach is potentially useful occur in a variety of problems in physics, applied physics and other disciplines. Besides the biophysical example which we discussed, the following come to the speculative mind. In cosmology the field equations in Einstein's theory of general relativity are nonlinear differential equations that relate the metric and curvature of space-time to the material mass that acts as a source term (stress-energy tensor). Traveling wave solutions (gravitational waves) occur naturally in the linearized equations and persist also in the nonlinear equations, needed to describe the powerful cosmic events observed recently by LIGO <cit.>. In electromagnetism the charge density features as the source term in Maxwell's equations. Nonlinear optics takes into account higher-order terms in electric or magnetic fields necessary for describing accurately the properties of metamaterials with (near-)zero refractive index <cit.>. For these systems the perturbative point of view of nonlinear optics is currently challenged <cit.>.

More examples of relevant and active research areas in which the method may be useful, can be given. In particular, our interest also goes to condensed matter systems with propagating liquid or solid fronts, relevant mainly to physics and materials engineering. These systems differ from the active matter example treated in this paper, because in condensed matter there are conservation laws (mainly mass conservation) which constrain to an important extent the form of the applicable differential equations. Also the role of sources and sinks is now different and these concepts are replaced by, for example, force terms that correspond to intermittent pinning of a moving three-phase contact line <cit.>. 




§ ACKNOWLEDGEMENTS
 
JOI thanks Universiteit Stellenbosch for hospitality and KU Leuven for a travel grant, in the frame-work of the bilateral agreement.  This work is based on the research supported in part by the National Research Foundation of South Africa (Grant No. 99116). 


Frasca  M. Frasca, 
Green function method for nonlinear systems,
Mod. Phys. Lett. A 22, 1293 (2007); M. Frasca, 
Green functions and nonlinear systems: Short time expansion,
Int. J. Mod. Phys. A, 23, 299 (2008).
DMM
G. De Magistris and D. Marenduzzo, 
An introduction to the physics of active matter,
Physica A 34, 65 (2015)
Fisher R.A. Fisher,
The advance of advantageous genes,
Ann. of Eugenics 7, 355 (1937)
Murray  J.D. Murray,  Mathematical biology: I. An Introduction, Third Edition,  Springer (2002)

Mishra2010
A. Mishra, and R. Kumar,
Exact solutions of variable coefficient nonlinear diffusion–reaction equations
with a nonlinear convective term,
Phys. Lett. A 374, 2921 (2010)

Mishra2012
A. Mishra, and R. Kumar,
Memory effects in Fisher equation with nonlinear convection term,
Phys. Lett. A 376, 1833 (2012)

IndS J.O. Indekeu and R. Smets, Traveling wavefront solutions to nonlinear reaction-diffusion-convection equations, J. Phys. A: Math. Theor. 50, 315601 (2017)
LIGO B. P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), GW150914: The advanced LIGO detectors in the era of first discoveries, Phys. Rev. Lett. 116, 131103 (2016)
Mazur Y. Li, S. Kita, Ph. Munoz, O. Reshef, D. I. Vulis, M. Yin, M. Loncar, and E. Mazur, On-chip zero-index metamaterials, Nature Photonics 9, 738–742 (2015)
Boyd O. Reshef, E. Giese, M. Z. Alam, I. De Leon, J. Upham, and R. W. Boyd, Beyond the perturbative description of the nonlinear optical response of low-index materials, Opt. Lett. 42, 3225 (2017)
review D. Bonn, J. Eggers, J.O. Indekeu, J. Meunier, E. Rolley, Wetting and spreading, Rev. Mod. Phys. 81, 739 (2009)

