NISTIR 8027 

Stable Explicit Time Marching in 

Well-Posed or Ill-Posed Nonlinear 

Parabolic Equations 

Alfred S. Carasso 

Applied and Computational Mathematics Division 

Information Technology Laboratory 

This publication is available free of charge from: 

http://dx.doi.org/10.6028/NIST.IR.8027
NISTIR 8027 

Stable Explicit Time Marching in 

Well-Posed or Ill-Posed Nonlinear 

Parabolic Equations 

Alfred S. Carasso 

Applied and Computational Mathematics Division 

Information Technology Laboratory 

This publication is available free of charge from: 

http://dx.doi.org/10.6028/NIST.IR.8027 

October 2014 

U.S. Department of Commerce 

Penny Pritzker, Secretary 

National Institute of Standards and Technology 

Willie May, Acting Under Secretary of Commerce for Standards and Technology and Acting Director 
STABLE EXPLICIT TIME MARCHING IN WELL-POSED OR
 

ILL-POSED NONLINEAR PARABOLIC EQUATIONS
 

ALFRED S. CARASSO∗

Abstract. This paper analyzes an effective technique for stabilizing pure explicit time differencing in the numerical computation of multidimensional nonlinear parabolic equations. The

method uses easily synthesized linear smoothing operators at each time step to quench the instability. Smoothing operators based on positive real powers of the negative Laplacian are helpful, and

(−Δ)p can be realized efficiently in rectangular domains using FFT algorithms. The stabilized explicit scheme requires no Courant restriction on the time step Δt, and is of great value in computing

well-posed parabolic equations on fine meshes, by simply lagging the nonlinearity at the previous

time step. Such stabilization leads to a distortion away from the true solution. However, that error

is often small enough to allow useful results in many problems of interest.

The stabilized explicit scheme is also stable when run backward in time. This allows for relatively easy and useful computation of a significant class of multidimensional nonlinear backward

parabolic equations, and complements the quasi-reversibility method. In the canonical case of linear

autonomous selfadjoint backward parabolic equations, with solutions satisfying prescribed bounds,

it is proved that the stabilized explicit scheme can produce results that are nearly best-possible. Such

backward reconstructions are of increasing interest in environmental forensics, where contaminant

transport is often modeled by advection dispersion equations.

The paper uses fictitious mathematically blurred 512 × 512 pixel images as illustrative examples.

Such images are associated with highly irregular jagged intensity data surfaces that can severely challenge ill-posed nonlinear reconstruction procedures. Instructive computational experiments demonstrate the capabilities of the method in 2D rectangular regions.

Key words. FFT Laplacian stabilization, forward or backward nonlinear parabolic equations;

non-integer power Laplacian; nonlinear image deblurring; quasi-reversibility method; stabilized explicit scheme; Van Cittert iteration.

AMS subject classifications. 35K55, 35R25, 65M06, 65M30, 68U10.

1. Introduction. This paper develops an effective technique for stabilizing pure

explicit time differencing in the numerical computation of multidimensional nonlinear

parabolic equations. As is well-known [21, p. 208], while explicit schemes are highly

desirable because of their simplicity, such schemes are seldom used as they entail prohibitive Courant stability restrictions on the time step Δt. Rather, unconditionally

stable implicit schemes are commonly used, even though such schemes require sophisticated and computationally intensive solutions of the resulting algebraic systems of

difference equations at each time step. The present stabilized explicit scheme requires

no Courant restriction on Δt, and hence would be of great value in computing multidimensional, well-posed, nonlinear parabolic equations on fine meshes, by simply

lagging the nonlinearity at the previous time step.

For ill-posed backward parabolic equations, classical stepwise time marching algorithms, whether explicit or implicit, are necessarily unconditionally unstable, [21, p.

59]. However, remarkably, the stabilized explicit scheme is stable even when run backward in time. As will be shown in the examples below, this scheme opens the door

to relatively easy and useful computation of multidimensional nonlinear backward

parabolic equations, a category of problems that is widely believed to be intractable.

The method complements the quasi-reversibility method [13], whose formulation is

primarily oriented toward linear problems. There is considerable interest in such

backward reconstructions, notably in the developing field of environmental forensics

∗Applied and Computational Mathematics Division, National Institute of Standards and Technology, Gaithersburg, MD 20899. (alfred.carasso@nist.gov).

1
2 A. S. CARASSO

[3], [4], [17], [23], where contaminant transport is often modeled by parabolic advection dispersion equations [14], [22]. However, the class of nonlinear problems that can

be handled by the present method is limited, although it is significant.

Stabilization of the explicit scheme is accomplished by applying an easily synthesized linear smoothing operator at each time step to quench the instability. Smoothing

operators based on positive non-integer powers of the negative Laplacian have been

found helpful. Inevitably, such intervention leads to a distortion away from the true

solution. This is the stabilization penalty. The scheme is useful on a given problem

only to the extent that the accompanying distortion is sufficiently small. When this

is the case, the method is an effective preliminary exploration tool that can produce

valuable information in quite difficult problems. Such information can subsequently

be confirmed using more elaborate computational methods. As noted below, the

method assumes the validity of Eq. (2.4). That hypothesis appears to be verified in

numerous computational experiments.

While the computation of nonlinear problems is the principal aim of this work,

as well as the focus of the illustrative examples, valuable insight is provided by the

analysis in Sections 2 through 5, which addresses a companion class of linear problems.

Consideration of the transparent linear case enables isolating the stabilization penalty,

evaluating its significance, and interpreting it as a form of regularization. In the

canonical case of linear autonomous selfadjoint backward parabolic equations with

solutions satisfying prescribed bounds [10], [15], a sharp estimate for the uncertainty in

reconstruction can be obtained using logarithmic convexity arguments. The analysis

in Eqs. (3.18–3.19) shows that the error in the stabilized explicit scheme differs from

that fundamental uncertainty only by the stabilization penalty. Moreover, it is shown

that that penalty can be expected to be small in a wide class of parabolic problems.

The present paper complements and extends earlier work discussed in [7], by

providing the analytical justification not included in [7], by considering the wellposed forward problem as well as the backward problem, and by exhibiting successful

computation of nonlinear problems for which the competing Van Cittert method,

previously developed and applied in [7], fails.

1.1. Use of fictitious images as illustrative test examples. Images can provide useful and challenging computational test examples for 2D nonlinear parabolic

equations. Using a sharp image as initial data in a well-posed forward parabolic

equation results in a fictitious blurred image, yet one in which the effects of various nonlinearities can be easily gauged with the naked eye. Using a mathematically

blurred image as input data in the ill-posed backward parabolic problem, poses a

severe test for an inverse reconstruction algorithm. As shown in Figure 1.1, highly

complicated intensity data are typically associated with images of familiar objects.

Such complex yet meaningful data are not easily generated analytically by combining

mathematical formulae. Also, given a familiar object, the degree of success in nonlinear backward reconstruction is often immediately apparent. On the other hand,

comparing the reconstructed intensity surface with the original intensity surface may

not always be helpful; a slight blurring of the image, one that does not affect object

recognition, can cause nonintuitive changes in the original surface.

Very little is known either analytically or computationally about nonlinear backward parabolic equations. Examples of physically plausible yet false reconstructions

have been found [6]. Given a specific nonlinear equation, preliminary exploration

of the associated fictitious image deblurring problem provides valuable information

about that equation, even when the intended scientific or engineering application is
STABLE EXPLICIT TIME-MARCHING IN NONLINEAR PARABOLIC EQUATIONS 3

Images are associated with complicated

 surface plots of intensity values

Shirley Temple Plot of intensity values

Pyramids at Giza Plot of intensity values

Fig. 1.1. Reconstruction of images requires recovery of highly irregular, jagged intensity data

surfaces that severely challenge ill-posed computational algorithms. Accordingly, fictitious nonlinearly blurred images provide instructive test examples in this paper.
4 A. S. CARASSO

unrelated to imaging. Accordingly, images will be used as test examples to illustrate

the merits of the stabilized explicit scheme.

2. A companion class of autonomous linear second order parabolic

equations. Useful insight into the nonlinear computations can be gained by considering the following class of problems. Let L be a linear autonomous second order

elliptic differential operator with smooth variable coefficients in a region Ω in Rn, with

homogeneous Dirichlet or Neumann boundary conditions on ∂Ω. Let 2 denote the

norm on L2(Ω). We assume that L has a complete set of orthonormal eigenfunctions

{φm}∞ in L2 m=1 (Ω), corresponding to eigenvalues {λm}m

∞=1 satisfying

(2.1) 0 ≤ |λ1| ≤ |λ2| ≤ · · · ≤ |λm| · · · ↑ ∞.

The eigenvalues λm of L need not all be of the same sign, and the evolution equation

(2.2) ut = −Lu, t > 0, u(0) = f,

may not be well-posed. However, we assume that if a solution exists in Eq. (2.2), it

is unique.

We shall also consider the auxiliary well-posed parabolic equation on Ω

(2.3) vt = Δv, t > 0, v(0) = f,

with the same homogeneous boundary conditions on ∂Ω as in Eq. (2.2). The Laplacian operator is selfadjoint and isotropic, whereas this may not be the case with the

differential operator L. Let L∗ denote the adjoint of L. With real p, q > 1, we

make the following assumption regarding the semigroups generated by −(−Δ)p and

−(L∗L)q/2 : Given any ω > 0, and q with 1 < q ≤ 3, there exist constants Γ, E > 0,

and p ≥ q, such that for all g ∈ L2(Ω) and sufficently small time step Δt,

(2.4) exp{−EΔt(−Δ)p}g 2≤ Γ exp{−ωΔt(L∗L)

q/2}g 2 .

 The linear operator H = (exp{−EΔt(−Δ)p}) exp{ωΔt(L∗L)q/2} is well-defined on

 the range of exp{−ωΔt(L∗L)q/2} , which is dense in L2(Ω). The inequality in Eq.

(2.4) would follow if it can be shown that H can be extended to a bounded operator

on all of L2(Ω), with H 2≤ Γ.

Eq. (2.4) appears to be validated in numerous computational experiments. Results of a similar nature are known in the theory of Gaussian estimates for heat kernels.

See e.g. [1], [2], [18], [19], and the references therein.

3. Stabilizing the pure explicit scheme. The discussion here will focus exclusively on time discretizations, and the space variables will remain continuous in

all the equations below. For the equation ut = −Lu, 0 ≤ t ≤ T, in Eq. (2.1), the n+1 pure explicit scheme, u = un − ΔtLun, n = 0, 1, 2, · · · , NT − 1, Δt = T/NT , is

a first order accurate scheme, requiring a stringent stability restriction on the time

step Δt, when Eq. (2.1) represents a forward, well-posed, multidimensional parabolic

equation [21]. For this reason, the pure explicit scheme is seldom used. If u(x,t) n+1 n is the exact unique solution of ut = −Lu, then u = u − ΔtLun + τ n, where

τ n = 1/2(Δt)2(L2u)(t

 ̃), with the intermediate time value t ̃ lying between nΔt and

(n + 1)Δt. The quantity τ n is the truncation error.

With fixed ω, Δt, and q, 1 < q ≤ 3, let S and A denote the linear operators in

L2(Ω) defined by

(3.1) S = exp{−ωΔt(L∗L)

q/2}, A = S(I − ΔtL).
STABLE EXPLICIT TIME-MARCHING IN NONLINEAR PARABOLIC EQUATIONS 5

We shall consider a compensated pure explicit scheme for ut = −Lu, t > 0, u(0) = f,

defined as follows

w

0 = f,

w

n+1 = Awn ≡ S{w

n − ∆tLwn}, n = 0, 1, 2, · · · , NT − 1, ∆t = T /NT .

(3.2)

Lemma 1. Let {λm}∞

m=1 be the eigenvalues of L satisfying Eq. (2.1). With real

q > 1, fix a positive integer Q such that |λQ| > 0, and choose ω ≥ |λQ|

1−q

. Then,

(1 + ∆t|λm|) exp{−ω∆t|λm|

q

(3.3) } ≤ 1 + ∆t|λQ|, m ≥ 1.

Proof : The inequality in Eq. (3.3) is valid for 1 ≤ m ≤ Q in view of Eq. (2.1). For

m > Q,

exp{−ω∆t|λm|

q

} ≤ exp{−ω∆t|λm||λQ|

q−1

(3.4) } ≤ exp{−∆t|λm|},

since ω|λQ|

q−1 ≥ 1. Also, exp{−∆t|λm|} ≤ (1 + ∆t|λm|)

−1

, since 1 + x ≤ e

x

for real

x. Hence, for m > Q, (1 + ∆t|λm|) exp{−ω∆t|λm|

q} ≤ 1. QED.

Lemma 2. Let q, Q, ω, be as in Lemma 1. Then, k A k2≤ 1 + ∆t|λQ|, the

compensated pure explicit scheme Eq. (3.2) is unconditionally stable, and w

n satisfies

k w

n

(3.5) k2≤ exp{|λQ|n∆t} k f k2, n = 0, 1, 2, · · · , NT .

Proof : Expanding in the orthonormal eigenfunctions {φm(x, y)}∞

m=1 of L, let

w

n

m =< wn, φm >. Then,

w

n+1 =

X∞

m=1

exp{−ω∆t|λm|

q

} (1 − ∆tλm) w

n

(3.6) mφm.

Hence, from Parseval’s relation and Lemma 1,

k w

n+1 k2 ≤ sup

m

{(1 + ∆t|λm|) exp[−ω∆t|λm|

q

]} k w

n

k2,

≤ (1 + ∆t|λQ|)

n+1 (3.7) k f k2≤ exp{|λQ|(n + 1)∆t} k f k2 .

QED.

Lemma 3. Let |||h|||2,∞ denote

|||h|||2,∞ = sup

0≤t≤T

(3.8) {k h(., t) k2} ,

and let u(t) be the unique solution of ut = −Lu, u(0) = f, on 0 ≤ t ≤ T. Then, with

S as in Eq.(3.1),

k τ

n

k2 ≤ 1/2(∆t)

2

|||L

2u|||2,∞,

k (I − S)u

n

k2 ≤ ω∆t |||(L

∗L)q/2u|||2,∞,

∆t k (I − S)Lunk2 ≤ ω(∆t)

2

|||(L

∗L)q/2

(3.9) Lu|||2,∞.
6 A. S. CARASSO

Proof : The inequality for the truncation error τ

n in Eq. (3.9) follows from Taylor’s

formula. Expanding in the orthonormal eigenfunctions {φm(x, y)}∞

m=1 of L, and using

1 − e

x ≤ −x for real x, we have, with un

m =< un, φm >

k (I − S)u

n

k

2

2 =

X∞

m=0



1 − e

−ω∆t|λm|

q

2

|u

n

m|

2 ≤

X∞

m=0

(ω∆t|λm|

q

)

2

|u

n

m|

2

= (ω∆t)

2

k (L

∗L)q/2un

k

2

2≤ (ω∆t)

2

|||(L

∗L)q/2u|||2

(3.10) 2,∞.

This proves the second inequality in Eq. (3.9), as well as the last inequality which is

a corollary of the second. QED.

Theorem 1. Let u

n be the unique solution of ut = −Lu, u(0) = f, at t = n∆t.

Let w

n be the corresponding solution of the S-compensated scheme in Eq. (3.2), and

let ES(t) ≡ u

n − wn, denote the error at t = n∆t, n = 0, 1, 2, · · · , NT . Then, with

q, Q, ω, as in Lemma 1, and |||h|||2,∞ as in Eq. (3.8),

k ES(t) k2 ≤ e

|λQ|t

k ES(0) k2 +

n

(e

|λQ|t − 1)/|λQ|

o

×

n

ω |||(L

∗L)q/2u|||2,∞ + ω∆t |||(L∗L)q/2Lu|||2,∞ + (∆t/2) |||L2u|||2,∞

o

(3.11) .

Proof : Let g

n = τn + (I − S)un + ∆t(I − S)Lun. Then, un+1 = Aun + gn, while

w

n+1 = Awn. Therefore

ESn+1 = u

n+1 − wn+1 = An

e

n + gn = An+1ES0 + ∆t

Xn

j=0

A

n−j

g

j

(3.12) /(∆t).

Hence, using Lemma 2 with t = (n + 1)∆t,

k ES(t) k2 ≤ e

|λQ|t

k ES(0) k2 + {|||g|||2,∞/∆t} ∆t

Xn

j=0

k A

n−j

k2,

≤ e

|λQ|t

k ES(0) k2 + {|||g|||2,∞/∆t}

Z t

0

e

|λQ|(t−u)

du

= e

|λQ|t

k ES(0) k2 + {|||g|||2,∞/∆t} (e

|λQ|t

(3.13) − 1)/|λQ|.

Next, using Lemma 3 to estimate {|||g|||2,∞/∆t}, we obtain Eq. (3.11) from Eq.

(3.13). QED.

3.1. The stabilization penalty in the well-posed forward problem. Consider the case of the well-posed forward parabolic problem ut = −Lu, t > 0, u(0) = f,

with L a non-negative autonomous selfadjoint elliptic operator, and the given initial

data f(x) known exactly, or with high accuracy. Then, ES(0) = 0 in Eq. (3.11).

Choosing ω = |λQ|

1−q

, there remains a non-vanishing residual in Eq. (3.11), the

stabilization penalty, as ∆t ↓ 0,

k ES(t) k2≤ (e

|λQ|T − 1) |||(L∗L)q/2u|||2,∞/|λQ|q

(3.14) + O(∆t), 0 ≤ t ≤ T.

Clearly, if T is large, the accumulated distortion may become unacceptably large as

t ↑ T, and the stabilized explicit scheme is not useful in that case. On the other
STABLE EXPLICIT TIME-MARCHING IN NONLINEAR PARABOLIC EQUATIONS 7

hand, if T is small, as is the case in many multidimensional problems of interest

involving small values of t, it may be possible to choose q > 2 and large |λQ|, with

small enough |λQ|T such that (e

|λQ|T − 1)/|λQ|q

is quite small. In that case, the

stabilization penalty remains acceptably small on 0 ≤ t ≤ T. We give several such

examples below.

3.2. The stabilization penalty in the ill-posed backward problem. In the

ill-posed backward parabolic case where L is a non-positive self-adjoint operator, the

class of initial values for which there exists a solution is very restricted. Here, the given

initial data f is seldom a highly accurate approximation to that unknown restricted

data, and the input data error contribution in Eq. (3.11) cannot be neglected. With

ω = |λQ|

1−q

, this leads to

k ES(t) k2 ≤ e

|λQ|t

k ES(0) k2

+ {(e

|λQ|t − 1)/|λQ|q

} |||(L

∗L)q/2

(3.15) u|||2,∞ + O(∆t), 0 ≤ t ≤ T.

Again, with large T, the non-vanishing residual in Eq. (3.15) will lead to large errors

as t ↑ T, and the backward scheme is not useful in such cases. For problems with

small T, Eq. (3.15) is more interesting. Let the given initial data f(x) approximate

the unknown true data u(x, 0) to within δ > 0 in the L

2 norm, and let the exact

solution u(., t) of the backward problem satisfy a known prescribed L

2 bound M at

t = T. Thus,

(3.16) k ES(0) k2≤ δ, k u(., T) k2≤ M.

A judicious choice for |λQ| is then

(3.17) |λQ| = (1/T) log(M/δ)

This leads to

k ES(t) k2 ≤ Mt/T δ

(T −t)/T

+ {(e

|λQ|t − 1)/|λQ|q

} |||(L

∗L)q/2

(3.18) u|||2,∞ + O(∆t), 0 ≤ t ≤ T.

The first term on the right in Eq. (3.18) is best-possible, and represents the fundamental uncertainty in solutions of autonomous selfadjoint backward parabolic equations,

satisfying a prescribed L

2 bound M, given initial data with L2

error at most δ, [10],

[15]. The stabilization penalty augments this basic logarithmic convexity error bound

by the addition of the second term in Eq. (3.18). In many problems of interest

this added term can be quite small. For example, with parameter values such as

T = 10−3, M = 102, δ = 10−4, q = 2.5, we have M/δ = 106, log(M/δ) = 13.816.

From Eq. (3.17), e

|λQ|T = 106

, |λQ| = 13816. Hence, |λQ|

q = 2.244 × 1010, ω =

|λQ|

1−q = 6.16 × 10−7

. We then obtain from Eq. (3.18),

k ES(t) k2 < Mt/T δ

(T −t)/T

+ (4.46 × 10−5) |||(L

∗L)q/2

(3.19) u|||2,∞ + O(∆t), 0 ≤ t ≤ T.

Remark 1. The above analysis is valid in general domains Ω ∈ Rn, and is helpful

in the computation of nonlinear problems in the following way. If L is a useful linear

representative for the nonlinear spatial differential operator Le, one may synthesize the
8 A. S. CARASSO

compensating operator S in Eq. (3.1) by precomputing a sufficiently large number K

of eigenpairs {λm, φm} of the linear elliptic operator L, spanning a subspace VK in

L

2

(Ω). See e.g. [5], [8], [20]. Such a time-consuming task may pay large dividends if

the resulting stabilized explicit scheme can actually be successfully applied to multidimensional nonlinear parabolic equations, equations with solutions well-approximated

by elements of VK. Expanding in the K precomputed eigenfunctions of L, we may

approximate the compensating operator S in Eq. (3.1) as follows for v ∈ VK ⊂ L

2

(Ω),

Sv =

X

K

m=1

exp{−ω∆t|λm|

q

}

Z

Ω

v(x)φm(x)dx(3.20) φm(x).

In a practical application, the parameters ω > 0 and q > 1 in Eq. (3.20) must first

be selected and adjusted interactively, based on prior knowledge about the solution.

Such adjustment is similar to the manual tuning of an FM station, or the focusing of

binoculars. A priori knowledge about the solution is a fundamental requirement in illposed reconstruction, and such knowledge is often found valuable even in well-posed

computations.

Having located a successful pair (ω, q), define |λQ| = ω

1/(1−q)

. In the well-posed

forward problem, one may then use Eq. (3.14) to estimate the stabilization penalty

as ∆t ↓ 0, in the successfully completed computation. However, that penalty is the

one that would have been incurred with the linear representative operator L, rather

than the actual nonlinear operator Le.

In the ill-posed backward problem, interactive adjustment of the parameter pair

(ω, q) based on prior knowledge, is even more crucial in obtaining useful reconstructions as the values of M and δ in Eq. (3.16) are seldom known accurately. Again,

after locating an optimal pair, define |λQ| = ω

1/(1−q)

. One may then use the second

term in Eq. (3.18) to evaluate the stability penalty deviation from the best-possible

result as ∆t ↓ 0. As before, this is the deviation that would have been incurred with

L rather than the actual Le.

4. Laplacian stabilization of the explicit scheme for ut = −Lu. Closed

form expressions for the eigenfunctions of the Laplacian are known for specific domains

that are important in applications, including rectangles, circles, and spheres [11], [16].

On such domains, the use of smoothing operators based on the Laplacian, rather than

L

∗L, to stabilize the explicit scheme for ut = −Lu, is of considerable interest. Such

stabilization is feasible for parabolic problems ut = −Lu for which the inequality in

Eq. (2.4) holds.

With fixed constants Γ, , ∆t, and p, 1 < p ≤ 3, let R and B denote the linear

operators in L

2

(Ω) defined by

R = exp{−∆t(−∆)p(4.1) }, B = R(I − ∆tL)/Γ.

As in Eq. (3.20), let {σm, ψm}

K

m=1 denote a sufficiently large set of known eigenpairs

of (−∆) on Ω, with the same boundary conditions on ∂Ω as prescribed for L. Let the

subspace VK ⊂ L

2

(Ω) spanned by {ψm}

K

m=1 be large enough to well-approximate the

solutions of ut = −Lu. We may then approximate the compensating operator R on

VK by

Rv =

X

K

m=1

exp{−∆t|σm|

p

}

Z

Ω

v(x)ψm(x)dx(4.2) ψm(x), v ∈ VK.
STABLE EXPLICIT TIME-MARCHING IN NONLINEAR PARABOLIC EQUATIONS 9

Using R and B leads to a new compensated explicit scheme for ut = −Lu, u(0) =

f, on 0 ≤ t ≤ T,

v

0 = f,

v

n+1 = Bvn ≡ R{vn − ∆tLvn}/Γ, n = 0, 1, 2, · · · , NT − 1, ∆t = T /NT .

(4.3)

Lemma 4. Let q, Q, ω, be as in Lemma 1, and let R and B be the operators in

Eq. (4.1). Choose Γ,  > 0, and p ≥ q, such that for g ∈ L

2

(Ω)

k exp{−∆t(−∆)p}g k2≤ Γ k exp{−ω∆t(L

∗L)q/2

(4.4) }g k2,

as postulated in Eq. (2.4). Then, k B k2≤ 1 + ∆t|λQ|, the compensated explicit

scheme Eq. (4.3) is unconditionally stable, and v

n satisfies

k v

n

(4.5) k2≤ exp{|λQ|n∆t} k f k2, n = 0, 1, 2, · · · , NT .

Proof : From Eq. (4.4) and Lemma 2, with A as in Eq. (3.1),

k v

n+1 k2 = k Bvn

k2≤k S(v

n − ∆tLvn

) k2=k Avnk2

(4.6) ≤ exp{|λQ|(n + 1)∆t} k f k2, n = 0, 1, 2, · · · , NT .

QED. Using the analog of Lemma 3 for the compensating operator R, together with

the argument in Theorem 1, leads to the corresponding result for the explicit scheme

in Eq. (4.3).

Theorem 2. Let u

n be the unique solution of ut = −Lu, u(0) = f, at t = n∆t.

Let v

n be the corresponding solution of the R-compensated scheme in Eq. (4.3), and

let ER(t) ≡ u

n − vn, denote the error at t = n∆t, n = 0, 1, 2, · · · , NT . Then, with

p, q, Q, Γ,  as in Lemma 4, and ω = |λQ|

1−q

,

k ER(t) k2 ≤ e

|λQ|t

k ER(0) k2 +

n

(e

|λQ|t − 1)/|λQ|

o

×



 |||(−∆)pu|||2,∞ + ∆t |||(−∆)pLu|||2,∞ + 1/2(∆t) |||L

2u|||2,∞

	

= e

|λQ|t

k ER(0) k2 +{(e

|λQ|t − 1)/|λQ|q

} (/ω) |||(−∆)p(4.7) u|||2,∞ + O(∆t).

Analogously to Eqs. (3.14), (3.18), we have the following Corollaries to Theorem 2.

Corollary 1. In Theorem 2, let ut = −Lu, u(0) = f, be a well-posed forward

problem with data f known exactly. Then, Eq. (4.7) reduces to

k ER(t) k2≤ {(e

|λQ|t − 1)/|λQ|q

} (/ω) |||(−∆)p(4.8) u|||2,∞ + O(∆t).

Corollary 2. In Theorem 2, let ut = −Lu, u(0) = f, 0 ≤ t ≤ T, be an

ill-posed backward problem with approximate input data f and solution u(t) satisfying

(4.9) k ES(0) k2≤ δ, k u(., T) k2≤ M,
10 A. S. CARASSO

with known M, δ. Then, choosing |λQ| = (1/T) log(M/δ), Eq. (4.7) reduces to

k ER(t) k2 ≤ Mt/T δ

(T −t)/T

+ {(e

|λQ|t − 1)/|λQ|q

} (/ω) |||(−∆)p(4.10) u|||2,∞ + O(∆t), 0 ≤ t ≤ T.

Remark 2. Laplacian stabilization and the operator R requires the selection of

 > 0, p > 1 in Eq. (4.2), and Γ > 0 in Eq. (4.3), together with subsequent interactive adjustment of these parameters based on prior knowledge about the solution.

Having located a useful triple (, p, Γ), an estimate of the stabilization penalty as

∆t ↓ 0 is not immediate in Eqs. (4.8), (4.10). Knowledge of the key parameters ω

and q in the governing inequality in Eq. (2.4) is needed to define |λQ| = ω

1/(1−q)

,

and evaluate the necessary terms in Eqs. (4.8), (4.10).

5. Fourier space Laplacian stabilization on rectangular regions. For simplicity of exposition we restrict attention to R2 and functions h(x, y), but the discussion is valid in higher dimensions. Define the 2D Fourier transform of an L

1

(R2)

function h(x, y) by

hˆ(ξ, η) ≡

Z

R2

h(x, y)e

−2πi(ξx+ηy)

(5.1) dxdy.

For fixed p > 1, and small  > 0, consider the well-posed forward diffusion initial

value problem on the whole space R2, involving possibly non-integer powers of the

negative Laplacian,

wt = −(−∆)p

(5.2) w, t > 0, w(x, y, 0) = h(x, y).

This has the unique Fourier space solution

wˆ(ξ, , t) = exp{−t[(2πξ)

2 + (2πη)2

]

p

}hˆ (5.3) (ξ, η), t > 0,

from which w(x, y, t) can be found by inverse Fourier transformation

w(x, y, t) = Z

R2

exp{2πi(ξx + ηy)} exp{−t[(2πξ)

2 + (2πη)2

]

p

}hˆ (5.4) (ξ, η)dξdη.

Next, consider the parabolic initial value problem ut = −Lu, t > 0, u(0) = f, on all

of R2, with the explicit scheme u

n+1 = un−∆tLun, n ≥ 0. Analogously to Eq. (4.2),

we may define the compensating operator R acting on g(x, y) ∈ L

1

(R2) ∩ L

2

(R2) by

Rg =

Z

R2

exp{2πi(ξx + ηy)} exp{−∆t[(2πξ)

2 + (2πη)2

]

p

(5.5) }gˆ(ξ, η)dξdη.

Assume the inequality in Eq. (2.4) holds for g(x, y) ∈ L

1

(R2) ∩ L

2

(R2). With Γ > 0

as in Eq. (2.4), this leads to the compensated explicit scheme on 0 ≤ t ≤ T,

v

0 = f,

v

n+1 = R{vn − ∆tLvn}/Γ, n = 0, 1, 2, · · · , NT − 1, ∆t = T /NT .

(5.6)
 

STABLE EXPLICIT TIME-MARCHING IN NONLINEAR PARABOLIC EQUATIONS 11

5.1. FFT periodized Laplacian stabilization and IM ×IM images. Consider now the parabolic problem ut = −Lu, t > 0, u(0) = f, defined on the unit

square Ω in R2, with u(x, y,t) an 8-bit IM × IM pixel image. We may discretize

the elliptic spatial operator L using centered finite differencing, with homogeneous

Dirichlet or Neumann conditions on ∂Ω. With pre-selected Γ, E, p, we then apply

the compensated explicit scheme in Eq. (5.6) above as follows. At each time step,

direct and inverse Fast Fourier Transform (FFT) algorithms are used to synthesize

Rg in Eq. (5.5), where g is the IM × IM pixel image g = (vn − ΔtLvn}/Γ. However, the FFT algorithm assumes g to be extended by periodicity to all of R2, and n+1 returns an IM ×IM array v = Rg satisfying periodic boundary conditions. At the

next time step, application of the discretized operator L to vn+1 restores the original

Dirichlet or Neumann boundary conditions. In practice, in many problems of interest,

such alternating erroneous boundary conditions are found to cause spurious artifacts

at the very edges of the image, without impairing the results away from the edges.

Such incovenience is a small price to pay for the highly efficient FFT synthesis of the

compensator R, which enables consideration of 512 × 512 or 1024 × 1024 arrays. This

technique can clearly be extended to nonlinear parabolic equations on hyperrectangles

in Rn.

6. Van Cittert iteration and backward continuation in ut = −Lu. As

originally formulated [24], the Van Cittert iterative procedure was intended for solving

1D linear convolution integral equations Pw = f with explicitly known kernels, by

means of the iterative procedure

(6.1) wm+1(x) = wm(x) + λ {f(x) − P[wm(x)]} , m ≥ 1.

Here, λ > 0 is a fixed relaxation parameter, w1(x) = λf(x), and the expectation is m that w → w. In reality, ill-posedness of the integral equation leads to error amplification at high frequencies which prevents convergence. However, in spectroscopy

and image processing applications [9], [12], the Van Cittert method is often found to

produce useful results after finitely many iterations, even though it may not converge.

Application of the Van Cittert iteration to nonlinear parabolic backward continuation is discussed in [7]. Consider the well-posed forward nonlinear problem

ut = −Lu, u(0) = f, 0 ≤ t ≤ T, in a domain Ω ∈ Rn. Let ΛT be the nonlinear

solution operator at time T so that ΛT [u(0)] = u(T). The Van Cittert procedure

would seek to recover the initial value f = u(0), from approximate data g for u(T),

through the iteration

m+1 (6.2) w = wm + λ g − ΛT [wm] , m ≥ 1,

1 with w = λg, and fixed relaxation parameter λ, 0 < λ < 1. By discretizing the

forward problem ut = −Lu, u(0) = f, 0 ≤ t ≤ T, using a stable computational

method, we may define a corresponding discrete nonlinear solution operator ΛT , and d

use it in Eq. (6.2) to implement the procedure computationally. Several examples of

successful nonlinear backward recovery using this technique are presented in [7].

7. 2D nonlinear computational examples. The illustrative examples below

involve two nonlinear second order parabolic equations on the unit square Ω ≡ {0 ≤

x, y ≤ 1} in R2, with homogeneous Neumann boundary conditions. The initial values

w(x, y, 0) = g(x, y), are 8 bit grey scale images with intensity values ranging from 0

to 255, leading to sizeable nonlinearities involving w. All images are of size 512 × 512
 

 

12 A. S. CARASSO

pixels, and are associated with complicated intensity surface plots such as are shown

in Figure 1.1. One equation has the form

t wt = κs(w)\.{q(x, y,t)\w} + c |w|wx + d(w cos2 w)wy, Ω × (0, T),

(7.1)

w(x, y, 0) = g(x, y),

where κ = 8.5 × 10−4, c > 0, d > 0, and

s(w) = exp(0.002w),

(7.2)

q(x, y,t) = exp(10t)(1 + 5 sin πx sin πy) ≥ 1, Ω × (0, T).

The other equat

t

ion used differs from the above by having the term cwwx in place

of the term c |w|wx . Clearly, the nonlinear elliptic operator LL in Eq. (7.1) is very

far from the linear L contemplated in the assumed inequality Eq. (2.4). Numerous

other interesting nonlinear parabolic equations remain to be explored.

Discretization of Eq. (7.1) uses centered finite differencing for the space variables

with Δx = Δy = 1/512, together with pure explicit time differencing, with Δt > 0

chosen so as to yield useful results.

8. Well-posed forward experiment. The first experiment illustrates the usefulness of the stabilized explicit scheme in forward solving Eq. (7.1) up to time

T = 0.004, with c = 0.0625, d = 0.075. With Δt = 0.0002, the uncompensated pure

explicit scheme applied to the 512×512 pixel Sydney and B2 Bomber images develops

instabilities, as shown in the middle images in Figure 8.1. However, applying FFT

Laplacian stabilization as discussed in Section 5, using Γ = 1.0, E = 4.0 × 10−6, p =

2.5, produces useful results with the same value of Δt. This is shown in the rightmost

images in Figure 8.1. This particular (E, p, Γ) triple was arrived at after very few

trials. Note the FFT-induced edge effects at top and bottom of the stabilized images.

For comparison, an uncompensated stable pure explicit calculation was also performed, with Δt chosen ten times smaller. Slightly different results were obtained at

T = 0.004. These differences amounted to an L1 relative error of 9.26% for the Sydney

image, and 3.87% for the B2 Bomber image. Such L1 errors indicate relatively small

stabilization penalties in these two examples.

It is noteworthy that the linear analysis provided in Sections 2 through 5, together

with the assumed inequality in Eq. (2.4), appear to be applicable in this nonlinear

case. The FFT compensated explicit scheme becomes ever more useful in well-posed

forward 3D problems on finer spatial meshes.

9. Ill-posed backward recovery experiments. While the forward computational results are encouraging, successful backward reconstruction in nonlinear problems is considerably more difficult, as the accuracy of the input data at time T

plays a crucial role. For linear autonomous selfadjoint forward parabolic equations

ut = −Lu, 0 ≤ t ≤ T, approximate data g at time T such that u(T) − g 2≤ δ,

for a solution u(t) known to satisfy u(0) 2≤ M, leads to an L2 uncertainty of

M(T −t)/T δt/T in backward reconstruction at time t < T. See [10], [15].

However, if L is nonlinear, the corresponding L2 uncertainty at time t < T may

be larger [10]. It is given by M1−μ(t)

δμ(t), where the H ̈older exponent μ(t) satisfies

0 ≤ μ(t) ≤ 1, with μ(t) > 0, t > 0, μ(T) = 1, μ(0) = 0, and μ(t) ↓ 0 monotonically

as t ↓ 0. Rapid exponential decay of μ(t to zero as t ↓ 0, is possible in some cases.
STABLE EXPLICIT TIME-MARCHING IN NONLINEAR PARABOLIC EQUATIONS 13

 STABILIZED EXPLICIT SCHEME IN FORWARD NONLINEAR PARABOLIC EQUATIONS

ORIGINAL SYDNEY IMAGE UNSTABLE PURE EXPLICIT SOLUTION STABILIZED EXPLICIT SOLUTION

ORIGINAL B2 BOMBER IMAGE UNSTABLE PURE EXPLICIT SOLUTION STABILIZED EXPLICIT SOLUTION Fig. 8.1. Useful computation of well-posed forward nonlinear parabolic problems on 512 × 512 mesh, using stabilized explicit scheme with a value of

Δt that

is unstable in the uncompensated pure explicit scheme. Stabilized results differ only slightly from uncompensated stable explicit computations with

Δt chosen

ten times smaller. Such differences, amounting to an

L1 relative error of 9.26% for the Sydney image, and 3.87% for the B2 Bomber image, indicate relatively

small stabilization penalties. Note FFT-induced edge effects at top and bottom of rightmost images.
14 A. S. CARASSO

For this reason, high quality nonlinear reconstructions from imprecise data are not

generally feasible. See [6].

In the examples below, the input data at time T is obtained by a preliminary forward solution of the nonlinear parabolic equation, using a stable numerical method.

That numerically obtained blurred image may not be a sufficiently accurate approximation to the true solution of Eqs. (7.1) or (9.1), at time T. In some cases, recovering

the finest details in the sharp image at t = 0 from such modestly accurate data at

t = T, may lie beyond the capabilities of any method on theoretical grounds alone.

Sagittal MRI brain scan. The original sharp MRI brain image in Figure 9.1

was blurred by using it as initial data in Eq. (7.1) with c = 0.05, d = 0.06, and

applying the stable uncompensated pure explicit scheme with Δt = 1.75 × 10−6, to

march forward 2000 time steps up to the final time T = 3.5 × 10−3. This produced

the blurred image in Figure 9.1. While this is a reasonably accurate computation, a

great deal of fine detail in the sharp image has been smoothed out by the nonlinear

evolution, as becomes evident at high magnification.

The numerically blurred image was used as input data in the FFT Laplacianstabilized explicit scheme discussed in Section 5, using Γ = 1.0, E = 2.0 × 10−8, p =

2.75. With Δt = 1.75 × 10−6, the compensated scheme was marched 2000 time steps

backward from T = 3.5 × 10−3, to produce the deblurred image shown at the bottom

left of Figure 9.1. Clearly, significant deblurring was achieved, and the linear analysis

in Sections 2 through 5, together with Eq. (2.4), again appear to be germane in this

nonlinear situation. Here, the small stabilization penalty is added to the nonlinear

L2 uncertainty, M1−μ(t)

δμ(t), with an unknown μ(t).

The Van Cittert iterative procedure discussed in Section 6 was succesfully applied

to several nonlinear examples in [7], although at substantially higher computational

costs. Significantly, this procedure fails in the present example, as shown at the

bottom right image in Figure 9.1. Inverse reconstruction methods generally work well

only in limited cases, and alternative approaches, when available, are valuable.

It is instructive to view the above reconstruction process as it operates on the actual image intensity data that is fed into the computational algorithm. This is shown

in Figure 9.2. Clearly, quite complex data surfaces are involved, and accurate recovery would seem quite difficult. Inspection of the bottom two reconstructed surfaces

in Figure 9.2 reveals how easily the process can go awry.

Elizabeth Taylor image. A familiar face image makes a valuable test example.

The original 512×512 pixel Elizabeth Taylor image in Figure 9.3 was blurred by using

it as initial data in the following equation on the unit square Ω

wt = κs(w)\.{q(x, y,t)\w} + cwwx + d(w cos2 w)wy, Ω × (0, T),

(9.1)

w(x, y, 0) = g(x, y),

with κ and q(x, y,t) as in Eq. (7.1), c = 0.05 and d = 0.06. The stable uncompensated

pure explicit scheme was applied with Δt = 7.0 × 10−6 and 600 time steps up to the

final time T = 4.2 × 10−3. With this larger Δt, the resulting numerical solution at

time T is a less accurate approximation to the true solution than was the case in

Figure 9.1.

The stronger nonlinearity associated with the wwx term leads to the striking

transformation shown in the resulting blurred image in Figure 9.3. Such distortion
STABLE EXPLICIT TIME-MARCHING IN NONLINEAR PARABOLIC EQUATIONS 15

ORIGINAL MRI BRAIN SCAN NONLINEARLY BLURRED IMAGE

 STABILIZED EXPLICIT DEBLUR VAN CITTERT ITERATIVE DEBLUR

STABILIZED EXPLICIT SCHEME IN BACKWARD

 NONLINEAR PARABOLIC EQUATIONS

Fig. 9.1. Backward recovery from nonlinearly blurred data in 512×512 MRI brain image, using

stabilized explicit scheme with Γ = 1.0, E = 2.0 × 10−8, p = 2.75. Van Cittert procedure develops

instabilities after 30 iterations.
16 A. S. CARASSO

ORIGINAL MRI BRAIN SCAN NONLINEARLY BLURRED IMAGE

 STABILIZED EXPLICIT DEBLUR VAN CITTERT ITERATIVE DEBLUR

STABILIZED EXPLICIT RECONSTRUCTION OF

 MRI BRAIN SCAN INTENSITY SURFACE

Fig. 9.2. Reconstruction of highly irregular original intensity surface from blurred data in MRI

brain image, using stabilized explicit scheme and Van Cittert procedure.
STABLE EXPLICIT TIME-MARCHING IN NONLINEAR PARABOLIC EQUATIONS 17

ORIGINAL LIZ TAYLOR IMAGE NONLINEARLY BLURRED IMAGE

STABILIZED EXPLICIT DEBLUR VAN CITTERT ITERATIVE DEBLUR

STABILIZED EXPLICIT SCHEME IN BACKWARD

 NONLINEAR PARABOLIC EQUATIONS

Fig. 9.3. Backward recovery from nonlinearly blurred data in 512×512 Liz Taylor image, using

stabilized explicit scheme with Γ = 1.0, E = 7.5 × 10−8, p = 2.735. Van Cittert procedure develops

instabilities after 8 iterations

.
18 A. S. CARASSO

 ORIGINAL LIZ TAYLOR IMAGE NONLINEARLY BLURRED IMAGE

STABILIZED EXPLICIT DEBLUR VAN CITTERT ITERATIVE DEBLUR

 STABILIZED EXPLICIT RECONSTRUCTION OF

LIZ TAYLOR IMAGE INTENSITY SURFACE

Fig. 9.4. Reconstruction of highly irregular original intensity surface from blurred data in Liz

Taylor image, using stabilized explicit scheme and Van Cittert procedure.

.
STABLE EXPLICIT TIME-MARCHING IN NONLINEAR PARABOLIC EQUATIONS 19

is not easily achieved using traditional convolution with a Gaussian point spread

function. Recovery of the original from such a strongly blurred image, computed

with only modest accuracy, seems less feasible than was possible in the MRI image.

Surprisingly, using Γ = 1.0, E = 7.5 × 10−8, p = 2.735, and marching back 600 time

steps Δt = 7.0 × 10−6 from t = T in the stabilized explicit scheme, produced the

reasonably good reconstruction shown in Figure 9.3. This is a significant example.

Clearly, the stabilization penalty must be quite small to permit such useful recovery

in the presence of the nonlinear L2 uncertainty M1−μ(t)

δμ(t).

The poor performance of the Van Cittert iteration in Figure 9.3 is a reflection of

the particularly difficult backward continuation problem associated with the blurred

Elizabeth Taylor image, due to the nonlinearity in Eq. (9.1). The reconstruction

process on the corresponding intensity surface data is displayed in Figure 9.4.

10. Relation to the quasi-reversibility method. With its emphasis on solving multidimensional nonlinear problems using explicit schemes, the present method

is a significant new approach that complements the quasi-reversibility (QR) method

[13]. Applications of the QR method to backward parabolic equations are developed in [13, Chapter 1]. There, linear equations in the form wt + A(t)w = 0, are

considered, with A(t) a linear elliptic partial differential operator of order 2m with

coefficients depending smoothly on space and time, and such that the forward problem is well-posed. Backward continuation from given noisy data at time T > 0, is

accomplished by numerically marching back from t = T in the modified equation

wt +A(t)w −ωA∗(t)A(t)w = 0, with suitably pre-selected small ω > 0. This modified

equation is well-posed backward in time, and a variety of stable difference schemes may

be used. As formulated in [13, Chapter 1], the QR method is not applicable to nonlinear problems in an obvious way, due to the difficulty of defining a suitable adjoint

operator A∗(t). Indeed, linear problems are the primary focus in [13]. In addition,

implicit Crank-Nicolson time differencing is contemplated for the numerical solution

of the modified higher order equation involving A∗(t)A(t). In multidimensional problems, this requires computationally intensive solution of the algebraic problem at each

time step. This is avoided in the present method.

A valuable feature of the stabilized explicit scheme is that it allows for efficient

and simultaneous exploration of the parameter values (Γ, E, p), including non integer

positive values of p in (−Δ)p, all within the same computational code. In the QR

method for second order parabolic equations, the modified equation involves a fixed

fourth order stabilizing elliptic operator, and there is no option for exploring other

stabilizers. In this connection, it is noteworthy that in both the well-posed and illposed test examples in Sections 8 and 9, values of p such that 2.5 ≤ p ≤ 2.75, were

arrived at interactively. Apparently, a compensating pseudo-diffusion equation of

order five or higher, seems to be necessary to control explicit scheme instability in

computing second order nonlinear parabolic equations.

11. Concluding Remarks. The unexpected success of the FFT Laplacianstabilized explicit scheme on 2D nonlinear parabolic problems in rectangular regions,

indicates that the use of Laplace stabilization in other special domains in Rn where

Laplacian eigenfunctions are known in closed form, is worthy of active consideration.

More generally, in arbitrary domains Ω ∈ Rn as envisaged in Eq. (3.20), preliminary

computation of a sufficient number of the characteristic pairs of a variable coefficient

linear elliptic operator L, believed to be representative of the given nonlinear elliptic operator LL, may provide the key to useful stabilized explicit solution of difficult

nonlinear parabolic problems.
20 A. S. CARASSO

The apparent validity of Eq. (2.4) in the numerical experiments in Sections 8 and

9, raises interesting analytical research problems.

The failure of the Van Cittert method in the two examples in Section 9, despite

its successful use in [7], indicates a continuing need for alternative reconstruction

procedures, based on different analytical strategies.

The use of fictitious images in exploring 2D nonlinear parabolic equations, together with displays of the corresponding intensity surfaces, was particulary instructive.

Future work will explore the use of compensated explicit schemes in non-parabolic

irreversible systems occurring in viscoelasticity, thermoelasticity, coupled sound and

heat flow, and other areas of application discussed in [13].

REFERENCES

[1]	 W. Arendt and A. F. M. ter Elst, Gaussian estimates for second order elliptic operators

with boundary conditions, J. Operator Theory 38 (1997), pp. 87-130.

[2]	 D. G. Aronson, Bounds for the fundamental solution of a parabolic equation, Bull. Amer.

Math. Soc. 73 (1967), pp. 890-896.

[3]	 J. Atmadja and A. C. Bagtzoglou, State of the art report on mathematical methods for

groundwater pollution source identification, Environmental Forensics 2 (2001), pp. 205­

214.

[4]	 J. Atmadja and A. C. Bagtzoglou, Pollution source identification in heterogeneous porous

media, Water Resources Research 37 (2001), p. 2113-2125. DOI 10.1029/2001WR000223

[5]	 I. Babuˇska and J. E. Osborn, Finite element-Galerkin approximation of the eigenvalues and

eigenvectors of selfadjoint problems, Math. Comp. 52 (1989), pp. 275-297.

[6]	 A. S. Carasso, Reconstructing the past from imperfect knowledge of the present: Effective

non-uniqueness in solving parabolic equations backward in time, Mathematical Methods in

the Applied Sciences (2012), DOI: 10.1002/mma.2582

[7]	 A. S. Carasso, Compensating operators and stable backward in time marching in nonlinear

parabolic equations, Int J Geomath 5 (2014), pp. 1–16.

[8]	 V. Heuveline, On the computation of a very large number of eigenvalues for selfadjoint elliptic

operators by means of multigrid methods, Journal of Computational Physics 184 (2003),

pp. 321-337.

[9]	 P. A. Jansson, Deconvolution with Application in Spectroscopy, Academic Press, New York,

1984.

[10]	 R. J. Knops, Logarithmic convexity and other techniques applied to problems in continuum

mechanics, in Symposium on Non-Well-Posed Problems and Logarithmic Convexity, R. J.

Knops, ed., Lecture Notes in Math. 316, Springer-Verlag, New York, 1973.

[11]	 J. R. Kuttler and V. G. Sigillito, Eigenvalues of the Laplacian in two dimensions, SIAM

Review 26 (1984), pp. 163–193.

[12]	 R. L. Lagendijk and J. Biemond, Iterative Identification and Restoration of Images, Kluwer

Academic Publishers, Boston, 1991.

[13] R. Lattes` and J. L. Lions, The Method of Quasi-Reversibility, American Elsevier, New York,

1969.

[14] J. D. Logan, Transport Modeling in Hydrogeochemical Systems, Springer, New York, 2001.

[15]	 K. Miller, Least squares methods for ill-posed problems with a prescribed bound, SIAM J.

Math. Anal. 1 (1970), pp. 52–74.

[16]	 P. M. Morse and H. Feshbach, Methods of Theoretical Physics, Mc Graw-Hill, New York,

(1953).

[17]	 R. M. Neupauer, B. Borchers, and J. L. Wilson, Comparison of inverse methods of reconstructing the release history of a groundwater contamination source, Water Resources

Research 36 (2000), pp. 2469–2475. DOI 10.1029/95WR02383.

[18]	 E. M. Ouhabaz, Gaussian estimates and holomorphy of semigroups, Proc. Amer. Math. Soc.

123 (1995), pp. 1465-1474.

[19]	 E. M. Ouhabaz, Gaussian upper bounds for heat kernels of second order elliptic operators with

complex coefficients on arbitrary domains, J. Operator Theory 51 (2004), pp. 335-360.
STABLE EXPLICIT TIME-MARCHING IN NONLINEAR PARABOLIC EQUATIONS 21

[20]	 R. B. Platte and T. A. Driscoll, Computing eigenmodes of elliptic operators using radial

basis functions, Computers and Mathematics with Applications 48 (2004), pp. 561-576.

[21]	 R.D. Richtmyer and K. W. Morton, Difference Methods for Initial Value Problems, Second

Edition, Wiley, New York, 1967.

[22]	 R. N. Singh, Advection diffusion equation models in near-surface geophysical and environmental sciences, J. Ind. Geophys. Union 17 (2013), pp. 117–127.

[23]	 T. H. Skaggs and Z. J. Kabala, Recovering the history of a groundwater contaminant plume:

Method of Quasi-Reversibility, Water Resources Research, 31 (1995), pp. 2669–2673. DOI

10.1029/95WR02383.

[24]	 P. H. Van Cittert, Zum Einfluss der Spaltbreite auf die Intensittsverteilung in Spektrallinien

II, Z. Physik 69 (1931), pp. 298-308.