A. Lapin and S. Lapin

Lobachevskii Journal of Mathematics Vol. 14, 2004, 69–84

©Alexander Lapin and Serguei Lapin

Alexander Lapin and Serguei Lapin
IDENTIFICATION OF NONLINEAR COEFFICIENT IN A TRANSPORT EQUATION

________________

2000 Mathematical Subject Classification. 65M32, 65N21.

Key words and phrases. nonlinear coefficient identification, transport equation, final observation, finite difference scheme, multilevel algorithm.

This work was supported by RFBR, grant 01-01-00070.

ABSTRACT. Considered a problem of identification a nonlinear coefficient in a first order PDE via final observation. The problem is stated as an optimal control problem and solved numerically. Implicit finite difference scheme is used for the approximation of the state equation. A space of control variables is approximated by a sequence of finite-dimensional spaces with increaing dimensions. Finite dimensional problems are solved by a gradient method and numerical results are presented.

1. Introduction

In this paper we consider the following nonlinear initial boundary-value problem

∂c ∂t + ∂c ∂x + ∂a(c) ∂t = 0,x ∈ (0, 1),0 < t ≤ T, c(0,t) = 1, c(x, 0) = 0 (1)

which models a convective transport of sorption chemical through a porous medium. Here c is the dissolved concentration of a chemical and a(c) is a so-called sorption isotherm. Function a(c) is the unknown of the problem, so we consider a structure identification problem. For physical reasons we assume a(c) to be continuous and non-decreasing and a(0) = 0. Under these assumptions there exists a unique solution c(x,t) to the problem (1) which takes its values from segment [0, 1]. In order, to define a(c) on [0, 1], we use a final observation φ(x): we try to choose a(c) in such a way that a(c(x,T)) = φ(x),(x ∈ [0, 1]), and where φ(x) is a non-negative continuous function. It is obvious that the formulated inverse problem is ill-posed, because it is not solvable for the arbitrary function φ(x).

We set up the above problem as an optimal control problem [5]. This approach is well-known in parameter identification problems (e.g. [1], [2]) and it is often called the output least squares method (cf. [3], [4]). In this paper we will concentrate on the numerical solution to the problem rather than its theoretical study.

In order to solve the optimal control problem we approximate the set of admissible coefficients a(c) by a finite dimensional set. We also approximate problem (1) by a finite difference scheme. The existence of an unique solution to the finite-dimensional optimal control problem is shown. We use a gradient-type method for its solution, where the gradient information is calculated via the solution of an adjoint state problem. Due to the highly ill-posedness of the problem we chose an approach which is characterized by increasing the dimension of the set for admissible coefficient (cf.[4]). The results of the numerical experiments are presented.

2. Formulation of the problem and its discretization

We consider problem (1) with non-linear ”coefficient” a(c) which belongs to the following subset of the Lipshitz-continuous on [0, 1] functions:

A = {a(c) ∈ C(0,1)[0, 1] : a(0) = 0, da dc ≥ 0 for a.a. c ∈ [0, 1]}. (2)

Let us define the cost functional for the control optimization problem by

J(c) = 1 2 ∫ 01(a(c(x,T)) − φ(x))2dx. (3)

The problem under consideration can be posed as the following optimal control problem:

find a(c) ∈ A such that it minimizes J(c) when c(x,t) satisfies the initial boundary-value problem (1).

Now in order to solve and stabilize the above stated optimal control problem we approximate the set A by a finite-dimensional set Ah. Ah is constructed by discretization of the coefficient space. Namely, let ah (c) ≡ a(u,c) = ∑ i=1Nu uiψi(c), where the functions ψi (c) compose a basis of a finite-dimensional space Uh, containing Ah, while u = {u1,…,uNu}T belongs to a set K of admissible parameters: u ∈K⇔ ah(c) ∈ Ah. After an approximation of set A we derive the problem of minimization to the functional

I(u) ≡ J(u,c) = 1 2 ∫ 01(a(u,c(x,T)) − φ(x))2dx (4)

with c = c(u) satisfying (1) and u ∈K.

For the differentiable function a(u,c) ∈ Ah the functional I(u) is also differentiable and

∇I(u) = ∫ 01(a(u,c(x,T)) − φ(x))a u′(u,c(x,T))dx + ∫ 0T ∫ 01 ∂ ∂tau′(u,c(x,t))λ(x,t)dxdt. (5)

Here c(x,t) is the solution of the state problem (1) and λ(x,t) is the solution of the corresponding adjoint problem:

∫ 0T ((1 + a c′(u,c(x,t))∂λ ∂t + ∂λ ∂x)μ(x,t)dxdt = ∫ 01(1 + a c′(x,T)λ(x,T)μ(x,T))dx + ∫ 0T λ(1,t)μ(1,t)dt + ∫ 01(a(u,c(x,T)) − φ(x))a c′(u,c(x,T)μ(x,T)dx, ∀μ(x,t) : μ(0,T) ≡ μ(x, 0) ≡ 0. (6)

We can rewrite (6) in pointwise form:

(1 + ac′(u,c))∂λ ∂t + ∂λ ∂x = 0,0 < x < 1,0 < t ≤ T, λ = 0,x = 1, 0 < t ≤ T (1 + ac′(u,c))λ + (a(u,c) − φ(x))a c′(u,c) = 0,0 < x < 1,t = T.

In order to solve numerically the optimal control problem we discretize the state problem (1) and construct corresponding discrete adjoint state problem. Let us divide the segment [0, 1] by the uniform mesh with stepsize hu = 1∕Nu, Nu ∈N, and define finite-dimensional spaces

Uh = {a(c) : a(c) is piecewise linear and continuous on [0, 1],a(0) = 0}

and

Ah = {a(c) ∈ Uh : a(ihu) ≤ a((i + 1)hu) for all i = 0, 1,...,Nu}.

Now we consider set {ψi(c)}i=1Nu which is the basis of Uh. This set consists of the piecewise polynomial functions which are defined the following way:

ψi(c) = 0 for 0 ≤ c ≤ (i − 1)hu Nu(c − (i − 1)hu)for(i − 1)hu < c ≤ ihu 1 for ihu < c ≤ 1 (7)

Use of the basis (7) instead of usual Courant basis allows us to obtain the simplest form of the set Ah via the nodal parameters of the functions a(c) ∈ Ah:

a(c) ≡ a(u,c) = ∑ i=1Nu uiψi(c) iff u ∈K≡ (RNu )+.

Now, let ω̄ = {xi = ih,i = 0, 1,..,N} be an uniform mesh on [0, 1] with mesh step-size h and ω+ = ω̄ ∖{x 0}. We also introduce ω̄τ = {tk = kτ,k = 0, 1,..,Nτ} being uniform mesh on [0,T] with step-size τ, ωτ + = ω̄ τ ∖{t0}, ωτ − = ω̄ τ ∖{t = T}.

We define the finite differences in time as:

vt̄ = 1 τ(v(x,t) − v(x,t − τ)),vt = 1 τ(v(x,t + τ) − v(x,t))

and in space as:

vx̄ = 1 h(v(x,t) − v(x − h,t)),vx = 1 h(v(x + h,t) − v(x,t)).

State equation (1) is approximated by the following implicit scheme:

(ch + ah)t̄(x,t) + chx̄(x,t) = 0 for x ∈ ω+,t ∈ ω τ+, ch(0,t) = 1, for t ∈ω̄τ, ch(x, 0) = 0 for x ∈ ω+, (8)

where ah = a(u,ch). Applying the quadrature formula to (4), we derive the following finite-dimensional cost functional:

I(u) ≡ Jh(u,ch) = ∑ i=1N(a(u,c i) − φ(xi))2, (9)

where ci = ch(xi,T).

Futher we consider the optimal control problem (OCP):

find minJh(u,ch) for u ∈K, when equation (8) is fulfilled. (10)

Lemma 1. Discretization scheme (8) has a unique solution ch(x,t) for any u ∈K, and 0 ≤ ch(x,t) ≤ 1 for all x,t.

Proof. For a fixed time level t problem (8) can be written as

(1 + τ h)ch(x,t) + a(u,ch(x,t)) = ch(x,t−τ) + a(u,ch(x,t−τ)) + τ hch(x−h,t).

Because ch(0,t) = 1, we solve recurrently the equations with monotone increasing and continuous functions to find ch(x,t) for all x ≥ h, whence the unique solvability follows. Further, owing to the non-negativeness of the initial and boundary conditions we obviously have ch(x,t) ≥ 0. Now, let ch(x,t − τ) ≤ 1 and ch(x,t) ≤ 1 for x ≤ xi. If we suppose that ch(xi+1,t) > 1, then from (8) it follows (ch + ah)t̄(xi+1,t) < 0, thus, ch(xi+1,t) < ch(xi,t) and we get a contradiction. □

In addition, we prove that ch is a Lipschitz-continuous function of u. To do this we rewrite equation (8) for a fixed time level in an algebraic form. Let c = c(t) = c(t,u) ∈RN be the vector of nodal values to ch(x,t) : ci = ch(xi,t) for i = 1,…,N, while f = f(u) ∈RN be the vector with coordinates fi = ch(xi,t − τ) + a(u,ch(xi,t − τ)) for i = 2,…,N and f1 = ch(x1,t − τ) + a(u,ch(x1,t − τ)) + τ∕h. Let further B ∈RN×N be the two-diagonal matrix with diagonal elements 1 + τ∕h and off-diagonal ones −τ∕h. Then (8) for a fixed time level t has the form

Bc + a(u,c) = f(u). (11)

Lemma 2. Let c(u) ≡ c(t,u) and c(v) ≡ c(t,v) be the solutions of (11), corresponding to u,v ∈K. Then there exists a constant M = M(t,u) such that

∣∣c(t,u) − c(t,v)∣∣∞ ≤ M∣∣u − v∣∣1, (12)

where ∣∣.∣∣∞ and ∣∣.∣∣1 are maximum and L1-norm of vectors, respectively.

Proof. Let I be unit matrix and σ = σ(u) = ∑ j=1Nu ujNu. From equation (11) we obtain

c(u) = (B + σI)−1(σc(u) − a(u,c(u)) + f(u)),

c(v) = (B + σI)−1(σc(v) − a(v,c(v)) + f(v))

and estimate the difference

c(u) − c(v) = (B + σI)−1(σ(c(u) − c(v)) − a(u,c(u)) + a(u,c(v)) + (B + σI)−1(f(u) − f(v) + (a(v,c(v)) − a(u,c(v))) = (B + σI)−1(∑ j=1Nu uj(Nu(c(u) − c(v)) + ψj(c(v)) − ψj(c(u))) + (B + σI)−1(f(u) − f(v) + ∑ j=1Nu (vj − uj)ψj(c(v))).(13)

Let Lj be the diagonal matrix with entries, defined the following way: if ci (u) ⁄= ci(v) then the i-th diagonal entry of Lj is ψj (ci(u)) − ψj(ci(v)) ci(u) − ci(v) ; otherwise it is equal to 0. With this definition one has

ψj(c(u)) − ψj(c(v)) = Lj(c(u) − c(v)). (14)

Owing to the definition of the functions ψj

0 ≪ Lj ≪ NuI, (15)

where by ≪ the componentwise inequality is denoted. From (14) and (15) we obtain

∣∑ j=1Nu uj(Nu(c(u) − c(v)) + ψj(c(v)) − ψj(c(u))∣≪

∑ j=1Nu uj∣(Lj − NuI)(c(u) − c(v))∣≪ σ∣(c(u) − c(v))∣. (16)

Direct calculations show that

∣∣(B + σI)−1∣∣ ∞ ≤ 1 1 + σ. (17)

Now, (13), (16) and (17) yield

∣∣c(u)−c(v)∣∣∞ ≤ (1+σ)∣∣(B+σI)−1(f(u)−f(v)+∑ j=1Nu (vj−uj)ψj(c(v)))∣∣∞

≤∣∣f(u) − f(v)∣∣∞ + ∑ j=1Nu ∣vj − uj∣.

Because f(u) − f(v) = c(t − τ,u) − c(t − τ,v) + a(u,c(t − τ,u)) − a(v,c(t − τ,v)) and similarly to the previous estimates we have

∣a(u,c(t − τ,u)) − a(v,c(t − τ,v))∣≪ σ∣c(t − τ,u) − c(t − τ,v)∣ + ∣∣u − v∣∣1,

then

∣∣f(u) − f(v)∣∣∞ ≤ (1 + σ)∣∣c(t − τ,u) − c(t − τ,v)∣∣∞ + ∣∣u − v∣∣1.

Thus,

∣∣c(t,u) − c(t,v)∣∣∞ ≤ (1 + σ)∣∣c(t − τ,u) − c(t − τ,v)∣∣∞ + 2∣∣u − v∣∣1,

whence inequality (12) follows with M(t,u) = 2(1 + σ(u))t∕τ σ(u) . □

Theorem 1. Problem (10) has a solution for any u ∈K

Proof. Owing to the previous Lemma 2 and definition of a(u,c), the cost function I(u) is continuous, and obviously it is coercive: ∣∣u∣∣→∞⇒ I(u) → +∞, whence the result. □

3. Iterative method

A function a(u,ch) ∈ Ah is not differentiable in ch. However, it is Lipcshitz-continuous and has left and right derivatives ac ′ . Below we use the piecewise constant function ac′ (by taking either left or right derivative) to construct adojnt state and to receive ”gradient” information. To construct the adjoint state problem we define the following Lagrange function

Lh(u,ch,λh) = Jh(u,ch)+τ ∑ t∈ωτ+h∑ x∈ω+λh(x,t)((ch+ah)t̄(x,t)+chx̄(x,t)), (18)

where the mesh function λh(x,t) vanishes when x = 0 or t = 0. Stationary points of Lagrange function Lh are defined from the following system

∇λLh ⋅ δλh ≡ τ ∑ t∈ωτ+h∑ x∈ω+((ch + ah)t̄(x,t) + chx̄(x,t))δλh(x,t) = 0; (19)

∇cLh ⋅ δch ≡ h∑ x∈ω+(a(u,ch(x,T)) − φh(x))ac′(u,c h(x,T))δch(x,T)

+τ ∑ t∈ωτ+h∑ x∈ω+((δch + ac′δc h)t̄(x,t) + (δch)x̄(x,t))λh(x,t) = 0; (20)

∇uLh ⋅ δu ≡ h∑ x∈ω+(a(u,ch(x,T)) − φh(x))au′(u,c h(x,T))δu

+τ ∑ t∈ωτ+h∑ x∈ω+λh(x,t)(au′) t̄(x,t)δu = 0, (21)

where δu ∈Rm and trial mesh functions δλh and δch vanish when x = 0 or t = 0.

Equation (20) gives us the adjoint state problem which has the following pointwise form

−(1 + ac′)λ t − λx = 0, 0 < x < 1,t ∈ ωτ−, −(1 + ac′)λ t + 1 hλx = 0, x = 1, t ∈ ωτ−, (1 + ac′)λ − τλ x + ac′(a − φ) = 0,0 < x < 1, t = T, (1 + ac′)λ + τ hλx + ac′(a − φ) = 0, x = 1, t = T. (22)

From the equation (21) we derive the formula for calculation the gradient ∇I(u) = ∇uLh :

∇I(u) = h∑ x∈ω+(a(u,ch(x,T)) − φh(x))au′(u,c h(x,T))

+τ ∑ t∈ωτ+h∑ x∈ω+λh(x,t)(au′) t̄(x,t) (23)

Now to minimize the functional I(u) we use the gradient method:

uk+1 = uk − ρ k+1∇I(uk), (24)

where an initial guess u0 ∈K, and iterative parameter ρk+1 is defined via the line search technique. To improve the convergence of the method we use “multilevel” implementation; namely, we first solve the problem with the dimension Nu = 2 of the space Uh and then increase Nu. This approach gives us the possibility to achieve better results than if using fixed dimension Nu.

4. Numerical results

In this section, we will describe numerical experiments that we have performed for the solution of the above defined optimal control problem.

We discretize domain Ω = (0,L) × T by uniform mesh with step size h and time step Δt. Then, we define discrete control space Uh with the dimension Nu. For the numerical experiments we have used the following objective functions φ:

φ1(x) = 0.5(1.0 − x),

φ2(x) = 1,0 ≤ x ≤ L∕2, 0,L∕2 < x ≤ L,

φ3(x) = 0.5 − x,0 ≤ x ≤ L∕2, 0, L∕2 < x ≤ L.

We have performed calculations for variable Nu = 2, 4, 8, 16, 32, where on each step after increasing Nu we used the previous solution as an initial data.

Figures 1 - 6 show evaluated coefficient a(c) and comparison between objective function φ and computed values of a(c(x,T)) for three considered functions φ and for different mesh and time step sizes.

Remark 1. We also performed calculations ”directly” for Nu = 32 without multilevel approach. The calculated results were similar to the previous ones, but the number of iterations and especially the time of calculation were much greater than for the multilevel method.

PIC PIC

Figure 1: Objective function φ1 and calculated a(c(x,T)) (top) and computed nonlinear coefficient a(c) (bottom) for h = τ = 0.01

PIC PIC

Figure 2: Objective function φ1 and calculated a(c(x,T)) (top) and computed nonlinear coefficient a(c) (bottom) for h = τ = 0.005

PIC PIC

Figure 3: Objective function φ2 and calculated a(c(x,T)) (top) and computed nonlinear coefficient a(c) (bottom) for h = τ = 0.01

PIC PIC

Figure 4: Objective function φ2 and calculated a(c(x,T)) (top) and computed nonlinear coefficient a(c) (bottom) for h = τ = 0.005

PIC PIC

Figure 5: Objective function φ3 and calculated a(c(x,T)) (top) and computed nonlinear coefficient a(c) (bottom) for h = τ = 0.01

PIC PIC

Figure 6: Objective function φ3 and calculated a(c(x,T)) (top) and computed nonlinear coefficient a(c) (bottom) for h = τ = 0.005

References

[1]   Chavent G., Lemonnier P., Identification de la Non-Linearité D’Une Équation Parabolique Quasilineaire Applied Mathematics and Optimization, Springer Verlag, 1974, vol.1, No.2,pp.121-162

[2]   DuChateau P. An Inverse Problem for the Hydraulic Properties of Porous Media, SIAM J. Math. Analysis, 1997, Vol.28, No.3, pp.611-632

[3]   Igler B., Totsche K.U., Knabner P. Unbiased Identification of Nonlinear Sorption Characteristics by Soil Column Breakthrough Experiments, Preprint IAM Erlangen, 1997, No. 224

[4]   Knabner P., Igler B. Structural Identification of Nonlinear Coefficient Functions in Transport Processes through porous media, in Lectures on Applied Mathematics, Springer Verlag, 2000, pp.157-178

[5]   Lions J.L.Optimal Control of Systems Governed by Partial Differential Equations, Springer Verlag, New York, 1971

KAZAN STATE UNIVERSITY, KAZAN, 420008, RUSSIA

E-mail address: alapin@ksu.ru

UNIVERSITY OF HOUSTON, HOUSTON, TX, 77204, USA

E-mail address: slapin@math.uh.edu