Diffusion Equation

Chapter 2
Diffusion equations
2.1 Separation of variables: intervals
Diffusion equation is a linear partial differential equation, since the functions related to u in the
equations (u
t
and ∆u) are both linear. Recall for the linear ordinary differential equations:
dP
dt
= kP, and D
d
2
P
dx
2
= kP, (2.1)
it is well-known that e
kt
and Acos(

k/Dtx) + Bsin(

k/Dx) are the solutions of the respective
equations. Indeed in the second case, the exponential function e
kitx/D
is a complex solution. Thus
we can make a lucky guess for the solution of the diffusion equation:
∂P
∂t
= D
∂
2
P
∂x
2
, (2.2)
that the solution is an exponential function
P(t, x) = e
at+bx
, (2.3)
for some constants a, b, which could be complex numbers. By substitution, we find that the function
in (2.3) is a solution of (2.2) if a = Db
2
, and depending on the values of b, we obtain three families
of solutions:
P(t, x) = e
Db
2
t
e
bx
, b > 0, (2.4)
P(t, x) = e
Db
2
t
e
−bx
, b > 0, (2.5)
P(t, x) = e
−Db
2
t
e
−bix
, b > 0. (2.6)
The first two families of solutions are not reasonable solutions, since when x → ∞ (or −∞)
P(t, x) →∞, which implies an unlimited growth at x = ∞ (or −∞). The third family of solutions
are complex, but their real and imaginary parts are both solutions of (2.2) too:
P(t, x) = e
−Db
2
t
cos(bx), and P(t, x) = e
−Db
2
t
sin(bx). (2.7)
19
20 CHAPTER 2. DIFFUSION EQUATIONS
These solutions are little more “reasonable” as they are bounded as x → ∞, but still they do not
satisfy natural boundary conditions P, P
x
→0 as x →∞.
Nevertheless solutions with above forms are solutions of diffusion equations, and we notice that
they are in form of P(t, x) = U(t)V (x). This motivates us to consider the solutions which is
separable on the two variables t and x:
u(t, x) = U(t)V (x). (2.8)
With the form of u in (2.8), the diffusion equation (2.2) becomes
U
′
(t)V (x) = DU(t)V
′′
(x), and
U
′
(t)
U(t)
= D
V
′′
(x)
V (x)
= Dk,
where k is a constant independent of both t and x. Thus U satisfies an equation:
U
′
(t) = DkU(t), t > 0, (2.9)
and
U(t) = C
1
e
Dkt
; (2.10)
and V satisfies
V
′′
(x) = kV (x). (2.11)
The solutions of (2.11) are C
2
e
√
kx
with any complex number k. Thus we obtain solution of (2.2):
Ce
Dkt
e
√
kx
, (2.12)
just like the three families of solutions we obtain by making lucky guess. To obtain more specific
solutions, we must take the boundary and initial conditions into considerations.
We start with the one-dimensional diffusion equation on an interval with length L:
∂P
∂t
= D
∂
2
P
∂x
2
, x ∈ (0, L). (2.13)
Here the individuals of the the species inhabits a one dimensional patch (0, L), and the diffusion
coefficient of the population is D > 0. To get a definite solution of the problem, we need to have
an initial population distribution:
P(0, x) = P
0
(x), x ∈ (0, L), (2.14)
and a boundary condition. We assume that the exterior is hostile so homogeneous Dirichlet bound-
ary condition is satisfied:
P(t, 0) = P(t, L) = 0. (2.15)
We look for solution with form P(t, x) = U(t)V (x). Then we obtain (2.9) and (2.11), and (2.15)
implies that
V (0) = V (L) = 0. (2.16)
(2.11) and (2.16) become a boundary value problem for an ordinary differential equation, which is
a special case of Sturm-Liouville problems. We shall discuss this kind of problems in a separate
section next.
2.2. EIGENVALUES AND EIGENFUNCTIONS 21
2.2 Eigenvalues and eigenfunctions of boundary value problems
We consider the Dirichlet boundary value problem:
y
′′
= ky, x ∈ (0, L), y(0) = y(L) = 0. (2.17)
We can see that y(x) = 0 is always a solution of the equation. So we shall look for solutions other
than the trivial zero solution. From elementary ODE knowledge, we know that the general solution
of y
′′
= ky is one of the following three forms:
c
1
e
−
√
kt
+c
2
e
√
kt
, if k > 0, (2.18)
c
1
+c
2
t, if k = 0, (2.19)
c
1
cos(
√
−kt) +c
2
sin(
√
−kt), if k < 0. (2.20)
One important characteristic of a boundary value problem is that it may not have a solution, while
the initial value problem always has a solution and the solution is unique if all terms in the equation
are differentiable. For (2.17), there is no solution when k > 0. Indeed if there is a solution, it must
be y(t) = c
1
e
−
√
kt
+c
2
e
√
kt
for some constants c
1
and c
2
, and we can solve the constants by
0 = y(0) = c
1
+c
2
, 0 = y(L) = c
1
e
−
√
kL
+c
2
e
√
kL
.
But this system of equations is unsolvable: if one gets c
1
= −c
2
from the first equation and plugs
it into the second one, one has (e
−
√
kL
−e
√
kL
)c
2
= 0, and one must have c
2
= 0 and c
1
= 0. Thus
(2.17) has only zero solution when k > 0. It is also easy to check that (2.17) has only zero solution
when k = 0. So the only meaningful case is when k < 0.
In this case, (2.17) may still have only the zero solution, but for some number k < 0, it can
have a solution y(t) which is not the zero solution. (Note that from the linear principle, if y(t) is a
solution of (2.17), so is c · y(t) for any constant c, and if y
1
(t) and y
2
(t) are solutions of (2.17), so
is y
1
(t) +y
2
(t).) So for k which (2.17) has non-zero solutions, we call these numbers eigenvalues,
and the solutions of (2.17) eigenfunctions. (In fact, historically the solutions of (2.17) were first
called eigenvalues and eigenfunctions.)
Now let’s search for the eigenvalues and eigenfunctions of (2.17): the eigenfunction must be
y(t) = c
1
cos(
√
−kt) + c
2
sin(
√
−kt) for some constants c
1
and c
2
; by the boundary conditions, we
have
0 = y(0) = c
1
, 0 = y(L) = c
1
cos(
√
−kL) +c
2
sin(
√
−kL);
so c
1
= 0 and if c
2
is not zero, sin(
√
−kL) = 0; thus
√
−kL = mπ, the eigenvalue k must be
k = −
m
2
π
2
L
2
, (2.21)
and the eigenfunctions associated with k = −
m
2
π
2
L
2
are c sin

mπt
L

. Therefore (2.17) has a
sequence of eigenvalues: k
m
= −
m
2
π
2
L
2
, and the set of eigenfunctions are the sine functions which
22 CHAPTER 2. DIFFUSION EQUATIONS
can “fit in” the interval (0, L). When the interval is (0, π) (so L = π), the expression of eigen-pairs
is much simpler:
k
m
= −m
2
, f
m
(t) = sin(mt). (2.22)
Back to (2.13) and (2.15), we find a sequence of solutions:
u
m
(t, x) = exp

−D
m
2
π
2
L
2
t

sin

mπx
L

, m is a positive integer. (2.23)
Because of the linear principle, the general solution of (2.13) is
u(t, x) =
∞
¸
m=1
c
m
exp

−D
m
2
π
2
L
2
t

sin

mπx
L

, t > 0, x ∈ (0, L). (2.24)
Example 2.1.
y
′′
+ky = 0, y(0) +y
′
(0) = y(π) = 0. (2.25)
This problem is same as previous example except the boundary conditions. So we still consider
the cases of k < 0, k = 0 and k > 0. When k < 0, y(t) = c
1
e
−
√
−kt
+ c
2
e
√
−kt
, from the boundary
conditions, we get y(0)+y
′
(0) = c
1
+c
2
−
√
−kc
1
+
√
−kc
2
= 0 and y(π) = c
1
e
−
√
−kπ
+c
2
e
√
−kπ
= 0,
So (c
1
, c
2
) satisfies the equations

1 −
√
−k 1 +
√
−k
e
−
√
−kπ
e
√
−kπ

·

c
1
c
2

=

0
0

(2.26)
The equation has non-zero solution only if det

1 −
√
−k 1 +
√
−k
e
−
√
−kπ
e
√
−kπ

= 0. The determinant equals
to (1 −p)e
pπ
−(1 +p)e
−pπ
if we denote p =
√
−k, thus det = 0 is equivalent to e
2pπ
=
1 +p
1 −p
. The
functions f
1
(p) = e
2pπ
and f
2
(p) =
1 +p
1 −p
have no common points when p > 0 (in fact, they only
intersect at p = 0. thus e
2pπ
=
1 +p
1 −p
is not possible for p > 0, and there is no eigenvalue such that
k < 0. When k = 0, y(t) = c
1
+c
2
t, from the boundary conditions, we get y(0) +y
′
(0) = c
1
+c
2
= 0
and y(π) = c
1
+c
2
π = 0, since the determinant of the matrix

1 1
1 π

is not zero, then only solution
for (c
1
, c
2
) is (0, 0). So 0 is not an eigenvalue either.
For k > 0, y(t) = c
1
cos(
√
kt) + c
2
sin(
√
kt). The boundary conditions: y(0) + y
′
(0) = c
1
+
c
2
√
k = 0, y(π) = c
1
cos(
√
kπ) +c
2
sin(
√
kπ) = 0. The matrix is

1
√
k
cos(
√
kπ) sin(
√
kπ)

, and the
determinant of the matrix is 0 if sin(
√
kπ) −
√
k cos(
√
kπ) = 0, or tan(πp) = p for p =
√
k. The
graphs of f
1
(p) = tan(πp) and f
2
(p) = p intersect at a sequence of points, p
m
, m = 1, 2, 3, · · · , and
0.5 < p
1
< 1.5, 1.5 < p
2
< 2.5, m − 0.5 < p
m
< m + 0.5. k
m
= p
2
m
are the eigenvalues of the
problem (2.25), and the eigenfunctions are y
m
(t) =
√
k
m
cos(
√
k
m
t) −sin(
√
k
m
t).
2.3. FOURIER SERIES 23
2.3 Fourier series
To completely solve (2.13)-(2.15), we need to determine the constants c
m
in (2.24) by setting t = 0
in (2.24):
u
0
(x) =
∞
¸
m=1
c
m
sin

mπx
L

, x ∈ (0, L). (2.27)
The expression in the right-hand side of (2.24) is the Fourier series of u
0
(x) on (0, L). Fourier
series is widely used in mathematical physics and signal process. All sine and cosine functions are
considered to be regular smooth waves, and Fourier series convert any noise (or signal) into a sum
of more regular waves, and that is a big help for the signal processing. Here we collect a few facts
of Fourier series.
Theorem 2.2. Let f(x) and f
′
(x) be piecewise continuous functions on the interval [0, L]. Then
f(x) can be expanded in either a pure sine series
f(x) =
∞
¸
m=1
b
m
sin

mπx
L

, (2.28)
where b
m
=
2
L

L
0
f(x) sin

mπx
L

dx, m = 1, 2, · · · , or a pure cosine series
f(x) =
a
0
2
+
∞
¸
m=1
a
m
cos

mπx
L

, (2.29)
where a
m
=
2
L

L
0
f(x) cos

mπx
L

dx, m = 0, 1, 2, · · · .
The Fourier series is similar to the Taylor series
f(x) =
¸
m=0
f
(n)
(0)
n!
x
n
, (2.30)
but instead of expanding the function in polynomials, Fourier series expands a function in sine or
cosine functions. But Taylor series converges to the original function only in a small neighborhood
of x = 0 (though some time it is large), and the neighborhood depends on the function f. Fourier
series converges to the original functions for all functions and on the whole interval (0, L). (But it
may not converge to f(x) at x = 0 and x = L.)
Theorem 2.2 shows all sound wave can be decomposed into several regular waves plus a small
noise. The approximation of Fourier series is sine or cosine Fourier polynomials:
f(x) ≈
N
¸
m=1
b
m
sin

mπx
a

, f(x) ≈
a
0
2
+
N
¸
m=1
a
m
cos

mπx
a

. (2.31)
24 CHAPTER 2. DIFFUSION EQUATIONS
Now we can complete the solution formula of (2.13)-(2.15):
u(t, x) =
∞
¸
m=1
c
m
exp

−D
m
2
π
2
L
2
t

sin

mπx
L

,
where c
m
=
2
L

L
0
u
0
(x) sin

mπx
L

dx, m ≥ 1, t > 0, x ∈ (0, L).
(2.32)
Example 2.3.





u
t
= u
xx
, t > 0, x ∈ (0, π),
u(t, 0) = u(t, π) = 0,
u(0, x) = sin x.
(2.33)
From the solution formula (2.32), we get
u(t, x) =
∞
¸
m=1
c
m
e
−m
2
t
sin(mx), and c
m
=
2
π

π
0
sin xsin(mx)dx. (2.34)
However

π
0
sin xsin(mx)dx = 0 for any m ≥ 2, and

π
0
sin
2
xdx = π/2. Therefore the solution is
simply
u(t, x) = e
−t
sin x. (2.35)
Example 2.4. For the homogeneous Neumann boundary value problem





u
t
= Du
xx
, t > 0, x ∈ (0, L),
u
x
(t, 0) = u
x
(t, L) = 0,
u(0, x) = u
0
(x),
(2.36)
the solution is
u(t, x) =
c
0
2
+
∞
¸
m=1
c
m
exp

−D
m
2
π
2
L
2
t

cos

mπx
L

,
where c
m
=
2
L

L
0
u
0
(x) cos

mπx
L

dx, m ≥ 0, t > 0, x ∈ (0, L).
(2.37)
The derivation of this formula is left as exercise.
For these examples, the asymptotic behavior of the solution can be read from the formulas. For
homogeneous Dirichlet problem (2.13)-(2.15), the asymptotic limit is
lim
t→∞
u(t, x) = 0, and the asymptotic profile is u(t, x) ∼ c
1
exp

−
π
2
t
L
2

sin

πx
L

. (2.38)
On the other hand, for homogeneous Neumann problem (2.36), the asymptotic limit is
lim
t→∞
u(t, x) =
c
0
2
, and the asymptotic profile is u(t, x) ∼
c
0
2
+c
1
exp

−
π
2
t
L
2

cos

πx
L

. (2.39)
2.4. RELAXATION TO EQUILIBRIUM SOLUTIONS 25
The constant c
0
above is (2/L)

L
0
u
0
(x)dx, thus the limit is
1
L

L
0
u
0
(x)dx, (2.40)
the average value of the initial value u
0
(x). The limits of both cases are also solutions of diffusion
equation, and in fact they are solutions which are independent of time t. Such solutions are
equilibrium solutions. For a general reaction-diffusion equation







∂u
∂t
= d∆u +f(x, u), t > 0, x ∈ Ω,
u(0, x) = u
0
(x), x ∈ Ω,
B(u) = 0, t > 0, x ∈ ∂Ω,
(2.41)
where B(u) is the boundary condition, the equilibrium solutions satisfy

d∆u +f(x, u) = 0, x ∈ Ω,
B(u) = 0, x ∈ ∂Ω.
(2.42)
A solution u
1
(x) of (2.42) is a solution of (2.41) in the sense that the function U
1
(t, x) = u
1
(x) is a
solution of (2.41) with u(0, x) = u
1
(x). It is not just a coincidence that for two examples we have
considered so far, the limit is always an equilibrium solution. We will show in next section that the
limit of a diffusion equation is usually an equilibrium solution.
Equilibrium solutions of diffusion equation are easier to find, and that will help us understand
the limit of solution when the solution itself is not so easy to find. The equilibrium equation for
(2.13)-(2.15) is
u
′′
= 0, u(0) = u(L) = 0. (2.43)
The general solution for u
′′
= 0 is u(x) = ax + b, and from the boundary conditions, we obtain
u(x) = 0 is the only equilibrium solution. The equilibrium equation for (2.36) is
u
′′
= 0, u
′
(0) = u
′
(L) = 0. (2.44)
For this problem u(x) = b for any b ∈ R is an equilibrium solution. However the limit in (2.36)
is the average value of the initial value u
0
(x) only. Why does the solution of diffusion equation
selects this particular equilibrium solution, but not the others? We will answer this question in
next section.
2.4 Relaxation to equilibrium solutions
We have obtained analytic solutions of 1-d diffusion equation in previous sections. The solutions
are in a form of infinite series. However in general, even such formulas cannot be obtained when
the spatial dimension is higher and the domain is arbitrary. (Later we will formulas of solutions
of diffusion equations for some special cases.) Here we use some integral formulas to show some
qualitative properties of solutions of general diffusion equation.
26 CHAPTER 2. DIFFUSION EQUATIONS
As a calculus warmup, we consider the diffusion equation with no-flux boundary condition







∂u
∂t
= D∆u, t > 0, x ∈ Ω,
u(0, x) = u
0
(x), x ∈ Ω,
∇u(t, x) · n(x) = 0, t > 0, x ∈ ∂Ω.
(2.45)
From the biological modelling point of view, the species only randomly moves in Ω, the growth rate
is zero, and the boundary is “sealed”. Thus the total population should not change since no one
moves in or out, and no one is born or dead. Indeed that is the basic character of the equation,
and it also explains why the limit of the solution of (2.36) is the average value of the initial value
function.
Proposition 2.5. Suppose that u(t, x) is the solution of (2.45), then

Ω
u(t, x)dx =

Ω
u
0
(x)dx. (2.46)
Proof. Integrating the left hand side of the equation in (2.45) over Ω, we obtain

Ω
∂u
∂t
(t, x)dx =
d
dt

Ω
u(t, x)dx

, (2.47)
and integrating the right hand side, we have
D

Ω
∆u(t, x)dx = D

∂Ω
∇u(t, x) · n(x)ds = 0, (2.48)
from the divergent theorem and the no-flux boundary condition. Therefore from (2.47) and (2.48),
we find that
d
dt

Ω
u(t, x)dx

= 0, (2.49)
which implies (2.46).
Next we show that the solution u(t, x) always tends to a limit constant c, and because of (2.46),
that constant must be
1
|Ω|

Ω
u
0
(x)dx. (2.50)
A rigorous proof of this fact needs many more advanced mathematical knowledge, so we only sketch
the ideas here. We multiply the equation in (2.45) by the function u(t, x),
u
∂u
∂t
= Du∆u. (2.51)
Similar to the proof of Proposition 2.5, we integrate the left hand side of (2.51):

Ω
u
∂u
∂t
dx =

Ω
∂
∂t

1
2
u
2

dx =
d
dt

Ω

1
2
u
2

dx =
1
2
d
dt

Ω
[u(t, x)]
2
dx, (2.52)
2.4. RELAXATION TO EQUILIBRIUM SOLUTIONS 27
and the right hand side:
D

Ω
u∆udx = D

∂Ω
u(∇u · n)ds −D

Ω
∇u · ∇udx. (2.53)
Here we use Green’s identity (a consequence of divergence theorem):

Ω
u∆vdx =

∂Ω
u(∇v · n)ds −

Ω
∇u · ∇vdx. (2.54)
From the boundary condition,

∂Ω
u(∇u · n)ds = 0, thus from (2.52) and (2.53), we find
1
2
d
dt

Ω
[u(t, x)]
2
dx = −D

Ω
|∇u(t, x)|
2
dx. (2.55)
Let E(t) =

Ω
[u(t, x)]
2
dx. Then (2.55) shows that E
′
(t) < 0 for all t > 0 and obviously E(t) > 0
for all t > 0. Thus lim
t→∞
E(t) ≥ 0 exists, and as t →∞, E
′
(t) →0, so
lim
t→∞

Ω
|∇u(t, x)|
2
dx = 0. (2.56)
(2.56) implies lim
t→∞
|∇u(t, x)| = 0 for any x ∈ Ω, thus the limit of u(t, x) when t → ∞ must
be a function with zero gradient, thus the limit function must be a constant u(x) = c. Since
c|Ω| =

Ω
u(x)dx = lim
t→∞

Ω
u(t, x)dx =

Ω
u
0
(x)dx, then c =

Ω
u
0
(x)dx/|Ω|. (Notice that we
implicitly assume u(t, x) has a limit u(x) when t → ∞. While such a limit indeed exists, the
proof of its existence is beyond the scope of this notes, that is why this is only a sketch of proof.)
Nevertheless, we put our findings into another proposition:
Proposition 2.6. Suppose that u(t, x) is the solution of (2.45), then u(t, x) tends to a constant
function u(x) =
1
|Ω|

Ω
u
0
(x)dx.
Similar result can also be shown for Dirichlet problem, in which the limit function is u = 0.
From this property of diffusion equation, we can get a conclusion that the diffusion equation will
make the initial value flatter and bring it eventually to an equilibrium state without spatial pattern.
This can also be indicated from the differential equation itself. Let the average of the initial value
be c. For the points where u(t, x) > c, the function is likely to be convex down, with ∆u < 0, thus
the function value will decrease over these points; and for the points where u(t, x) < c, the opposite
occurs. This is called smothering effect of diffusion. In the following we use a Maple program to
demonstrate the smothering effect and the asymptotic tendency to the equilibrium. We consider





u
t
= 4u
xx
, t > 0, x ∈ (0, 4),
u(t, 0) = u(t, 4) = 0,
u(0, x) = f(x),
(2.57)
28 CHAPTER 2. DIFFUSION EQUATIONS
where f(x) is given by the formula f(x) =

100, x ∈ [1, 3],
0, x ∈ [0, 1)
¸
(3, 4].
. From (2.32), we get the
formula of the solution of (2.57):
u(t, x) =
∞
¸
m=1
c
m
exp

−
m
2
π
2
4
t

sin

mπx
4

,
where c
m
=
1
2

4
0
f(x) sin

mπx
4

dx, m ≥ 1, t > 0, x ∈ (0, 4).
(2.58)
The coefficients c
m
in the Fourier series can easily be calculated, but we will let Maple to do that.
Here is the Maple program diffusion.mws and some of the results:
> restart;
> with(plots): n:=’n’:
> f:=x->100*Heaviside(x-1)-100*Heaviside(x-3);
> plot(f(x), x=0..4);
> c:=n->0.5*int(f(x)*sin(n*Pi*x/4),x=0..4): c(n);
200.0
cos(1/4 n Pi) - cos(3/4 n Pi)
n Pi
>uapprox1:=(x,t)->sum(c(n)*sin(n*Pi*x/4)*exp(-4*n^2*t), n=1..30);
uapprox1 :=(x,t) - >
30
¸
n=1
c(n) sin(1/4 n Pi x) exp(-4 n
2
t)
> tvals:=seq(0.07/20*i,i=0..20):
> toplot:=[seq(uapprox1(x,t),t=tvals)]:
> plot(toplot,x=0..4);
> densityplot(uapprox1(x,t),x=0..4,t=0..0.2,colorstyle=HUE);
> animate(uapprox1(x,t), x=0..4, t=0..0.1,frames=50);
This is an animation generated by approximating the solution by the partial sum with n = 30 of
the series solution. The approximation is accurate since the remainder is an exponential decaying
function with smaller exponent. Note that the initial condition f(x) can be written as the difference
of two Heaviside functions, which is defined as
H(x) =

1, x ≥ 0,
0, x < 0.
(2.59)
It is an interesting experiment to observe the changes of the spatial patterns of u(t, x) when t gets
forward. Here are a few snap shots of the animation:
2.4. RELAXATION TO EQUILIBRIUM SOLUTIONS 29
0
20
40
60
80
100
1 2 3 4
x
0
20
40
60
80
100
1 2 3 4
x
Figure 2.1: (a) f(x); (b) partial sum of Fourier series of f(x)
0
20
40
60
80
100
1 2 3 4
x
0
20
40
60
80
100
1 2 3 4
x
Figure 2.2: (a) u(0.002, x); (b) u(0.010, x)
Because of the truncation of the series solution, u(t, x) above indeed is the solution of (2.57)
with initial data
¯
f(x) =
30
¸
m=1
c
m
sin

mπx
4

dx, c
m
=
200
mπ

cos

mπ
4

−cos

mπ
4

. (2.60)
Figure 2.1 is a comparison of f and
¯
f. From the last section, we know that when t → ∞,
u(t, x) →0 exponentially fast. However the road to distinction is worth a more careful looking—it
can be dismantled into the following stages:
1. Forming a plateau From t = 0 (Figure 2.1(b)) to t = 0.002 (Figure 2.2 (a)), we can observe
that the initial data has been smoothed by the diffusion equation, and all the small wiggles
in u(0, x) are gone at t = 0.002. On the other hand, although u(0, x) is not a constant for
x ∈ (1, 3), u(0.002, x) is almost a constant on an interval slightly smaller than (1, 3). So in
the initial stage of the evolution of
¯
f under diffusion, a plateau (or flat core) of height 100
forms, with a sharp but smooth drop-off to 0 near x = 1 and x = 3.
30 CHAPTER 2. DIFFUSION EQUATIONS
0
20
40
60
80
100
1 2 3 4
x
0
20
40
60
80
100
1 2 3 4
x
Figure 2.3: (a) u(0.020, x); (b) u(0.040, x)
0
20
40
60
80
100
1 2 3 4
x
0
20
40
60
80
100
1 2 3 4
x
Figure 2.4: (a) u(0.100, x); (b) u(0.600, x)
2. Erosion of the plateau From t = 0.002 (Figure 2.2 (a)) to t = 0.01 (Figure 2.2 (b)), the flat
core keeps shrinking though the top of the core is still at the level near 100. At t = 0.01, the
erosion finally dissolves the whole flat top, and the maximum point of u(t, ·) at x = 2 looks
more like a regular local maximum point. The graph of u(0.01, x) is close to a bell-shape.
In that time interval, the inflection points of u(t, ·) are still near x = 1 and x = 3, but at
t = 0.01, the interface between u = 0 and u = 100 is no longer as sharp as when t = 0.002.
3. Bell getting round In the next stage, smoothing effect brings the concave part of the graph
down, and the convex part of the graph up. And at the same time, the convex part is shrinking
as population loss via boundary gets bigger. At t = 0.04, the graph becomes almost concave
for all x, and the bell-shape becomes an arch. (see Figure 2.3)
4. Collapse of the arch In the final stage, the arch is close to the leading term of the series:
c
1
exp(−π
2
t/4) sin(πx/4), and an exponential collapse is obvious. After t = 1, the graph can
hardly be seen. (see Figure 2.4)
2.5. ONE DIMENSIONAL CHEMICAL MIXING PROBLEM 31
2.5 One dimensional chemical mixing problem
Here we consider a simplified version of the chemical mixing problem which is discussed in Chapter
1: A tube with length L contains salt water. The cross-section (yz-direction) of the tube is so small
so we can assume that the concentration of the salt water is same for any point on a cross-section.
Let the concentration of the salt water in the tube be c(t, x) kg/m
3
, t > 0 and 0 < x < L. We
assume that the diffusion of the salt is one-dimensional in the tube, and the velocity of fluid is
ignorable, thus c(t, x) satisfies the diffusion equation:
∂c
∂t
= D
∂
2
c
∂x
2
. (2.61)
At the left hand side of the tube, a salt water solution with constant concentration C
0
kg/m
3
enters
the tube at a rate of R
0
m
3
/s, and on the right hand side of the tube, the mixed solution is removed
at the rate of R
0
m
3
/s. We also assume that the area of the cross section of the tube is A. Let A
1
and A
2
be the cross sections at x = 0 and x = L correspondingly. From the Fick’s law, the total
in-flux at x = 0 is

A
1
J · (−n)ds =

A
1
[−Dc
x
(t, 0)] · (1)ds = −DAc
x
(t, 0) = C
0
R
0
, (2.62)
and the total out-flux at x = L is

A
2
J · nds =

A
1
[−Dc
x
(t, L)] · (1)ds = −DAc
x
(t, 0) = C(t, L)R
0
, (2.63)
thus we obtain the boundary conditions
∂c(t, 0)
∂x
= −
C
0
R
0
DA
, and
∂c(t, L)
∂x
= −
c(t, L)R
0
DA
. (2.64)
To simplify the notations, we define
B =
C
0
R
0
DA
, and E =
R
0
DA
. (2.65)
Thus the initial-boundary value problem is









∂c
∂t
= D
∂
2
c
∂x
2
, t > 0, x ∈ (0, L),
∂c(t, 0)
∂x
= −B,
∂c(t, L)
∂x
+Ec(t, L) = 0,
c(t, x) = c
0
(x), t > 0, x ∈ (0, L).
(2.66)
Before attempting to solve the equation, we could find the equilibrium solutions, which satisfy
equation:





D
∂
2
c
∂x
2
= 0, x ∈ (0, L),
∂c(0)
∂x
= −B,
∂c(L)
∂x
+Ec(L) = 0.
(2.67)
32 CHAPTER 2. DIFFUSION EQUATIONS
From a simple calculation, we find the unique equilibrium solution:
c
∗
(x) =
B
E
+B(L −x) = C
0
+
C
0
R
0
DA
(L −x), (2.68)
which is a linear positive function with negative slope on [0, L], and c
∗
(L) = C
0
. To solve the
equation, we define b(t, x) = c(t, x) −c
∗
(x). Then b(t, x) satisfies









∂b
∂t
= D
∂
2
b
∂x
2
, t > 0, x ∈ (0, L),
∂b(t, 0)
∂x
= 0,
∂b(t, L)
∂x
+Ec(t, L) = 0,
b(t, x) = c
0
(x) −c
∗
(x), t > 0, x ∈ (0, L).
(2.69)
The reason we consider (2.69) instead of (2.66) is that the boundary conditions in (2.69) is homoge-
neous, thus we can use the technique of separation of variables. We assume that b(t, x) = U(t)V (x).
Then
U
′
(t)
U(t)
= D
V
′′
(x)
V (x)
= Dk, (2.70)
U(t) = e
Dkt
, and V (x) satisfies
V
′′
−kV = 0, x ∈ (0, L), V
′
(0) = 0, V
′
(L) +EV (L) = 0. (2.71)
This is a half-Neumann and half-Robin boundary condition. We use the method in Section 2.2 to
solve (2.71).
If k > 0, then V (x) = c
1
e
−
√
kx
+c
2
e
√
kx
. From the boundary conditions, we get
det

−
√
k
√
k
(E −
√
k)e
−
√
kL
(E +
√
k)e
√
kL

= 0, (2.72)
and k must satisfy
e
2
√
kL
=
√
k −E
√
k +E
. (2.73)
But the function on the left is greater than 1 when k > 0, and the one on the right is always less
than 1 when k > 0. So (2.73) has no solution, thus there are no positive eigenvalues. Similarly
k = 0 is not eigenvalue either.
If k < 0, V (x) = c
1
cos(
√
−kx) + c
2
sin(
√
−kx). From the boundary conditions, we get c
2
= 0
and c
1
[cos(
√
−kL) −
√
−k sin(
√
−kL)] = 0. Thus k satisfies
cot(
√
−kL) =
√
−k. (2.74)
From the graph of the functions f
1
(x) = cot(Lx) and f
2
(x) = x, there are exactly one intersection
point x
m
of f
1
and f
2
in ((m−1)π/L, mπ/L) for any positive integer m > 0. Then this k
m
= −x
2
m
is the m-th eigenvalue of (2.71), and the corresponding eigenfunction is cos(
√
−k
m
x). Thus the
solution of (2.69) is
b(t, x) =
∞
¸
m=1
b
m
e
−Dx
2
m
t
cos(x
m
x), where x
m
∈ (
(m−1)π
L
,
mπ
L
) satisfies cot(x
m
L) = x
m
. (2.75)
2.6. CRITICAL PATCH SIZE AND BIFURCATIONS 33
Hence the solution is the sum of the equilibrium solution and an exponential decaying Fourier
series, and the solution converges to the equilibrium solution exponentially fast. It shows that in
equilibrium state, the exiting solution also has the concentration C
0
, same as the incoming solution,
thus the amounts of salt entering and exiting the tube are balanced, but the concentration of solution
in the tube is higher than C
0
, and the concentration at x = 0 is the highest.
The example in this section can be thought as a model of a continuously polluted river. The
left end point x = 0 is the source of the pollution where the pollutant enters the river at a constant
rate, and the right end point x = L is where the river merges to a bigger river or the ocean. Our
mathematical result shows that the concentration of pollutant in a continuously polluted river is
even higher than that in the source. As a consequence, we can see that at the confluence of two
polluted rivers, the pollution can be most serious, since the concentration there is much higher than
the ones of either branches. We should be cautious that we ignore the effect of drifting here, which
could make a big difference for rivers.
2.6 Critical patch size and bifurcations
The simplest population growth model is the Malthus model dP/dt = aP for a positive constant
a > 0. In this section we consider Dirichlet boundary value problem of diffusive Malthus equation







∂P
∂t
= D
∂
2
P
∂x
2
+aP, t > 0, x ∈ (0, L),
u(t, 0) = u(t, L) = 0,
u(0, x) = u
0
(x), x ∈ (0, L).
(2.76)
We have seen that population are at a loss from the hostile boundary, but on the other hand, the
population will increase through linear reproduction. So the loss and gain are in a competition.
Let’s see what factor will determine the outcome of this battle.
Equation (2.76) can be solved in a similar way as (2.13)-(2.15). In fact, the form of the solution
only changes slightly: adding a term e
at
to reflect the exponential growth. The solution of (2.76)
is
u(t, x) =
∞
¸
m=1
c
m
exp

at −D
m
2
π
2
L
2
t

sin

mπx
L

,
where c
m
=
2
L

L
0
u
0
(x) sin

mπx
L

dx, m ≥ 1, t > 0, x ∈ (0, L).
(2.77)
However the fate of the solution may be different: since at is added to the exponent, now the
exponent may be a positive one, then u(t, x) will have a part which is exponentially increasing. In
fact, when the index m in the sum increases, the exponent gets more negative, thus whether u(t, x)
is exponentially increasing is determined by the first exponent coefficient a −(Dπ
2
/L
2
). We have
when a >
Dπ
2
L
2
or L >

D
a
π, lim
t→∞
u(t, x) = ∞,
when a <
Dπ
2
L
2
or L <

D
a
π, lim
t→∞
u(t, x) = 0.
(2.78)
34 CHAPTER 2. DIFFUSION EQUATIONS
If we assume that the growth rate a and D are intrinsic parameters determined by the nature of
species and environment, then the size of the habitat L will play the deciding role here. The number
L
0
=

D/aπ is called the critical patch size, as the population cannot survive if the habitat is too
small (the outgoing flux wins over the growth), but the population will not only survive but thrive
to an exponential growth if the habitat is large enough (growth outpaces the emigration).
We say that a bifurcation occurs at L = L
0
as the equilibrium solution u = 0 changes stability
type from stable for L < L
0
to unstable for L > L
0
. We notice that at L = L
0
, the equation has
other equilibrium solutions u(x) = sin(πx/L), the eigenfunction associated with first eigenvalue
k
1
= −π
2
/L
2
. We will have more discussion on bifurcation problems when the growth rate is
nonlinear. From (2.78), we can also use a or D as bifurcation parameter instead of L, sometime
that is more convenient since we do not need to change the domain, but we get equivalent results.
The concept of critical patch size is related to habitat fragmentation in ecological studies.
Habitat fragmentation
1
is the breaking up of a continuous habitat, ecosystem, or land-use type
into smaller fragments, which is considered to be one of several spatial processes in land transfor-
mation. It is commonly used in relation to the fragmentation of forests. Habitat fragmentation is
mainly caused by human activities such as logging, conversion of forests into agricultural areas and
suburbanization, but can also be caused by natural processes such as fire. From our mathematical
result in this section, if the the fragmentation of habitat limits the spatial movement of the plants
and animals, then the species can become extinct.
2.7 Separation of variables: rectangles
Most of this chapter, we have exclusively considered one dimensional spatial domain, for the sim-
plicity of the mathematical analysis. However any interesting application happens in 2-D or 3-D
spaces. In this section, we discuss the separation of variables in a two-dimensional rectangle.
We first consider a linear diffusion reaction equation:









∂u
∂t
= D

∂
2
u
∂x
2
+
∂
2
u
∂y
2

+λu, t > 0, (x, y) ∈ R = (0, a) ×(0, b),
∇u · n = 0, (x, y) ∈ ∂R,
u(0, x, y) = u
0
(x, y), (x, y) ∈ R,
(2.79)
where D, λ, a, b > 0, and ∂R is the boundary of rectangle R. From the method of separation of
variables, we assume that
u(t, x, y) = U(t)V (x, y). (2.80)
Then similar to one-dimensional case, we obtain
U
′
(t) = (Dk +λ)U(t), (2.81)
1
From TEMS, Terrestrial Ecosystem Monitoring Sites, http://www.fao.org/gtos/tems/
2.7. SEPARATION OF VARIABLES: RECTANGLES 35
and



∂
2
V
∂x
2
+
∂
2
V
∂y
2
= kV, t > 0, (x, y) ∈ R = (0, a) ×(0, b),
∇V · n = 0, (x, y) ∈ ∂R.
(2.82)
For (2.82), we need a further separation of variables:
V (x, y) = W(x)Z(y), (2.83)
Then we obtain that
W
′′
(x)
W(x)
+
Z
′′
(y)
Z(y)
= k, (2.84)
so both W
′′
(x)/W(x) and Z
′′
(y)/Z(y) must be constants. On the other hand, the Neumann
boundary condition implies that W
′
(0) = W
′
(a) = 0 and Z
′
(0) = Z
′
(b) = 0. Therefore W and Z
satisfy
W
′′
(x) = k
1
W(x), x ∈ (0, a), W
′
(0) = W
′
(a) = 0, (2.85)
Z
′′
(y) = k
2
Z(y), y ∈ (0, b), Z
′
(0) = Z
′
(b) = 0, (2.86)
and
k = k
1
+k
2
. (2.87)
Both (2.85) and (2.86) are familiar eigenvalue problems in one dimensional space which have been
studied in Chapter 2. The eigenvalues and eigenfunctions of (2.85) are
k
1n
= −
n
2
π
2
a
2
, W
n
(x) = cos

nπx
a

, n = 0, 1, 2, · · · , (2.88)
and the eigenvalues and eigenfunctions of (2.86) are
k
2m
= −
m
2
π
2
b
2
, Z
m
(y) = cos

mπy
b

, m = 0, 1, 2, · · · . (2.89)
Therefore the eigenvalues and eigenfunctions of (2.82) are
k
n,m
= −
n
2
π
2
a
2
−
m
2
π
2
b
2
, V
n,m
(x, y) = cos

nπx
a

· cos

mπy
b

, (2.90)
where n, m = 0, 1, 2, · · · . Thus we find that the solution of (2.79) is
u(t, x, y) =
∞
¸
n=0,m=0
c
n,m
e
(Dk
n,m
+λ)t
cos

nπx
a

· cos

mπy
b

, (2.91)
where c
n,m
can be determined by the initial conditions.
The spatial patterns of the eigenfunctions for rectangle are much more complex than those
of one-dimensional. Here we consider a square R = (0, π) × (0, π), then the first five distinctive
eigenvalues are (see the contour graph of eigenfunctions below)
k
0,0
= 0, k
0,1
= k
1,0
= −1, k
1,1
= −2, k
2,0
= k
0,2
= −4, k
2,1
= k
1,2
= −5, (2.92)
36 CHAPTER 2. DIFFUSION EQUATIONS
0
0.5 1
1.5 2
2.5 3
y
0
0.5
1
1.5
2
2.5
3
x
0
0.5 1
1.5 2
2.5 3
y
0
0.5
1
1.5
2
2.5
3
x
0
0.5 1
1.5 2
2.5 3
y
0
0.5
1
1.5
2
2.5
3
x
0
0.5 1
1.5 2
2.5 3
y
0
0.5
1
1.5
2
2.5
3
x
Figure 2.5: (From left to right) (a) V
0,0
= 1; (b) V
0,1
= cos y; (c) V
1,0
= cos x; (d) V
1,1
= cos x· cos y.
0
0.5 1
1.5 2
2.5 3
y
0
0.5
1
1.5
2
2.5
3
x
0
0.5 1
1.5 2
2.5 3
y
0
0.5
1
1.5
2
2.5
3
x
0
0.5 1
1.5 2
2.5 3
y
0
0.5
1
1.5
2
2.5
3
x
0
0.5 1
1.5 2
2.5 3
y
0
0.5
1
1.5
2
2.5
3
x
Figure 2.6: (From left to right) (a) V
0,2
= cos(2y); (b) V
2,0
= cos(2x); (c) V
1,2
= cos x · cos(2y); (d)
V
2,1
= cos(2x) · cos y.
We notice that many eigenvalues k
m,n
have multiplicity more than 1, for example k
m,n
=
k
n,m
= −(m
2
+n
2
). In such case, the eigenspace is of higher dimensional, so more possible spatial
patterns are generated. For example, k
1,0
= k
0,1
= −1, and the eigenspace is span{cos x, cos y},
thus k
1
cos x +k
2
cos y is a spatial pattern for any k
1
and k
2
. In particular, the patterns generated
by cos x ±cos y are symmetric with respect to one of the diagonal lines.
0
0.5 1
1.5 2
2.5 3
y
0
0.5
1
1.5
2
2.5
3
x
0
0.5 1
1.5 2
2.5 3
y
0
0.5
1
1.5
2
2.5
3
x
Figure 2.7: (From left to right) (a) cos x + cos y; (b) cos y −cos x.
The solution (2.91) is always an exponential growth if λ > 0. Indeed, from the equation, the
no-flux boundary condition prevents the emigration through the boundary, and λ > 0 implies a
positive growth rate, thus the total population in R has an exponential growth. We can consider
corresponding Dirichlet boundary problem:









∂u
∂t
= D

∂
2
u
∂x
2
+
∂
2
u
∂y
2

+λu, t > 0, (x, y) ∈ R = (0, a) ×(0, b),
u(t, x, y) = 0, (x, y) ∈ ∂R,
u(0, x, y) = u
0
(x, y), (x, y) ∈ R.
(2.93)
2.8. STEADY STATES IN RADIALLY SYMMETRIC DOMAINS AND TRANSIT TIMES 37
Similar to the method above, the solution is
u(t, x, y) =
∞
¸
n=1,m=1
b
n,m
e
(Dk
n,m
+λ)t
sin

nπx
a

· sin

mπy
b

, (2.94)
where b
n,m
can be determined by the initial conditions, and k
n,m
is defined in (2.90). Similar to
last section, the persistence/extinction of the population can be determined:
when
λ
Dπ
2
>
1
a
2
+
1
b
2
, lim
t→∞
u(t, x) = ∞,
when
λ
Dπ
2
<
1
a
2
+
1
b
2
, lim
t→∞
u(t, x) = 0.
(2.95)
2.8 Steady states in radially symmetric domains and transit times
Chapter 2 Exercises
1. Find the eigenvalues and eigenfunctions of homogeneous Neumann boundary problem
y
′′
= ky, x ∈ (0, L), y
′
(0) = y
′
(L) = 0. (2.96)
2. Find the eigenvalues and eigenfunctions of periodic boundary problem
y
′′
= ky, x ∈ (0, π), y(0) = y(π), y
′
(0) = y
′
(π). (2.97)
3. Find the eigenvalues and eigenfunctions of Robin boundary problem (b > 0)
y
′′
= ky, x ∈ (0, L), y
′
(0) = by(0), y
′
(L) = −by(L). (2.98)
4. Find the eigenvalues and eigenfunctions of the problem
y
′′
+ 2y
′
+ky = 0, x ∈ (0, π), y(0) = y(π) = 0. (2.99)
5. In Example 2.1, use Maple to find the numerical values of the first three eigenvalues, and plot
the graphs of the corresponding eigenfunctions.
6. Show that the solution of





u
t
= Du
xx
, t > 0, x ∈ (0, L),
u
x
(t, 0) = u
x
(t, L) = 0,
u(0, x) = u
0
(x), x ∈ (0, L),
(2.100)
is
u(t, x) =
c
0
2
+
∞
¸
m=1
c
m
exp

−D
m
2
π
2
L
2
t

cos

mπx
L

,
where c
m
=
2
L

L
0
u
0
(x) cos

mπx
L

dx, m ≥ 0, t > 0, x ∈ (0, L).
(2.101)
38 CHAPTER 2. DIFFUSION EQUATIONS
7. (a) Find the solution of





u
t
= 4u
xx
, t > 0, x ∈ (0, 1),
u
x
(t, 0) = u
x
(t, 1) = 0,
u(0, x) = f(x),
(2.102)
where f(x) is given by the formula f(x) =

1, x ∈ [0, 0.5],
−1, x ∈ (0.5, 1].
(b) Modify the Maple program in Section 2.4 to simulate the solution of (2.102).
8. Find the solutions of equilibrium equation
u
′′
= 0, x ∈ (0, 1), u(0) + 3u
′
(0) = 5, u(1) −5u
′
(1) = 7. (2.103)
9. Find the solutions of equilibrium equation
u
′′
+ 3u
′
= 0, x ∈ (0, 1), u(0) = 5, u
′
(1) = 4. (2.104)
10. Suppose that u(t, x) is the solution of







∂u
∂t
= D∆u, t > 0, x ∈ Ω,
u(0, x) = u
0
(x), x ∈ Ω,
∇u(t, x) · n(x) +a · u(t, x) = 0, t > 0, x ∈ ∂Ω.
(2.105)
where a > 0. Show that the mean-squared value I(t) =

Ω
[u(t, x)]
2
dx is strictly decreasing
unless u(t, x) ≡ 0.
11. Suppose that u(t, x) is the solution of the convective-diffusion equation: (D > 0 and V > 0)









∂u
∂t
= D
∂
2
u
∂x
2
−V
∂u
∂x
, x ∈ (0, π),
u(0, x) = u
0
(x), x ∈ (0, π),
∂u
∂x
(t, 0) =
∂u
∂x
(t, π) = 0.
(2.106)
Show that the total population

1
0
u(t, x)dx is a constant if u
0
(x) = sin x but not a constant
if u
0
(x) = cos(x).
12. Suppose that in the chemical mixing problem in Section 2.5, the salt solution drifts to the
right with velocity V . Then c(t, x) satisfies the convective-diffusion equation:
∂c
∂t
= D
∂
2
c
∂x
2
−V
∂c
∂x
, (2.107)
with D, V > 0. Suppose that c(t, x) satisfies the same initial and boundary conditions in
(2.66).
2.8. STEADY STATES IN RADIALLY SYMMETRIC DOMAINS AND TRANSIT TIMES 39
(a) Find the equilibrium solution of (2.107).
(b) Show that the equilibrium solution is positive, decreasing and concave down, and explain
the effect of the drifting to the geometric properties of the equilibrium solution.
(c) Compare the exiting flux −Dc
x
(t, L) of the equilibrium solution with that of (2.67).
13. Suppose that in the chemical mixing problem in Section 2.5, the salt is consumed in a chemical
reaction at a rate k. Then c(t, x) satisfies a linear diffusion equation:
∂c
∂t
= D
∂
2
c
∂x
2
−kc, (2.108)
with D, k > 0. Suppose that c(t, x) satisfies the same initial and boundary conditions in
(2.66).
(a) Find the equilibrium solution of (2.108). (You can use Maple to find the general formula.)
(b) Show that the equilibrium solution is positive, decreasing and concave up, and explain
the effect of the reaction to the geometric properties of the equilibrium solution.
(c) Compare the exiting flux −Dc
x
(t, L) of the equilibrium solution with that of (2.67).
(d) Plot the graphs of the equilibrium solutions of (2.67), (2.107) and (2.108) with same
parameters in Maple.
14. Determine the critical patch for a diffusive Malthus equation with Robin boundary conditions:







∂P
∂t
= D
∂
2
P
∂x
2
+aP, t > 0, x ∈ (0, L),
u
x
(t, 0) = bu(t, 0), u
x
(t, L) = −bu(t, L),
u(0, x) = u
0
(x), x ∈ (0, L).
(2.109)
15. Consider the doubly periodic boundary value problem:













∂u
∂t
= D

∂
2
u
∂x
2
+
∂
2
u
∂y
2

+λu, t > 0, (x, y) ∈ R = (0, a) ×(0, b),
u(0, y) = u(a, y), u
x
(0, y) = u
x
(a, y),
u(x, 0) = u(x, b), u
y
(x, 0) = u
y
(x, b),
u(0, x, y) = u
0
(x, y), (x, y) ∈ R.
(2.110)
Find the eigenvalues and eigenfunctions of the problem, and the series representation of the
solution.
16. Consider a convective-diffusion equation







∂u
∂t
=
∂
2
u
∂x
2
+ 4
∂u
∂x
+ 8u, t > 0, x ∈ (0, L),
u(t, 0) = 0, u(t, L) = 0,
u(0, x) = u
0
(x), x ∈ (0, L).
(2.111)
Find the solution of the equation with the following steps:
40 CHAPTER 2. DIFFUSION EQUATIONS
(a) Use separation of variables method to show that if u(t, x) = U(t)V (x) is a solution, then
for some constant k, U and V satisfy
U
′
(t) = kU(t), V
′′
(x) + 4V
′
(x) + 8V (x) = kV (x), V (0) = V (L) = 0.
(b) Find the eigenvalues and eigenfunctions of
V
′′
(x) + 4V
′
(x) + 8V (x) = kV (x), V (0) = V (L) = 0.
(Hint: treat the cases of k > 4, k = 4 and k < 4 separately.)
(c) Find the solution of the equation in a series form.
(Hint: c
n
=
2
L

L
0
e
2x
u
0
(x) sin

nπx
L

dx.)
(d) Determine the critical patch L
0
of the problem.
(e) Describe the population distribution qualitatively when L > L
0
.
(f) Suppose that the species lives in a river with length L, and the convection is due to the
drifting of the river. Explain your mathematical results in this context.