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 nonlinearcoefficient in a first order PDE via final observation. The problem isstated as an optimal control problem and solved numerically. Implicitfinite difference scheme is used for the approximation of the stateequation. A space of control variables is approximated by a sequence offinite-dimensional spaces with increaing dimensions. Finite dimensionalproblems are solved by a gradient method and numerical results arepresented.
1. Introduction
In this paper we consider the following nonlinear initial boundary-value problem
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,dadc≥0 for a.a. c∈[0,1]}.
(2)
Let us define the cost functional for the control optimization problem by
J(c)=12∫01(a(c(x,T))−φ(x))2dx.
(3)
The problem under consideration can be posed as the following optimal control
problem:
find a(c)∈Asuch thatit 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=1Nuuiψ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
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:
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=1Nuuiψ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 ω+=ω̄∖{x0}. 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̄=1h(v(x,t)−v(x−h,t)),vx=1h(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,ci)−φ(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)≤1for all x,t.
Proof.For a fixed time level
t
problem (8) can be written as
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)andc(v)≡c(t,v)be the solutions of (11),corresponding to u,v∈K.Thenthere exists a constant M=M(t,u)such that
∣∣c(t,u)−c(t,v)∣∣∞≤M∣∣u−v∣∣1,
(12)
where ∣∣.∣∣∞and∣∣.∣∣1are maximumand L1-normof vectors, respectively.
Proof.Let I
be unit matrix and σ=σ(u)=∑j=1NuujNu.
From equation (11) we obtain
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
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
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,ch(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” forNu=32without multilevel approach. The calculated results were similar to the previousones, but the number of iterations and especially the time of calculation weremuch greater than for the multilevel method.
Figure 1:
Objective function
φ1
and calculated
a(c(x,T))
(top) and computed nonlinear coefficient
a(c)
(bottom) for
h=τ=0.01
Figure 2:
Objective function
φ1
and calculated
a(c(x,T))
(top) and computed nonlinear coefficient
a(c)
(bottom) for
h=τ=0.005
Figure 3:
Objective function
φ2
and calculated
a(c(x,T))
(top) and computed nonlinear coefficient
a(c)
(bottom) for
h=τ=0.01
Figure 4:
Objective function
φ2
and calculated
a(c(x,T))
(top) and computed nonlinear coefficient
a(c)
(bottom) for
h=τ=0.005
Figure 5:
Objective function
φ3
and calculated
a(c(x,T))
(top) and computed nonlinear coefficient
a(c)
(bottom) for
h=τ=0.01
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 ÉquationParabolique 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, SIAMJ. Math. Analysis, 1997, Vol.28, No.3, pp.611-632
[3]Igler B., Totsche K.U., Knabner P. Unbiased Identification of Nonlinear SorptionCharacteristics by Soil Column Breakthrough Experiments, Preprint IAM Erlangen,1997, No. 224
[4]Knabner P., Igler B. Structural Identification of Nonlinear Coefficient Functions inTransport 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