NISTIR 7967
 

Compensating Operators and 

Stable Backward in Time Marching 

in Nonlinear Parabolic Equations 

Alfred S. Carasso 
 

 

 

 

 

 

NISTIR 7967 

Compensating Operators and 

Stable Backward in Time Marching 

in Nonlinear Parabolic Equations 

Alfred S. Carasso 

Applied and Computational Mathematics Division 

Information Technology Laboratory 

November 2013 

U.S. Department of Commerce 

 Penny Pritzker, Secretary 

National Institute of Standards and Technology 

Patrick D. Gallagher, Under Secretary of Commerce for Standards and Technology and Director 
COMPENSATING OPERATORS AND STABLE BACKWARD IN

TIME MARCHING IN NONLINEAR PARABOLIC EQUATIONS.

ALFRED S. CARASSO∗

Abstract. Step by step time-marching schemes are fundamental tools in the numerical exploration of well-posed nonlinear evolutionary partial differential equations. However, when the initial

value problem is ill-posed, such stepwise numerical schemes are necessarily unconditionally unstable

and result in explosive noise amplification. This paper outlines a novel stabilized time-marching

procedure for computing nonlinear parabolic equations on 2D rectangular regions, backward in time.

Very little is known either analytically, or computationally, about this class of exponentially ill-posed

problems. To quench the instability, the procedure uses easily synthesized FFT-based compensating

operators at every time step. A fictitious nonlinear image deblurring problem is used to evaluate

the effectiveness of this computational approach. The method is compared with a previously introduced, global in time nonlinear Van Cittert iterative procedure. The latter is significantly more time

consuming, and impractical on large problems.

Key words. nonlinear backward parabolic equations; ill-posed initial value problem; backward

time-marching scheme; FFT compensating operators; non-integer power Laplacian; nonlinear Van

Cittert method; nonlinear image deblurring.

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

1. Introduction. Reconstructing the past behavior of a physical system from

knowledge of its current state, is of interest in numerous scientific contexts. However, when the physical process involves diffusion, as is the case in groundwater contamination studies for example, the resulting backward parabolic reconstructions are

exponentially ill-posed and present serious computational issues [1], [2], [3]. Similar

exponential ill-posedness characterizes elliptic downward continuation of satellite measurements of gravity and geomagnetism [4], [5], [6]. In recent years, significant progress

has been made in successfully solving linear time-reversed parabolic equations, leading to useful reconstructions of valuable scientific data in astronomy, medicine, and

other areas of application [7]. However, little is known about backward continuation

in multidimensional nonlinear parabolic equations, either analytically or computationally.

This paper outlines a novel time-marching procedure for computing time-reversed

nonlinear parabolic equations on 2D rectangular regions. As is well-known, stepwise marching schemes consistent with ill-posed initial value problems are necessarily

computationally unstable. Here, easily synthesized compensating operators, based on

FFT algorithms, are applied at every time step to quench the instability. A fictitious

image deblurring problem is used to evaluate the effectiveness of this computational

approach, as well as to probe the behavior of various types of nonlinearities. Such canary in the coalmine synthetically blurred images are useful indicators of the type and

quality of information that might be recoverable, when the backward computation is

applied to real data in engineering contexts unrelated to imaging. Further exploration

and refinement of this time-marching methodology will be developed in subsequent

reports. Similar ideas may be applicable in ill-posed initial value problems associated

with other types of partial differential equations.

Recently [8], a global in time iterative procedure was successfully applied to a class

of 2D nonlinear backward parabolic equations. This method is an adaptation to the

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

1
2 A. S. CARASSO

nonlinear case of the Van Cittert iteration [9], originally used in spectroscopy in the

1930’s to solve 1D linear convolution integral equations, with explicitly known kernels

having positive Fourier transforms. It is remarkable that the Van Cittert method

works as well as it does on the limited, but significant class of nonlinear problems

considered in [8]. However, the nonlinear Van Cittert process becomes prohibitively

time consuming on large problems, such as 512 × 512 or 1024 × 1024 pixel images.

As will be shown below, when backward continuation is feasible, the new stepwise

marching scheme can produce results with almost the same visual quality as in the

Van Cittert iteration, but almost 80 times faster.

2. Limited backward recovery in nonlinear parabolic problems. Most of

the useful computational results in backward in time parabolic continuation have generally been obtained in the canonical case of linear, autonomous, selfadjoint problems.

As shown in [3], [8], nonlinear parabolic equations deviating too strongly from this

canonical model can present additional challenges. The reason for this marked difference in behavior can be traced back to the H ̈older-continuity stability inequalities

that underlie backward reconstruction from noisy data.

Let Ω be a bounded domain in Rn with smooth boundary ∂Ω. Let L be a linear or

nonlinear elliptic operator in Ω, acting on smoothly differentiable functions satisfying

homogeneous Dirichlet or Neumannn conditions on ∂Ω. Let L be such that the

forward initial value problem wt = Lw, t > 0, w(0) = g(x), is well-posed in L2(Ω).

Let w1(x,t) and w2(x,t) be any two solutions, and let F(t) = w1(.,t) − w2(.,t) 2

2,

for 0 ≤ t ≤ T.

Using logarithmic convexity techniques [10], the folowing sharp inequality can be

established for a wide class of parabolic equations wt = Lw,

(2.1) F(t) ≤ {F(0)}1−μ(t)

{F(T)}μ(t), 0 ≤ t ≤ T.

Here, 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.

Given noisy data f(x) ∈ L2(Ω) at time T, satisfying w(., T)−f 2≤ δ for small

δ > 0, we restrict attention to backward continuations w(x,t) satisfying a prescribed

bound at t = 0, i.e., w(., 0) 2≤ M. Eq. (2.1) then implies the following stability

inequality for the difference of any two possible continuations w1(x,t), w2(x,t), from

the data f(x) at t = T,

(2.2) w1(.,t) − w2(.,t) 2≤ 2M1−μ(t)

δμ(t), 0 ≤ t ≤ T.

In the linear autonomous selfadjoint case, the H ̈older exponent μ(t) in Eq. (2.2) decays

linearly to zero, as t ↓ 0. In the nonlinear case, rapid exponential decay of μ(t) to zero

is possible [3], [10]. Therefore, a level of accuracy δ at time T that can produce useful

backward reconstruction in an autonomous selfadjoint equation, may be insufficient in

a nonlinear equation. As a rule, high quality nonlinear reconstructions from imprecise

data are not easily achieved. Examples of physically plausible, smooth, positive, yet

false nonlinear reconstructions are given in [3]. In applications involving real data with

limited but unknown precision, detailed prior information about the true solution may

be critical in evaluating the merits of a candidate reconstruction.

3. Two nonlinear parabolic examples wt = Lw. The computational experiments discussed below involve four images and two instructive parabolic equations.

However, numerous other equations can be studied, exhibiting a wide variety of unexpected phenomena. Let Ω be the unit square 0 < x, y < 1 in the (x, y) plane. With
STABLE TIME-MARCHING IN ILL-POSED EVOLUTION EQUATIONS 3

fixed T > 0, and homogeneous Neumann boundary conditions on ∂Ω, the following

initial value problem will be studied,

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

(3.1)

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

Here κ = 8.5 × 10−4, c = 1.25, d = 0.6, and

s(w) = exp(0.025w),

(3.2)

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

An equation with different nonlinearities will also be considered. This is

 wt = γr(w)\.{q(x, y,t)\w} + a |w|wx, +b(w cos2 w)wy, Ω × (0, T),

(3.3)

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

where γ = 0.04, a = b = 0.5, r(w) = exp{cos2(w)}, and q(x, y,t) is as previously

defined in Eq. (3.2).

Each of Eqs.(3.1), (3.3), is well posed in L2(Ω). Accordingly, given any initial

value w(x, y, 0) = g(x, y) ∈ L2(Ω), a unique solution w(x, y, T) exists at time T, and

the solution operator ΛT , where

(3.4) ΛT w(x, y, 0) = w(x, y, T),

is well-defined on L2(Ω). The nonlinear operator ΛT is not known explicitly. Rather,

ΛT w(x, y, 0) must be found by solving the appropriate initial value problem Eq.(3.1),

or Eq.(3.3), and obtaining the corresponding solution at time T. Note that w(x, y, T)

necessarily belongs to a very restricted class of smooth functions.

4. Explicit scheme for computing wt = Lw forward in time. Step by step

time marching explicit finite difference schemes are simple, convenient, and effective

tools for solving nonlinear evolutionary partial differential equations, as they allow

lagging the nonlinearities at the previous time step. Despite the necessary stability

restriction on the time step Δt, such schemes enable a user to quickly gain a preliminary indication of the salient characteristics of the unknown solution. Higher

precison numerics, along with unconditionally stable implicit schemes, may be applied at a subsequent stage. Here, using explicit time differencing and centered space

differencing in the nonlinear initial value problems in Eqs.(3.1) and (3.3), leads to

g i modest O Δt + (Δx)2 + (Δy)2 accuracy. However, the explicit stability condition

g i on Δt improves that accuracy to O (Δx)2 + (Δy)2 . This paper deals with 8 bit

gray scale IM ×IM pixel images, with pixel values ranging between 0 and 255. With

Δx = Δy = 1/IM, Δt = T /NT, the following difference approximation is used to

march the discrete mesh function Wn ≡ W(jΔx, kΔy, nΔt) in Eq. (3.1), NT time

steps Δt forward in time, up to time T = NT Δt,

Wn+1 = Wn + ΔtκS(Wn)\.{Qn\Wn} + cWnWn

x

+ d(Wn cos2 Wn)Wn, n = 0, NT − 1, y

(4.1) W0 = g(x, y).
�

4 A. S. CARASSO

Similar differencing is applied to Eq. (3.3). Homogeneous Neumann conditions are

applied on the boundary of the unit square. In this notation, Wx, Wy, \.{Q\W},

represent centered finite difference approximations in the space variables, W0 denotes

the original sharp image g(x, y), while WNT is the nonlinearly blurred image f(x, y),

using either Eq.(3.1) or Eq.(3.3). Analogously to the analytic problem in Eq. (3.4),

we define the discrete nonlinear parabolic blurring operator ΛT by d

= WNT (4.2) ΛT

d W0 ⇐⇒ ΛTd g(x, y) = f(x, y).

This nonlinear operator is defined on any 8 bit gray scale IM × IM pixel image

g(x, y) = W0. Applying ΛT to that image simply means applying the above explicit d

scheme for NT time steps to W0, and acquiring the resulting array f(x, y) = WNT .

We stress that the blurred image f(x, y) so obtained is only an approximation to the

true solution w(x, y, T) in Eq.(3.1) or Eq.(3.3).

5. Backward continuation in wt = Lw and the Van Cittert iteration.

In its original formulation, given the data f(x) and the explicitly known 1D linear ˆ convolution integral operator P with Fourier transform P(ξ) > 0, the Van Cittert

method solves P g = f for the unknown g(x) by means of the iterative procedure

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

Here, λ > 0 is a fixed relaxation parameter chosen so that 1 − λPˆ(ξ) > 0, h1(x) =

λf(x), and the expectation is that hm → g. Indeed, taking Fourier Transforms in Eq.

(5.1), we find

⎛ ⎞ m

(5.2) hˆm+1(ξ) = ⎝ {1 − λPˆ(ξ)}j

⎠ λfˆ(ξ), m ≥ 1.

j=0

 ∞ Since, |1 − λPˆ(ξ)| < 1, one has {(1 − λPˆ(ξ)})j = {λPˆ(ξ)}−1. Hence, hˆm(ξ) j=0

‘converges’ to gˆ(ξ) as m → ∞, and hm(x) → g(x). In reality, noise in the data f(x), ˆ together with rapid decay of P(ξ) at infinity in commonly occurring Gaussian-like

convolution operators P, lead to error amplification at high frequencies which prevents

convergence. The advantage of the Van Cittert iteration is that it provides a gradual

deconvolution process whereby low and moderately high frequency information can

be recovered in the first few iterations, while many more iterations are necessary to

obtain the remaining high frequency information. However, since such high frequency

information is quite likely contaminated by amplified data noise, it is wise to monitor

the iterative process and be able to terminate it before amplified noise begins to

degrade the reconstruction. In spectroscopy and image processing applications [11],

[12], the Van Cittert method is generally found to produce useful results after finitely

many iterations, even though it may not converge.

We now contemplate applying the Van Cittert iteration in the present parabolic

context to recover w(x, y, 0) = g(x, y) in Eqs.(3.1) and (3.3), given approximate values

f(x, y) for the true solution w(x, y, T) at time T. This amounts to inverting the

discrete nonlinear operator ΛT in Eq. (4.2) by solving for the unknown g(x, y) using d

the iterative procedure

 (5.3) hm+1(x, y) = hm(x, y) + λ f(x, y) − Λd

T hm(x, y) , m ≥ 1,

with some fixed λ such that 0 < λ < 1, and h1(x, y) = λf(x, y). While we may

hope that hm(x, y) → g(x, y), there are considerable theoretical and computational
 

STABLE TIME-MARCHING IN ILL-POSED EVOLUTION EQUATIONS 5

gaps between the iterative process in Eq. (5.1) and the one in Eq. (5.3). Unlike the

fixed, explicitly known, linear convolution operator P in Eq. (5.1), ΛT in Eq. (5.3) d

is unknown, and a new nonlinear initial value must be solved numerically to obtain

ΛT hm(x, y) for each m. In addition, the Neumann series expansion in Eq. (5.2) which d

underlies ‘convergence’ in Eq. (5.1), is not applicable in Eq. (5.3). Thus, the iterative

process in Eq. (5.3) is liable to be computationally intensive, yet with little or no

theoretical expectation of success. Nevertheless, remarkably, the Van Cittert iteration

is found to be a valuable tool in a wide variety of 2D nonlinear backward parabolic

equations. In many cases, this procedure generates iterates hm(x, y) such that the

L∞ norm of the residual, f −ΛT hm ∞, decays quasi-monotonically to a reasonably d

small value after a finite number N of iterations, and hN (x, y) is found to be a useful

approximation to w(x, y, 0).

6. FFT compensating operator B on rectangular regions. Define the 2D

Fourier transform of an L1(R2) function h(x, y) by

 

(6.1) hˆ(ξ, η) ≡ h(x, y)e−2πi(ξx+ηy)

dxdy,

R2

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

in R2 involving possibly non-integer powers of the negative Laplacian

(6.2) wt = −E(−Δ)pw, t > 0, w(x, y, 0) = h(x, y)

This has the unique Fourier space solution

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

2 + (2πη)2]p}hˆ(ξ, η), t > 0,

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

 

(6.4) w(x, y,t) = exp{2πi(ξx + ηy)} exp{−Et[(2πξ)

2 + (2πη)2]p}hˆ(ξ, η)dξdη.

R2

Now, let h(x, y) be an 8-bit IM × IM pixel image. Using direct and inverse FFT

algorithms, we can mimick the operations in Eqs. (6.3), (6.4), and apply the diffusion

equation Eq.(6.2) to the image h(x, y), to produce the smoothed image w(x, y,t). For

fixed t > 0, such smoothing will depend strongly on the chosen values of E and p in

Eqs. (6.3), (6.4).

Let Δt be the fixed time step in a step by step finite difference marching computation involving an IM × IM image. Let h(x, y) be the resulting image at some

discrete time tk = kΔt in the evolving computation. Using FFT algorithms, for given

E, p, we may apply the compensating linear operator B(E, p, Δt) to h(x, y), where Bh

is defined by

(6.5) B(E, p, Δt)h(x, y) = w(x, y, Δt),

with w(x, y,t) as in in Eq. (6.4).

Compensating operators in non rectangular regions, based on (−Δ)p, are discussed in Section 9 below.

7. Explicit scheme for computing wt = Lw backward in time. Given

the blurred image f(x, y) = ΛT

d g(x, y) in Eq. (4.2), we can attempt to recover the

original sharp image g(x, y) by marching backward from T = NT Δt in the explicit
6 A. S. CARASSO

finite difference scheme in Eq. (4.1). This becomes possible provided we apply a

properly chosen FFT compensating operator B to the computed image V at each

time step, to control the ill-posedness. With pre-assigned E, p, the proposed explicit

backward in time scheme corresponding to Eq. (4.1) is then the following

V = W m − ΔtκS(W m)\.{QNT −m\W m} + cW mW m

x

+ d(W m cos2 W m)W m, m = 0, NT − 1, y
 

W m+1
 = BV, m = 0, NT − 1.

(7.1) W0 = f(x, y).

Here, W0 is the given blurred image at time T = NT Δt, and the scheme marches

backward NT time steps to produce the sharp image WNT at t = 0. With obvious

modifications, this can also be applied in backward in time continuation in Eq.(3.3).

As may be expected, the choices of p and E are critical in obtaining good results, and

several interactive trials may be necessary. However, there are infinitely many useful

pairs (E, p), and such a pair can usually be located fairly quickly.

A different type of marching scheme for a class of ill-posed initial value problems

is discussed in [13].

7.1. Relation to the quasi-reversibility method. While close in spirit to

the quasi-reversibilty (QR) method [14], the compensated scheme in Eq. (7.1) differs

from QR in essential ways that render it more versatile, and more easily applicable to

nonlinear problems. The QR method deals with various ill-posed problems in partial

differential equations, including spatial continuation in elliptic equations, and boundary control and other problems in parabolic and non-parabolic evolution equations.

Application of QR to backward in time parabolic equations is developed in Chapter

1 of [14]. There, linear equations

(7.2) wt + A(t)w = 0

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

the space variables, with coefficients depending smoothly on space and time. The

differential operator A(t) need not be selfadjoint, but is such that the forward in time

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

(7.3) wt + A(t)w − EA∗(t)A(t)w = 0

with suitably preselected small E > 0. For the simple heat equation wt = Δw,

the modified equation is wt = Δw + EΔ2w. For that simple case, the authors also

contemplate the use of alternative modified equations, such as wt = Δw + E(−Δ)mw

with preselected positive integer m. Illustrative 1D linear computational examples,

using Crank-Nicolson time differencing, are also given in [14, Chapter 1]. However,

nonlinear backward problems are not considered. Indeed, the formulation in Eq. (7.3)

involving the adjoint A∗(t), is not applicable to nonlinear problems in an obvious way.

In the case of the nonlinear second order parabolic equations (3.1), (3.3), a modified

equation involving some appropriate nonlinear generalization of A∗(t)A(t)w, would

result in a fairly complex fourth order nonlinear equation, and present a challenging

computational task.

In the QR approach, a particular modified equation is selected and translated

into a specific computational code, using various trial values of E > 0. Exploring
STABLE TIME-MARCHING IN ILL-POSED EVOLUTION EQUATIONS 7

other possible modified equations within the same code would not be practical, and

is not typically contemplated. Also, only positive integer powers of various partial

differential operators enter the modified equations. A significant practical advantage

of the compensated scheme in Eq. (7.1), or its counterpart in more general domains

described in Section 9 below, is the ability to simultaneously explore non integer

positive values of p, as well as positive values of E, within the same computational code.

Interestingly, in the examples below involving the nonlinear second order parabolic

equations (3.1), (3.3), relatively simple compensating operators B based on (−Δ)p

with 2 < p < 3, were found fully capable of controlling instability in Eq. (7.1).

8. Some examples. Figure 8.1 shows four sharp test images. These images

are used to probe the behavior of the blurring nonlinear partial differential equations

wt = Lw in Section 3, as well as to evaluate the merits of the proposed backward in

time continuation algorithms. Each of these images features areas of zero intensity

near the edges. Such zero padding is important in applications involving FFT algorithms. All computations described below were carried out on an AMD Athlon 64 X2

4880+, running Centos 5 32 bit.1 The computation times mentioned are intended to

illustrate the possible advantages, or disadvantages, of various procedures.

Failure of nonlinear reconstruction. We begin by observing that useful

backward reconstruction may not be feasible in the absence of sufficient precision

δ in the given data at time T. Figure 8.2 is particularly instructive in the context of continuation backward in time in nonlinear parabolic equations. Each of

the three 256 × 256 degraded images on the left hand side, was obtained by blurring the corresponding sharp image in Figure 8.1, using the nonlinear initial value

problem Eqs.(3.1), (3.2). This was done using the explicit scheme in Eq.(4.1), with

Δx = Δy = 1/256, Δt = 3.0×10−7, and marching the solution 400 steps Δt forward,

up to the final time T = 1.2 × 10−4. These degraded images do not appear blurred

beyond remediation. Indeed, in linear image deconvolution experiments, similarly

blurred imagery can often be significantly improved. Here, useful recovery using the

Van Cittert method was not possible. Indeed, the original blurred images on the left

are more informative than the proposed reconstructions ! This example illustrates

the remarks in Section 2 regarding the significance of exponentially decaying H ̈older

exponents μ(t) in the backward stabilty estimate Eq. (2.2). Evidently, despite the

small value of Δt in Eq.(4.1), the blurred image WNT = f(x, y) is not a sufficiently

close approximation to the true solution w(x, y, T) in Eqs.(3.1), (3.2). Clearly, that

nonlinear initial value problem is a difficult test case. All remaining examples in this

paper will involve blurring using Eq.(3.3), an equation with milder nonlinearities.

US Air Force resolution chart. Figure 8.3 involves a 512 × 512 image,

and is designed to show the possible merits of using the compensated backward in

time explicit scheme in Eq. (7.1), versus the more orthodox nonlinear Van Cittert

method in Eq. (5.3). The scheme in Eq. (4.1) was applied to Eq. (3.3), with

Δx = Δy = 1/512, Δt = 1.5 × 10−7, and 600 time steps, to yield the blurred image at the final time T = 9.0 × 10−5. That blurring computation was accomplished

in 3 minutes. As seen in the leftmost image in Figure 8.3, due to the pronounced

anisotropic character of q(x, y,t) in Eq. (3.3), the vertical sequence of numbers near

1Mention of commmercial products or services in this report does not imply NIST appproval

or endorsement of these products or services, nor does it imply that such products or services are

necessarily the best available for the intended purpose.
8 A. S. CARASSO

the right edge has almost been erased. However, after 100 Van Cittert iterations in

the rightmost image in Figure 8.3, that vertical sequence is partially restored, along

with noticeable improvements in other parts of the image. As may be expected, that

Van Cittert reconstruction required almost 5 hours, as it necessitated 100 sequential

nonlinear blurring computations, each lasting 3 minutes.

Using E = 3.5 × 10−4, p = 2.5, in the compensated backward explicit scheme in

Eq. (7.1), produced the middle image in Figure 8.3. The vertical sequence of numbers

near the right edge appears better defined in the middle image than in the rightmost

Van Cittert image. However, certain other parts of the image are marginally better

recovered in the Van Cittert restoration. The above choice of (E, p) was arrived at

after several interactive trials, each lasting 4 minutes. Clearly, on a parallel computer,

several such pairs (E, p) can be explored simultaneously. No parameter choice other

than λ in Eq. (5.3) is necessary in the Van Cittert method.

Sagittal Brain MRI image. The 256 × 256 MRI Brain image in Figure 8.4

was blurred by applying the scheme in Eq.(4.1) to Eq. (3.3), with Δx = Δy =

1/256, Δt = 3.0 × 10−7, and capturing the resulting image after 400 time steps at

the final time T = 1.2 × 10−4. That computation required 20 seconds. As in Figure

8.3, there is substantially more blurring near the right edge of the blurred image. A

competent restoration results after 100 Van Cittert iterations and 2000 seconds, as

shown in the rightmost image in Figure 8.4. However, after a few interactive trials,

each requiring 25 seconds, the compensated backward explicit scheme in Eq. (7.1)

with E = 2.25 × 10−5, p = 2.875, produced the almost as good reconstruction shown

in the middle image in Figure 8.4.

Ben Franklin image. Identical mesh parameters were used to produce the

256×256 blurred Ben Franklin image in Figure 8.5. Recovery of the right half of that

image is again challenging. Quite useful reconstruction was achieved after 100 Van

Cittert iterations and 2000 seconds, as shown in the rightmost image in Figure 8.5.

However, after a few 25 second trials, reasonably good recovery was achieved with the

compensated backward explicit scheme in Eq. (7.1), using E = 4.0 × 10−5, p = 2.775.

As in the preceding MRI Brain image, the Van Cittert image is again somewhat

sharper.

Satellite image. Behavior in the 256 × 256 Satellite image differs from the

preceding two cases, although identical mesh parameters were applied to generate

the blurred image in Figure 8.6, using Eq. (3.3). That blurred image appears more

degraded than the blurred Satellite image in Figure 8.2, where a more severe nonlinear

process was used. In contrast to the case in Figure 8.2, very useful reconstruction is

achieved after 100 Van Cittert iterations and 2000 seconds, as shown in the rightmost

image in Figure 8.6. Here, interestingly, after a few 25 second interactive trials, equally

good results can be achieved using the compensated backward explicit scheme in Eq.

(7.1) with E = 6.5 × 10−4, p = 2.35. Behavior in this example is similar to that found

in the previously discussed 512 × 512 US Air Force resolution chart.

9. Compensating operators B in general domains. One may consider backward in time continuation in nonlinear second order parabolic problems on general

domains Ω in R2, with homogeneous Dirichlet or Neumann boundary conditions.

Stabilizing backward in time marching computations by means of a compensating

operator B can be accomplished as follows. For the selfadjoint spatial differential
STABLE TIME-MARCHING IN ILL-POSED EVOLUTION EQUATIONS 9

 Resolution chart Brain MRI

Ben Franklin Satellite

SHARP TEST IMAGES IN NONLINEAR

 DEBLURRING EXPERIMENTS

Fig. 8.1. A collection of test images is useful in exploring the feasibility of backward in time

continuation in nonlinear parabolic equations, even when the intended application is unrelated to

imaging.
10 A. S. CARASSO

BLURRED IMAGE USING EQ. (3.1) AFTER 20 VAN CITTERT ITNS.

AFTER 20 VAN CITTERT ITNS.

BLURRED IMAGE USING EQ. (3.1) AFTER 100 VAN CITTERT ITNS.

BLURRED IMAGE USING EQ. (3.1)

FAILED BACKWARD IN TIME RECONSTRUCTIONS

IN NONLINEAR PARABOLIC EQUATION (3.1−3.2) .

Fig. 8.2. Failed image recovery in Eqs.(3.1), (3.2), may indicate adverse H ̈older-continuous

behavior in stability estimate Eq. (2.2) associated with Eqs. (3.1), (3.2), and highlights possibility

of untrustworthy backward continuation in other intended applications of Eqs. (3.1), (3.2).
STABLE TIME-MARCHING IN ILL-POSED EVOLUTION EQUATIONS 11

USEFUL DEBLURRING OF 512x512 NONLINEARLY BLURRED IMAGERY

Blur with anisotropic nonlinearparabolic Eq. (3.3). (3 min).

Deblur with stabilized backward intime marching (4 min).

Deblur with 100 Van Cittertiterations (5 hrs).

Fig. 8.3. In 512 × 512 resolution chart image, stabilized backward time-marching in Eq. (3.3), with = 3.5 × 10−4, p = 2.5, competently recovers vertical

sequence of numbers near image right edge in relatively little computing time, while the Van Cittert iteration is computationally impractical on large problems.
 

12 A. S. CARASSO

USEFUL DEBLURRING OF 256x256 NONLINEARLY BLURRED IMAGERY

Blur with anisotropic nonlinearparabolic Eq. (3.3). (20 sec).

Deblur with stabilized backward intime marching (25 sec). 

Deblur with 100 Van Cittertiterations (2000 sec).

Fig. 8.4. Stabilized backward time-marching in Eq. (3.3), with = 2.25 × 10−5, p = 2.875, produces useful reconstruction in sagittal MRI brain image.

However, computationally intensive Van Cittert iteration results in moderately sharper image.
 

STABLE TIME-MARCHING IN ILL-POSED EVOLUTION EQUATIONS 13

USEFUL DEBLURRING OF 256x256 NONLINEARLY BLURRED IMAGERY

Blur with anisotropic nonlinearparabolic Eq. (3.3). (20 sec).

Deblur with stabilized backward intime marching (25 sec). 

Deblur with 100 Van Cittertiterations (2000 sec).

Fig. 8.5. Stabilized backward time-marching in Eq. (3.3), with = 4.0 × 10−5, p = 2.775, competently restores Ben Franklin image. However, with

significantly greater computational effort, Van Cittert iteration produces moderately sharper image.
 

14 A. S. CARASSO

USEFUL DEBLURRING OF 256x256 NONLINEARLY BLURRED IMAGERY

Blur with anisotropic nonlinearparabolic Eq. (3.3). (20 sec).

Deblur with stabilized backward intime marching (25 sec). 

Deblur with 100 Van Cittertiterations (2000 sec).

Fig. 8.6. As was the case in the resolution chart image in Figure 8.3, stabilized backward time-marching in Eq. (3.3), with = 6.5 × 10−4, p = 2.35,

produces a reconstructed Satellite image comparable in visual quality to that resulting from the computationally intensive Van Cittert iteration.
STABLE TIME-MARCHING IN ILL-POSED EVOLUTION EQUATIONS 15

operator −∆ on Ω, with the same homogeneous boundary conditions as in the original parabolic equation, precompute and store a reasonably large number of eigenvalues {σm}M

m=1, with corresponding orthonormal eigenfunctions φm(x, y), where

0 ≤ σ1 ≤ σ2 · · · ≤ σM. Analogously to the situation in Eq. (6.5), given , p, let

∆t be the fixed time step in a step by step finite difference marching computation

on Ω. Let v(x, y) be the resulting numerical solution at some discrete time tk = k∆t

in the evolving computation. Expanding in the precomputed eigenfunctions of −∆,

we may then define the compensating linear operator B(, p, ∆t) applied to v(x, y) as

follows

Bv(x, y) = X

M

m=1

exp{−∆t(σm)

p

}

Ω

(9.1) v(x, y)φm(x, y)dxdy φm(x, y).

Z 

As in Eq. (7.1), this may be incorporated in a backward time-marching algorithm on

the spatial domain Ω. By using (−∆)p with p > 1, and choosing p interactively, it

may be possible to successfully march backward in time in difficult nonlinear parabolic

equations where an adjoint spatial differential operator is ill-defined, and for which

the quasi-reversibility formulation in Eq. (7.3) is inapplicable.

10. Concluding remarks. Due to the lack of reliable computational tools,

backward in time continuation in multidimensional nonlinear parabolic equations remains largely uncharted territory. Two approaches were presented in this paper and

applied to interesting nonlinear examples in 2D rectangular regions. The first method,

introduced in [8], is a global in time nonlinear Van Cittert iterative procedure. While

the linear convergence analysis is inapplicable in the nonlinear case, the Van Cittert iteration can produce useful reconstructions in a significant class of parabolic

equations. However, the computational cost becomes prohibitive for large problems.

The second method is based on stepwise time-marching schemes, the universal

computational approach in well-posed nonlinear evolution equations. Unfortunately,

when the evolution equation is ill-posed, such schemes are necessarily computationally unstable and lead to explosive noise amplification. Here, an easily synthesized

compensating operator, based on (−∆)p with p > 1, is applied at each time step to

quench the instability. In 2D rectangular regions, such compensating operators can

be simulated using FFT algorithms. In more general domains Ω, compensating operators can be simulated using a sufficiently large number of precomputed eigenpairs

(σm, φm) for −∆ on Ω, with the same homogeneous boundary conditions on ∂Ω as

in the original nonlinear parabolic equation. An easily implemented compensated explicit scheme, marching backward in time, was constructed in Section 7 and applied

in Section 8.

A useful idea in exploring possible backward in time continuation in a given

nonlinear parabolic equation, is to associate a fictitious image deblurring problem

with that equation. Numerically blurring a sharp image by marching forward in

time up to some fixed time T > 0, creates interesting noisy data that can be used

to attempt reconstruction of the sharp image by backward continuation to t = 0.

Such preliminary exploration may be instructive even if the intended engineering

application is unrelated to imaging. The failed reconstructions in Figure 8.2 are useful

examples that signal possible untrustworthy backward continuation in Eqs.(3.1), (3.2).

The experiments in Figures 8.3 through 8.6 indicate that the compensated backward time-marching explicit scheme in Section 7 can produce reconstructions that are

often comparable in visual quality to those obtained using the Van Cittert method,
16 A. S. CARASSO

but at greatly reduced computational cost. A major advantage of the scheme in Section 7 is the ability to explore various non integer trial values of p within the same

program. Such versatility is not present in the quasi-reversibility method [14], and, in

fact, the QR formulation in Eq. (7.3) is not readily applicable to nonlinear problems.

Significantly, in the above nonlinear reconstruction experiments, simple compensating operators based on (−Δ)p with 2 < p < 3, were found adequate in controlling

instability.

REFERENCES

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

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

[2]	 A. C. Bagtzoglou and J. Atmadja, Mathematical methods for hydrologic inversion: The

case of pollution source identification, Handb Environ Chem Vol. 5, Part F,Vol. 3 (2005):

pp. 6596. DOI 10.1007/b11442.

[3]	 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

[4]	 W. Freeden and H. Nutz, Satellite gravity radiometry as tensorial inverse problem, Int J

Geomath (2011), pp. 177-218. DOI 10.1007/s13137-011-0026-x.

[5]	 W. Freeden, Z. Nashed, and T. Sonar (editors), Handbook of Geomathematics, 2nd edition,

Springer, Heidelberg, 2013.

[6]	 R. J. Blakely, Potential Theory in Gravity and Magnetic Applications, Cambridge University

Press, Newv York, 1995.

[7]	 A. S. Carasso, Bochner subordination, logarithmic diffusion equations, and blind deconvolution of Hubble space telescope imagery and other scientific data, SIAM J. Imaging Sciences,

3 (2010), pp. 954–980.

[8]	 A. S. Carasso, Hazardous continuation backward in time in nonlinear parabolic equations,

and an experiment in deblurring nonlinearly blurred imagery, J. Res. NIST, 118 (2013),

http:/dx.doi.org/10.6028/jres.118.010.

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

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

[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]	 P. A. Jansson, Deconvolution with Application in Spectroscopy, Academic Press, New York,

1984.

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

Academic Publishers, Boston, 1991.

[13]	 A. S. Carasso, A stable marching scheme for an ill-posed initial value problem, in Improperly Posed Problems and their Numerical Treatment, (Oberwolfach Conference, September

1982). G. H ̈ammerlin and K. H. Hoffman, eds. Birkh ̈auser-Verlag, Boston, 1983.

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

1969.