Journal of Nonlinear Mathematical Physics

Volume 28, Issue 2, June 2021, Pages 171 - 181

Exact Solutions of a Nonlinear Diffusion Equation with Absorption and Production

Robert Conte*, ORCID
Université Paris-Saclay, ENS Paris-Saclay, CNRS, Centre Borelli, Cachan F-94235, France
Department of Mathematics, The University of Hong Kong, Pokfulam Road, Hong Kong
Corresponding Author
Robert Conte
Received 29 June 2020, Accepted 19 August 2020, Available Online 10 December 2020.
10.2991/jnmp.k.200922.008How to use a DOI?
Nonlinear diffusion equation; turbulent kinetic energy; exact solutions

We provide closed form solutions for an equation which describes the transport of turbulent kinetic energy in the framework of a turbulence model with a single equation.

© 2020 The Author. Published by Atlantis Press B.V.
Open Access
This is an open access article distributed under the CC BY-NC 4.0 license (


The Partial Differential Equation (PDE) [10]

xreal,t>0:vt+(vm)xx+cvpdvq=0,m>1,p>0,c>0,q>0,d>0,v(x,0)=v0(x)0withcompactsupport, (1.1)
governs the transport of turbulent kinetic energy k in the framework of a turbulence model with a single equation. This minimal model retains the essential physical ingredients: time evolution (t), space diffusion ( x2), one space dimension only, nonlinearity (m ≠ 1), dissipation (cv p) and production (dv q). Our goal in the present paper is to find closed form solutions to serve as a validation test for the numerical schemes. In order to achieve that, we will not design new methods but apply existing methods to generate several new solutions of physical interest.

The so-called “model equation” [10] corresponds to m = p = 3/2, c = d = 1, q = 1/2, and the very special case c = d = 0 (of no interest to us) is known as the porous media equation [2,12].

In order to search for closed form solutions, it is better to first convert (1.1) to an algebraic equation,

vm1=w,wtmwwxxmm1wx2+(m1)cwp+m2m1(m1)dwq+m2m1=0, (1.2)
with the choice p = m, q = 2 − m, m arbitrary (since it contains the case of physical interest p = 3/2, q = 1/2 when m = 3/2).

Therefore the final algebraic equation which we will study is (after rescaling t) [10, Eq. (3.24)]

λvm1=w:E(w)wtbwwxxwx2+cw2d=0,b>0,c>0,d>0. (1.3)

Since there could exist physical systems not requiring the positivity of b, c, d, we will also mention a few solutions with b, c, d not all positive.

The paper is organized as follows. In Section 2 we first recall the ingenious, however not generalizable, method which has allowed Maire to obtain a solution matching all the physical constraints in (1.3).

Section 3 is devoted, using the method of infinitesimal Lie point symmetries, to the construction of all the reductions of the PDE (1.3) to an Ordinary Differential Equation (ODE), then to their integration.

In Section 4, we study the local behaviour of the field w(x, t) near its movable singularities, a prerequisite to the search for closed form solutions.

The next Section 5 is devoted to a search for new exact solutions based on the singularity structure.


The PDE (1.3) belongs to the class

wtP(x,w)=0,Pdifferentialpolynomial, (2.1)
with the additional property of existence of a finite number of functions fk(x) whose linear combinations are stable under P(x, w),
(α1,,αK)(β1,,βK):P(x,k=1Kαkfk(x))=k=1Kβkfk(x). (2.2)

Then Galaktionov and Posashkov [7] observed that the (kind of “adiabatic”) assumption

w(x,t)=k=1Kγk(t)fk(x), (2.3)
amounts to solving a nonlinear system of ODEs (no more PDEs) for the functions γk(t). In the present case (1.3),
one thus obtains four solutions leaving the physical parameters b, c, d unconstrained and positive.

The most physically relevant solution [10, Eq. (3.37)],

{λvm1=w=Acosh(ωt)cosh(kx)sinh(ωt),ω=Ac2+b1+b,k2=c1+b,A2=dc, (2.4)
does not depend on any movable (i.e. function of the initial conditions) constant, excluding of course the arbitrary origins x0 and t0, and it provides a good description [10] of the transport of the turbulent kinetic energy.

The second solution is stationary,

w=Acosh(kx),k2=c1+b,A2=(1+b)dc, (2.5)
and the third one homogeneous in space
w=ω22ctanhω22t,ω 22=4cd. (2.6)

Finally, the fourth solution is characterized by a second order nonlinear ODE,

{w=g0(t)+g1(t)cosh(kx),k2=c1+b,g 12=g 02dc0,(b+1)g02(2b+3)cg0g02(b+2)cg0(cg 02d)=0, (2.7)
and, if one excludes the particular values g0 listed in the three previous solutions, its only physical (b positive) solutions are multivalued [3, Eq. (9.2)], but nevertheless expressible via quadratures [9]. The same kind of ODE will again be encountered hereafter, see Eq. (5.17) for F1.

If one relaxes the constraint b > 0, there also exists a fifth solution,

w=ω22ctanhω22t+Acos(kx),ω 22=4cd+4c2A2,k2=c,b=2,A=arbitrary. (2.8)

Since the above ingenious method only applies to the restricted class (2.1) and (2.2), let us in addition apply the two main classes of methods able to find explicit solutions of algebraic PDEs:

  • Those based on the symmetries of the PDE, which generate reductions to ODEs;

  • Those based on the movable singularities of the PDE, which (after a double study, local then global) generate closed form solutions w(x, t).


For generic values of b, c, d, the only Lie point symmetries of the PDE (1.3) are arbitrary translations of both x and t. The resulting characteristic system,

dxα=dtβ=dw0,α,βarbitraryconstants, (3.1)
admits the two invariants w and β xα t, with α, β not both zero, thus defining the unique reduction to an ODE,
w(x,t)=W(ξ),ξ=βxαt,αWβ2(bWW+W2)+cW2d=0. (3.2)

One must distinguish β ≠ 0 (the reduction preserves the differential order two) and β = 0 (the reduction lowers it to one).

For β ≠ 0, since the positive parameter b is never equal to −(1 − 1/n), with n some signed integer, the general solution W(ξ) is multivalued [11] and generically cannot be obtained in closed form. A notable exception is β ≠ 0 and α = 0 (stationary wave),

ξ=x,w(x,t)=W(x),bWWW2+cW2d=0, (3.3)
which admits the first integral [10]
K=[(1+b)(W2d)+cW2]W2b. (3.4)

For K = 0, the general solution is physically acceptable and has already been obtained, see Eq. (2.5), and for K ≠ 0 the solution is given implicitly by the quadrature

w(x,t)=W(x),x=x0+(1+b(1+b)dcW2KW2b)12dW. (3.5)

The only cases of invertibility of this quadrature are −2/b = 1, 2, 3, 4, yielding expressions W(x) trigonometric (b = −2, −1) or elliptic (b = −2/3, −1/2), which however all violate the physical requirement b > 0.

For β = 0 one obtains the front independent of x, Eq. (2.6).


Any closed form solution depends on arbitrary functions or constants, such as x0, t0 in (2.4), which may define movable singularities. For instance, the solution (2.4) definies two families of movable singularities: on one side the movable poles of w located at the points t=t0+niπ/ω,n, on the other side movable poles of 1/w (movable zeroes of w) located on the singular manifold φ (x, t) = 0 defined by

φ(x,t)coshω(tt0)coshk(xx0)=0. (4.1)

A prerequisite to the systematic search for solutions is therefore the determination of the structure of the movable singularities of (1.3),

wφ0w0φp,p𝒩. (4.2)

Since the highest derivative term wwxx displays the singularity w = 0, one must also study the movable zeros of w (movable poles of w−1 = f ),

λvm1=w=f1:E(f)f2ft+bffxx(2b+1)f x2df4+cf2=0. (4.3)

This is a classical computation [6, §4.4.1], whose results are the following.

The PDE (1.3) admits two types of movable singularities.

  1. (1)

    If φx = 0 (the singular manifold is then said characteristic), the highest derivative term wwxx does not contribute to the leading order, and w (as well as f ) presents one family of movable simple poles,

    φx=0:wφ0c1φtφ1, (4.4)
    φx=0:w1φ0d1φtφ1. (4.5)

  2. (2)

    If φx ≠ 0 (noncharacteristic singular manifold), w has no movable poles and w−1 = f presents two families of movable simple poles,

    fφ0f0φxφ1,df 02φtφxf0+1=0,f00. (4.6)

To finish this local analysis, one must then compute the Fuchs indices of the linearized equation of (1.3) in the neighborhood of φ = 0. Indeed, a necessary condition of singlevaluedness is that all Fuchs indices be integers of any sign. For the families (4.4) or (4.5), the unique Fuchs index is i = −1. For each of the two families (4.6), one finds

i=1,df 021b, (4.7)
and the noninteger value (even nonrational) of the second index reflects the high level of nonintegrability of the initial PDE.

In order to build solutions of this kind of nonintegrable PDE, the various methods based on the singularity structure are reviewed in summer school lecture notes [5], where the proper original references can be found. Let us now investigate a few of them.


They consist in requiring the Laurent series of a single family to terminate. For the respective local behaviours (4.4)(4.6), the corresponding possible solutions are defined by

(φx=0):w=c1χ1+w1, (5.1)
(φx=0):f=d1χ1+f1, (5.2)
(φx0):f=f0χ1+f1,f00, (5.3)
in which the expansion variable χ (x, t) is any homographic transform of φ (x, t) vanishing with φ. There exists an optimal choice of χ [4], characterized by its gradient ((ξ, η) represents either (x, t) or (t, x))
χξ=1+S2χ2,χη=C+Cξχ12(CS+Cξξ)χ2, (5.4)
and the constraint
Sη+Cξξξ+2CξS+CSξ=0. (5.5)

After substitution in Eqs. (1.3) and (4.3), the LHS E of these equations become Laurent series in χ which also terminate,

Ej=0qEjχj+q,q, (5.6)
and the coefficients Ej only depend on w0, w1 (or f0, f1) and (S, C). The solutions are then provided by solving the determining equations
j=0,...,q:Ej(w0,w1,S,C)=0. (5.7)

5.1. Characteristic One-family Truncation of w

The truncation (5.1) with χx = 0 defines the system

{E00,E1bc1w1,xx+2w1=0,E2bw1w1,xxw1,x2+w1,t+cw 12S2cd=0, (5.8)
whose general solution is
w1=Acos(kx),S=2cd2c2A2,K2=c,(b+2)A=0. (5.9)

The value of χ results by integrating the Riccati system (5.4),

χ1=ω22tanhω22t,ω 22=4cd+4c2A2. (5.10)

The bifurcation (b + 2)A = 0 defines two previously found solutions w, (2.6) and (2.8).

5.2. Characteristic One-family Truncation of f

The truncation (5.2) with χt = 0 defines the system

{E00,E1f1=0,E2cd+S25d2f 12df1,t=0,E32cdf1+Sf14d2f 132df1f1,t+bdf1,xx=0,E4E(f1)+S2df 12=0, (5.11)
whose general solution is
f1=0,S=2cd. (5.12)

After integration of the Riccati system (5.4),

χ1=ω22tanhω22t,ω 22=4cd, (5.13)
one obtains
f=d1ω22tanhω22t,ω 22=4cd, (5.14)
identical to (2.6).

To conclude, these characteristic truncations yield nothing new.

5.3. Noncharacteristic One-family Truncation of f

Let us finally consider one of the two families (4.6), whose Fuchs indices are (4.7), and let us assume bd ≠ 0. The truncation (5.3) defines the system

{E0df 02+Cf0+1=0,E12(2df 02Cf0+b)f1+f 02Cxf0f0,t+2(1+b)f0,x=0,EjEj(f0,f1,S,C)=0,j=2,3,E4b(2b+1)[f 02S+2f 12+2f1f0,x2f0f1,x]2=0,XSt+Cxxx+2CxS+CSx=0, (5.15)
and the factorization of E4 makes the resolution easy. It is even easier after the change of functions
(f0,f1)(F0,F1):f0=1/F0,f1=f0F112f0,x. (5.16)

In the case b ≠ −1/2, the algebraic (i.e. nondifferential) elimination of F0,t, F1,t, Sx yields the much simpler equivalent system,

b12:{E1F0,t(b+2)F0F0,x+2(b+1)F 02F12dF1=0,E2bF0,xx+2(b+2)F1F0,x+2[(3b+2)F1,x(4(b+1)F 12+c)]F02F1,t=0,XbF1,xx2(3b+2)F1F1,x+4(b+1)F 13cF1=0,S=F0,xxF032F 0,x2F 022F 12+2F1,x2F0,xF0F1,C=F0dF 01. (5.17)

This system (5.17) is solved in three steps:

  1. (i)

    integration of the ODE in F1 defined by X = 0, which introduces at most two arbitrary functions of t;

  2. (ii)

    determination of F0 by solving the overdetermined system (E1 = 0, E2 = 0);

  3. (iii)

    knowing the values of (S, C), integration of the Riccati system (5.4) for χ.

In the unphysical case b = −1/2, we could not find an equivalent system as simple as (5.17).

Let us now perform the above mentioned three steps.

5.3.1. Values of F1

As already mentioned about the similar ODE (2.7), we will discard the multivalued solutions of the ODE X = 0 for F1(t), refering the interested reader to Ref. [9].

All singlevalued solutions of the ODE for F1 are obtained by two methods: for the general solution, by looking in the classical exhaustive tables [8,11]; for particular solutions, by looking for Darboux polynomials. One thus finds exactly seven solutions, five of them with a negative b (unphysical for the diffusion problem, but possibly admissible for other systems) and two with an arbitrary value of b,

b=43:F1=xlog[cosh(kxg(t))cosh(K(t))],k2=34c, (5.18)
b=45:F1=xlog[(xg(t))e0],e0=5c48,g2=12e 02,g3(t), (5.19)
b=12:F1=12xlog[(xg(t))e0],e0=c6,g2=12e 02,g3(t), (5.20)
b=23:F1=(xg(t))e0,e0=c2,g2=12e 024K(t),g3=8e 03+4e0K(t), (5.21)
b+10:F1=k2tanhk2(xg(t)),k2=cb+1, (5.22)
b+10:F1=b2(b+1)k2tanhk2(xg(t)),k2=4(b+1)cb2, (5.23)
b=2:F1=0, (5.24)
in which (x,g2,g3) is the elliptic function of Weierstrass, and g, K two arbitrary functions of t. The first four depend on two arbitrary functions of t, the next two on one arbitrary function of t.

For b = −1/2, the obtained solution is a particular case of the solution, which we could not obtain, resulting from the system (5.15).

For b = −1 and c ≠ 0, the ODE for F1 has no singlevalued solution.

Let us next determine F0. By the elimination of F1 between E1 = 0 and E2 = 0, the value of F0 is the root of a sixth degree polynomial whose coefficients are polynomial in F1 and its derivatives. However, this computation is only tractable for the two tanh solutions (which only depend on g(t)), and one finds constant values for F0 and g′(t),

b+10:F 02=db+1,g(t)=+(b+2)k2F0, (5.25)
b+10:F 02=db+1,g(t)=(b+2)k2F0. (5.26)

In the four other cases b = −4/3, −4/5, −1/2, −2/3 of the list (5.18)(5.21), this technical difficulty is overcome by integrating the linear inhomogeneous ODE E2 = 0 for F0 as follows.

One first notices that, in its homogeneous part, the simple pole of F1 with residue r = −1 is a Fuchsian singularity for the EDO in F0, whose Fuchs indices i, the roots of

bi2(b+4r+2br)i+8(b+1)r2+2(3b+2)r=0, (5.27)
are irrational for the four values of b. Since F0 is necessarily an algebraic function of F1 and its derivatives, the only such algebraic solution of the homogeneous part of E2 = 0 is F0 = 0, see the example Eq. (5.34) hereafter. One then computes F0 as a particular solution of the inhomogeneous equation E2 = 0. Once this value of F0 obtained, the nonlinear equation E1 = 0 generates constraints on g(t) and K(t). This method is exemplified in Subsection 5.3.5.

5.3.2. Solution, case of the first tanh value of F1


b+10:F0=1a0,a 02=b+1d,g(t)=+(b+2)k2F0, (5.28)
one finds successively
S=k22,C=b+2a0, (5.29)
then, by integration of the Riccati system (5.4)
χ1=k2tanhk2(x+b+2a0t), (5.30)
the solution
w1=f=a0k2[tanhk2(x+b+2a0t)tanhk2(xb+2a0t)],k2=cb+1, (5.31)
which is just another representation of the solution (2.4) already found by the method of Galaktionov and Posashkov.

5.3.3. Solution, case of the second tanh value of F1

Solving (E1, E2) for (F0(x, t), g(t)) is again quite easy,

then one obtains

The ODE for χ−1 is a Lamé equation in its Riccati form, whose solution is singlevalued for b = −1 + 1/n, n, multivalued otherwise. For the present diffusion problem, this new solution leaves b, c, d unconstrained and can indeed be used to test the validity of numerical schemes.

5.3.4. Case b = −2

A computation similar to the above one yields

F1=0,F0=Acoshkx,S=c[132tanh2kx],C=A2+ddcoshkx,k2=c, (5.32)
χ1=xlog[cosh(kx)1/2(ΩAkcosh(Ωt)+sinh(Ωt)cosh(kx))],Ω2=(A2+d)c, (5.33)
i.e. a solution identical to (2.8) already found in Section 2.

5.3.5. Solutions, case b = −4/3

The resolution of the four cases b = −4/3, −4/5, −1/2, −2/3 follows the same pattern, so we only detail it for b = −4/3 [Eq. (5.18)].

For this value b = −4/3, the ODE E2 = 0 possesses the general solution

F0=g+(t)[coshK(t)cosh(kxg(t))+sinhK(t)sinh(kxg(t))1]2(cosh(kxg(t))coshK(t))21+g(t)[coshK(t)cosh(kxg(t))+sinhK(t)sinh(kxg(t))1]2(cosh(kxg(t))coshK(t))21+3[KsinhK(t)sinh(kxg(t))+gcoshK(t)cosh(kxg(t))g]2ksinh2K(t), (5.34)
in which g±(t) are two other arbitrary functions of t.

As already argued in Subsection 5.3.1, since the relation between F0 and ekxg(t) is necessarily algebraic, the two functions g+ and g must vanish. Equation E1 = 0 then yields the necessary and sufficient conditions,

g=0,K=0,gK=0,g2+K2=cd, (5.35)
solved as
g(t)=ω1t,K(t)=Ω1t+k0,ω1Ω1=0,ω 12+Ω 12=cd, (5.36)
in which ω1, Ω1, k0 are constant.

The system (5.17) therefore has for solution, in the first case ω1 ≠ 0, Ω1 = 0,

{F1=xlog[cosh(kxω1t)coshk0],k2=34c,ω 12=cd,F0=3ω12ksinh2k0[coshk0cosh(kxω1t)1],S=k22[13sinh2k0(coshk0cosh(kxω1t)1)2],C=ω12k3(coshk0cosh(kxω1t)1)2sinh4k0sinh2k0(coshk0cosh(kxω1t)1), (5.37)
and in the second case ω1 = 0, Ω1 ≠ 0,
{F1=xlog[coshkxcoshΩ1t],F0=3Ω1sinhkx2ksinhΩ1t,k2=34c,Ω 12=cd,S=k22(1+3sinh2kx),C=Ω12k[sinhΩ1tsinhkx3sinhkxsinhΩ1t]. (5.38)

The integration of the Riccati system (5.4)(5.4) introduces another arbitrary constant t0, then the solutions w are defined by (5.3) and (5.16). One thus obtains two solutions w, respectively

{χ1=k2a1sinhξ(coshk0sinhξ+cosh2k02)+a2(coshk0sinh2ξ2sinhξ+coshk0)[a1(sinhξcoshk0)+a2sinhξ][coshk0sinhξ1],a1a2=ω4tanh(ω4(tt0))2ω1,ξ=kxω1t,k2=34c,ω 12=cd,ω 42=(13sinh2k0)cd,w=3ω1(coshk0coshξ)2ksinh2k0[ω1ω42tanh(ω4(tt0))(coshk0coshξ)], (5.39)
{χ1=k2a1(cosh2kx2)a2coshkx[a1coshkx+a2]sinhkx,k2=34c,a1a2=3sinh(2Ω1t)6Ω1(tt0)6sinhΩ1t2sinh3(Ω1t)6Ω1(tt0)coshΩ1t),Ω 12=cd,w=3Ω12k2[coshΩ1tcoshkx][a1coshkx+a2]sinhΩ1t[a1coshΩ1t+a2]. (5.40)

Each of these two new solutions is outside the class (2.4) and depends on a single arbitrary constant (k0 in the first one, t0 in the second one).

5.3.6. Case b = −4/5

With the value (5.19) of F1, the function F0, algebraic in the derivatives of F1, is necessarily an affine function of ζ, and xg [1, §18.6],

F0=R0+R1+R2ζ+(xg(t))R3, (5.41)

whose coefficients Rj are rational in (xg(t)). One then proves that, since d is nonzero, the discriminant g 2327g 32 must vanish, thus reducing F1 and F0 to simply periodic functions,

F1=xlogψ(x,t),ψ=(ktanh(k(xg(t))))2(2/3)k2e0,k2=3ɛe0,ɛ2=1,F0=56g[1+(l31/8)e 03ψ3],(ɛ1)l3=0. (5.42)

Equation E1 = 0 then generates the constraints

ɛ=1,g2=365d, (5.43)
thus restricting this solution to the particular case b = −4/5 of (5.23).

5.3.7. Case b = −1/2

Since we could not solve the system (5.15), we only solve here its particular case (5.17). One first establishes the particular solution F0 of E2 = 0, namely, in the elliptic subcase g 2327g 320,

y=(xg(t))e0,F0=2y3+g3+8e 03y3[g33g32Δ((xg)g38e 02ζ)]3g32Δ(8e 02y2(g3+8e 03)(y2e0)),Δ=27(64e 06g 32), (5.44)

and in the trigonometric subcase g 2327g 32=0,

F1=12xlogψ(x,t),ψ=(ktanh(k(xg(t))))2(2/3)k2e0,k2=±3e0,F0=23g[1+4(1ɛ)e 03ψ3].(5.45)

Equation E1 = 0 then generates the constraints g′ = 0 in the elliptic subcase, and ε = −1, g2 = 9d/2 in the trigonometric subcase, restricting the solution to (5.23) with b = −1/2.

5.3.8. Case b = −2/3

This case also happens to be a particular case of (5.23).


By a systematic investigation, we have obtained several new solutions of this diffusion problem. Two of them [Subsection 5.3.3, Eq. (2.7)] match all the physical constraints b > 0, c > 0, d > 0 and can serve to calibrate and validate the numerical schemes.

The two other new solutions, Eqs. (5.39) and (5.40), only valid for a negative value of b, could be useful for other diffusion problems governed by (1.3) with b < 0.

Finally, in this short paper we have only investigated those solutions w which take into account one of the two movable poles. Taking account of both poles via the two-singular manifold method (see [6, §] and references therein) should certainly generate additional solutions.

Remark. As suggested by the two new solutions Eqs. (5.39) and (5.40), the class w equal to a second degree polynomial in cosh(kx) with time-dependent coefficients could also generate physically interesting solutions.


The author declares no conflicts of interest.


The author gratefully acknowledges the support of LRC Méso and is happy to thank B.-J. Gréa, A. Llor and R. Motte for suggesting this interesting and challenging problem.


[1]M Abramowitz and IA Stegun, Handbook of Mathematical Functions, tenth printing, Dover, New York, 1972.
[7]VA Galaktionov and SA Posashkov, Exact solutions and invariant spaces for nonlinear gradient diffusion equations, Comput. Math. Math. Phys., Vol. 34, 1994, pp. 313-321.
[10]PH Maire, Étude d’une équation de diffusion non-linéaire. Application à la discrétisation de l’équation d’énergie cinétique turbulente pour un modèle de turbulence à une équation, 2001, pp. 81. available from: Rapport CEA D01 03661.
Journal of Nonlinear Mathematical Physics
28 - 2
171 - 181
Publication Date
ISSN (Online)
ISSN (Print)
10.2991/jnmp.k.200922.008How to use a DOI?
© 2020 The Author. Published by Atlantis Press B.V.
Open Access
This is an open access article distributed under the CC BY-NC 4.0 license (

Cite this article

AU  - Robert Conte
PY  - 2020
DA  - 2020/12/10
TI  - Exact Solutions of a Nonlinear Diffusion Equation with Absorption and Production
JO  - Journal of Nonlinear Mathematical Physics
SP  - 171
EP  - 181
VL  - 28
IS  - 2
SN  - 1776-0852
UR  -
DO  - 10.2991/jnmp.k.200922.008
ID  - Conte2020
ER  -