# Identification of Flexural Rigidity in a Kirchhoff Plates Model Using a Convex Objective and Continuous Newton Method

Authors: B. Jadamba; R. Kahler; A. A. Khan; F. Raciti; B. Winkler
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

This work provides a detailed theoretical and numerical study of the inverse problem of identifying flexural rigidity in Kirchhoff plate models. From a mathematical standpoint, this inverse problem requires estimating a variable coefficient in a fourth-order boundary value problem. This inverse problem and related estimation problems associated with general plates and shell models have been investigated by numerous researchers through an optimization framework using the output least-squares (OLSs) formulation. OLS yields a nonconvex framework and hence it is suitable for investigating only the local behavior of the solution. In this work, we propose a new convex framework for the inverse problem of identifying a variable parameter in a fourth-order inverse problem. Existence results, optimality conditions, and discretization issues are discussed in detail. The discrete inverse problem is solved by using a continuous Newton method. Numerical results show the feasibility of the proposed framework.


---

## Body

## 1. Introduction

LetΩ be a nonempty bounded domain in R
2 with a sufficiently smooth boundary Γ. We consider the following fourth-order elliptic boundary value problem (BVP):(1a)
Δ
a
Δ
u
=
f
in
Ω
,
(1b)
u
=
0
on
Γ
,
(1c)
∂
u
∂
n
=
0
on
Γ
,where ∂
u
/
∂
n is the usual normal derivative.The above BVP models the pure bending (without twisting) of a Kirchhoff plate occupying the regionΩ. Here, the parameter a
(
x
,
y
) is the flexural rigidity and the solution u is the lateral deflection under the force f per unit area (see [1] for details). The parameter a is connected to the thickness of the plate, Young’s modulus of elasticity, and the Poisson ratio of the material. Boundary conditions (1b) and (1c) are the so-called clamped boundary conditions. However, all the results given below can easily be transferred to the pinned boundary conditions given by(2)
u
=
Δ
u
=
0
on
Γ
.Moreover, the proposed framework can be extended to the following more general plate model [2]:(3)
μ
Δ
a
Δ
u
+
1
-
μ
a
u
x
x
x
x
+
2
a
u
x
y
x
y
+
a
u
y
y
y
y
=
f
in
Ω
,
u
=
∂
u
∂
n
=
0
on
Γ
.Given the forcef and the flexural rigidity a, the direct problem is to find the lateral deflection u. In this work, our objective is to identify the flexural rigidity a
(
x
,
y
) so that the corresponding deflection u is the nearest to a deflection measurement z in a certain sense. This inverse problem and its dynamic analogue have been explored mostly through the output least-squares formulation. In the present context, the output least-squares (OLSs) approach seeks to minimize the functional:(4)
a
⟶
1
2
u
a
-
z
2
,defined by an appropriate norm. Here, z is the data (a measurement of u) and u
(
a
) is the unique solution of the weak form of (1a), (1b), and (1c) that corresponds to the flexural rigidity a
(
x
,
y
).One of the major difficulties associated with the OLS approach is the fact that the coefficient-to-solution mapu
:
a
→
u
(
a
) is nonlinear and hence, in general, the OLS functional is nonconvex. Consequently, the OLS-based approach only allows the study of local properties of the inverse problem. Moreover, the numerical algorithms designed through the necessary optimality conditions can only ensure convergence to a local optimizer and such conditions are not sufficient optimality conditions.In the following, we will give a quick overview of some related works. In an interesting paper, Lesnic et al. [3] investigate the coefficient identification problem in a fourth-order differential equation that governs the deflection u of the Euler-Bernoulli beam and give useful stability estimates. An optimization-based approach for the same problem was recently given in [4]. In a recent paper, Marinov and Marinova [5] developed a numerical method by using the so-called method of variational embedding for solving an inverse problem for bending stiffness estimation in a Kirchhoff-Love plate from overdetermined data. Ewing et al. [6] study a nonlinear analogue of this problem using the mixed finite element method.Manservisi and Gunzburger [7] studied a control problem associated with a plate model and give interesting applications to a simplified model for the “sag bending process” in the manufacturing of automobile windscreens. Inspired by this inverse problem and its useful application, Engl and Kügler [8] and Kügler [9] presented a regularized Landweber method for its numerical solution. Salazar and Westbrook [10] studied the identification of the flexural rigidity of a nonhomogeneous plate of uniform thickness under the assumption that Poisson’s ratio is constant and only the Young modulus varies pointwise. Lopes et al. [11] used identification techniques to detect holes in plates. In a related work, Kim and Kreider [12] conducted a numerical study of the inverse problem of parameter identification in nonlinear elastic and viscoelastic plates. They presented numerical results for a variety of constitutive models. In an intriguing paper by Yuan and Yamamoto [13], they studied Lipschitz stability estimates for both inverse problems of determining spatially varying Lamé coefficients and determining the mass density by a finite number of boundary observations. This inverse problem for the case of interior measurement was studied in [14–16]. A systematic study of inverse problems of parameter identification in plates and shells models was conducted in a series of excellent papers by White [2, 17, 18] where the author investigates various aspects of the inverse problem. An overview of various approaches for parameter identification problems in inverse problems can be found in [19].In this work, our primary objective is to employ a new convex functional for the identification of flexural rigiditya
(
x
,
y
). The convexity of the new functional circumvents one of the major drawbacks, nonconvexity, of the OLS functional. Furthermore, the proposed convex objective functional characterizes global optimality, and the associated variational inequality is a necessary as well as a sufficient optimality condition. The new objective functional also has a computational advantage in that to compute its gradient one does not need to compute the derivative of the coefficient-to-solution map.This paper is divided into six sections. Section2 begins with the essential background material and continues with the introduction of the related minimization problem which is subsequently shown to have a solution. The problem is discretized using finite elements in Section 3 and it is shown that the continuous minimization problem can be approximated by the discrete one. A description of the continuous Newton method is given in Section 4 and the details of numerical experiments and their results are provided in Section 5. The paper concludes with some remarks concerning the outlined approach and other related issues.


## 2. Modified Output Least-Squares (MOLSs)

The variational formulation of (1a), (1b), and (1c) plays an important role in formulating the MOLS approach. Assuming that the feasible parameters a
(
x
,
y
) are in L
∞
(
Ω
), f
∈
L
2
(
Ω
), and by taking V
=
H
0
2
(
Ω
), the weak form of (1a), (1b), and (1c) reads: find u
∈
V such that(5)
∫
Ω
a
Δ
u
Δ
v
=
∫
Ω
f
v
,
∀
v
∈
V
.In the following, we takeK
≔
{
a
∈
L
∞
∣
0
<
c
≤
a
(
x
,
y
)
} to be the set of admissible coefficients. A variety of other conditions can be imposed on the feasible coefficients a
(
x
,
y
) so that the following two inequalities hold for some constants α
>
0 and β
>
0:(6)
∫
Ω
a
Δ
u
Δ
v
≤
β
a
∞
u
V
v
V
,
∀
u
,
v
∈
V
,
∀
a
∈
K
,
(7)
α
u
V
2
≤
∫
Ω
a
Δ
u
Δ
u
,
∀
u
∈
V
,
∀
a
∈
K
.Motivated by [20, 21], we propose the following objective functional to identify the flexural rigidity:(8)
J
a
=
1
2
∫
Ω
a
Δ
u
a
-
z
Δ
u
a
-
z
,where z is the measured data and u
(
a
) solves the weak form (5).To give optimality conditions for the above minimization problem using (8) and to design some derivative-based numerical schemes, we need to compute the derivatives of J
(
·
). Since J depends on u
(
a
), the derivatives of the solution map are instrumental to evaluate the derivatives of the functional J. For this, we give the following result.Lemma 1.
For eacha in the interior of K, u is infinitely differentiable at a
. Given u
=
u
(
a
), the first derivative δ
u
=
D
u
(
a
)
δ
a is the unique solution of the variational equation:(9)
∫
Ω
a
Δ
δ
u
Δ
v
=
-
∫
Ω
δ
a
Δ
u
Δ
v
,
f
o
r
e
v
e
r
y
v
∈
V
,and the second derivative δ
2
u
=
D
2
u
(
a
)
(
δ
a
1
,
δ
a
2
) is the unique solution of the variational equation:(10)
∫
Ω
a
Δ
δ
2
u
Δ
v
=
-
∫
Ω
δ
a
2
Δ
D
u
a
δ
a
1
Δ
v
-
∫
Ω
δ
a
1
Δ
D
u
a
δ
a
2
Δ
v
,
f
o
r
e
v
e
r
y
v
∈
V
.Proof.
A direct proof relying on (6) and (7) can be extracted from [20]. A proof based on the implicit function theorem can be found in [2].The following result gives derivative formulas for (8).Theorem 2.
Leta be an arbitrary element in the interior of K. Then, consider the following:(1)
The first derivative ofJ is given by(11)
D
J
a
δ
a
=
-
1
2
∫
Ω
δ
a
Δ
u
a
+
z
Δ
u
a
-
z
.
(2)
The second derivative ofJ is given by(12)
D
2
J
a
δ
a
,
δ
a
=
∫
Ω
a
Δ
D
u
a
δ
a
Δ
D
u
a
δ
a
.Proof.
The first derivative ofJ is a direct consequence of the chain rule:(13)
D
J
a
δ
a
=
1
2
∫
Ω
δ
a
Δ
u
a
-
z
Δ
u
a
-
z
+
∫
Ω
δ
a
Δ
D
u
a
δ
a
Δ
u
a
-
z
.
We note that, from (9), we have(14)
∫
Ω
a
Δ
D
u
a
δ
a
Δ
u
a
-
z
=
-
∫
Ω
δ
a
Δ
u
a
Δ
u
a
-
z
,and hence(15)
D
J
a
δ
a
=
1
2
∫
Ω
δ
a
Δ
u
a
-
z
Δ
u
a
-
z
-
∫
Ω
δ
a
Δ
u
a
Δ
u
a
-
z
=
-
1
2
∫
Ω
δ
a
Δ
u
a
+
z
Δ
u
a
-
z
.
It now follows that(16)
D
2
J
a
δ
a
,
δ
a
=
-
1
2
∫
Ω
δ
a
Δ
D
u
a
δ
a
Δ
u
a
-
z
-
1
2
∫
Ω
δ
a
Δ
u
a
+
z
Δ
D
u
a
δ
a
=
-
∫
Ω
δ
a
Δ
D
u
a
δ
a
Δ
u
a
=
∫
Ω
a
Δ
D
u
a
δ
a
Δ
D
u
a
δ
a
,where in the last step we applied (9) again.The following result gives a useful feature of MOLS.Theorem 3.
The modified output least-squares functionalJ defined in (8) is convex in the interior of the set K
.Proof.
In view of the form of the second derivative ofJ, it follows from (7) that the following inequality holds for all a in the interior of K:(17)
D
2
J
a
δ
a
,
δ
a
≥
α
D
u
a
δ
a
2
,which proves the convexity of the functional J.The inverse problem of identifying parameters is ill-posed and hence we need to regularize the MOLS functional. In the following, we assume thatK is the set of admissible coefficients which we further assume to be closed and convex subset of B
≔
H
2
(
Ω
)
. We also assume that K is in the interior of the domain of the parameter-to-solution map so that we can use the convexity of MOLS. We now consider the following regularized minimization problem: find a
∗
∈
K by solving(18)
min
a
∈
K
⁡
J
ɛ
a
≔
1
2
∫
Ω
a
Δ
u
a
-
z
Δ
u
a
-
z
+
ɛ
2
R
a
,where ɛ
>
0 is the regularization parameter, z
∈
V is the data, and R
(
a
) is a regularization map. Throughout this work, we take R
(
a
)
≔
a
H
2
(
Ω
)
2. The associated inner product is denoted by ·
,
·
.We have the following useful information concerning a minimizer of the MOLS functional.Theorem 4.
The regularized minimization problem (18) is uniquely solvable. Moreover, a
∈
K is a minimizer if and only if the following variational inequality holds:(19)
-
1
2
∫
Ω
b
-
a
Δ
u
a
+
z
Δ
u
a
-
z
+
ɛ
DR
a
,
b
-
a
≥
0
,
f
o
r
e
v
e
r
y
b
∈
K
.Proof.
SinceJ
ɛ
(
a
)
≥
0, for every a
∈
K, we can find a minimizing sequence {
a
n
}
⊂
K such that(20)
lim
n
→
∞
⁡
J
ɛ
a
n
=
inf
a
∈
K
⁡
J
ɛ
a
.
The inequality(21)
ɛ
2
a
n
H
2
Ω
2
≤
1
2
∫
Ω
a
n
Δ
u
a
n
-
z
Δ
u
a
n
-
z
+
ɛ
2
a
n
H
2
Ω
2confirms that the sequence {
a
n
} is bounded in ·
H
2
(
Ω
). By the compact embedding of H
2
(
Ω
) into L
∞
(
Ω
), there exists a subsequence that converges weakly in H
2
(
Ω
) and strongly in L
∞
(
Ω
). By using the same notation for the subsequences as well, we obtain that a
n
⇀
a
~
∈
K in H
2
(
Ω
) and a
n
→
a
~ in L
∞
(
Ω
). Let u
n
=
u
(
a
n
) be the corresponding sequence of the solutions of the weak form (5). It follows from (6) and (7) that {
u
n
} is uniformly bounded in H
2
(
Ω
) and consequently there exists a weakly convergent subsequence. Once again by using the same notion for the subsequences, we assume that u
n
⇀
u
~. By the definition of u
n, we have(22)
∫
Ω
a
n
Δ
u
n
Δ
v
=
∫
Ω
f
v
,
∀
v
∈
V
,which, in view of the convergence properties of a
n and u
n, after passing to the limit implies that(23)
∫
Ω
a
~
Δ
u
~
Δ
v
=
∫
Ω
f
v
,where v
∈
V is arbitrary. Since the solution to the above variational problem is unique, we deduce that u
~
=
u
(
a
~
). In fact, u
n converges to u
~ strongly in H
2
(
Ω
). From the identities,(24)
∫
Ω
a
n
Δ
u
n
Δ
v
=
∫
Ω
f
v
,
∀
v
∈
V
,
∫
Ω
a
~
Δ
u
~
Δ
v
=
∫
Ω
f
v
,
∀
v
∈
V
.After a rearrangement of the terms, we deduce that(25)
∫
Ω
a
n
Δ
u
n
-
u
~
Δ
v
=
∫
Ω
a
~
-
a
n
Δ
u
~
Δ
v
,and the strong convergence follows at once from (6) and (7).
We next claim that(26)
∫
Ω
a
n
Δ
u
n
-
z
Δ
u
n
-
z
⟶
∫
Ω
a
~
Δ
u
a
~
-
z
Δ
u
a
~
-
z
.
For this, we notice that the following identities hold due to the definitions ofu
n, a
n, u
~, and a
~:(27)
∫
Ω
a
n
Δ
u
n
-
z
Δ
u
n
-
z
=
∫
Ω
f
u
n
-
z
-
∫
Ω
a
n
Δ
z
Δ
u
n
-
z
,
∫
Ω
a
~
Δ
u
a
~
-
z
Δ
u
a
~
-
z
=
∫
Ω
f
u
-
z
-
∫
Ω
a
~
Δ
z
Δ
u
a
~
-
z
.
For (26), we note that ∫
Ω
f
(
u
n
-
z
)
→
∫
Ω
f
(
u
~
-
z
) and ∫
Ω
a
n
Δ
z
Δ
(
u
n
-
z
)
→
∫
Ω
a
~
Δ
z
Δ
(
u
(
a
~
)
-
z
).
Finally, we have(28)
J
ɛ
a
~
=
1
2
∫
Ω
a
~
Δ
u
a
~
-
z
Δ
u
a
~
-
z
+
ɛ
2
a
~
H
2
Ω
2
=
lim
n
→
∞
⁡
1
2
∫
Ω
a
n
Δ
u
n
-
z
Δ
u
n
-
z
+
ɛ
2
a
~
H
2
Ω
2
≤
liminf
n
→
∞
⁡
J
ɛ
a
n
=
inf
⁡
J
ɛ
a
,
a
∈
K
,where we used the weak-lower semicontinuity of ·
H
2
(
Ω
)
. The uniqueness follows from the fact that the regularizer is strictly convex. Finally, since J
ɛ
(
a
) is a convex functional, the variational inequality is a necessary and sufficient optimality condition. The proof is complete.The continuous problem (18) has to be discretized to obtain a numerical solution. In this work, we use the finite element discretization on a triangulation T
h of Ω. We choose B
n to be the finite-dimensional space of the parameter space B
≔
H
2
(
Ω
) relative to T
h
. Analogously, we assume that {
V
n
} is a sequence of finite-dimensional subspaces of V and, for each n, P
n
:
V
→
V
n is a projection operator. We suppose that V
n and P
n have the property that, for every v
∈
V,(29)
v
-
P
n
v
V
⟶
0
as
n
⟶
∞
.We defineK
n
=
B
n
∩
K and assume that(30)
⋂
n
=
1
∞
K
n
≠
∅
.To keep the structure of finite-dimensional set of admissible sets quite general and also to allow other regularizers, we impose approximation properties given by the following assumption: for anya
∈
K, there exists a sequence {
a
n
} such that(31a)
a
n
∈
K
n
∀
n
,
(31b)
a
n
⟶
a
in
L
∞
Ω
,
(31c)
R
a
n
⟶
R
a
.The above property is an analogue of (29) for the set of feasible constraints.We now defineJ
(
n
)
:
K
n
→
R as the following discrete analogue of (8):(32)
J
n
a
=
1
2
∫
Ω
a
Δ
u
n
a
-
P
n
z
Δ
u
n
a
-
P
n
z
,where u
n
(
a
) is the solution of the following discrete variational problem:(33)
∫
Ω
a
Δ
u
n
a
Δ
v
=
∫
Ω
f
v
,
for  every
v
∈
V
n
.Here, the dataz is replaced by the projected data P
n
z for computational convenience. Since P
n
z converges to z strongly in V, all of the following results hold irrespective of whether z or P
n
z is used in the definition of J
(
n
). We then define the following regularized analogue of the above discrete MOLS(34)
J
ɛ
n
a
=
J
n
a
+
ɛ
R
aand solve(35)
min
a
∈
K
n
⁡
J
ɛ
n
a
.The solvability of the above regularized problems can easily be proved by arguments used above.We can now prove the following convergence result.Theorem 5.
Assume that, for eachn, a
n
∗
∈
K
n is a solution of (35). Then, any weak H
2
(
Ω
) limit point of the sequence {
a
n
∗
} is a solution of (18).Proof.
By (30), there exists a constant C such that J
ɛ
(
n
)
(
a
n
∗
)
≤
C, for all n. Therefore, the sequence {
a
n
∗
} is bounded in H
2
(
Ω
) and hence it has a subsequence that converges weakly in H
2
(
Ω
) and strongly in L
∞
(
Ω
)
. For simplicity of notation, we assume that a
n
∗
⇀
a
∗
∈
K in H
2
(
Ω
) and a
n
∗
→
a
∗ in L
∞
(
Ω
)
. We will show that a
∗ is a solution of (18). Let u
n
=
u
(
a
n
) be the corresponding solution of (33). Evidently, {
u
n
} is bounded and hence there is a subsequence that weakly converges to some u
∗ in V. We claim that u
∗
=
u
(
a
∗
)
. In view of the definition of u
n, we have(36)
∫
Ω
a
n
∗
Δ
u
n
Δ
v
n
=
∫
Ω
f
v
n
,
for  every
v
n
∈
V
n
.
Letv
∈
V be arbitrary. We set v
n
=
P
n
v in the above identity to get(37)
∫
Ω
a
n
∗
Δ
u
n
Δ
P
n
v
=
∫
Ω
f
P
n
v
,which after a rearrangement of terms implies that(38)
∫
Ω
f
P
n
v
=
∫
Ω
a
∗
Δ
u
n
Δ
v
+
∫
Ω
a
n
∗
-
a
∗
Δ
u
n
Δ
v
-
∫
Ω
a
n
∗
Δ
u
n
Δ
P
n
v
-
v
, and after passing to the limit n
→
∞, we obtain that (39)
∫
Ω
a
∗
Δ
u
a
∗
Δ
v
=
∫
Ω
f
v
,and because v was chosen arbitrarily, we obtain that u
∗ solves (5). That is, u
∗
=
u
(
a
∗
)
.
By using the arguments given in Theorem4, we can also show that(40)
∫
Ω
a
n
∗
Δ
u
n
-
z
Δ
u
n
-
z
⟶
∫
Ω
a
∗
Δ
u
a
∗
-
z
Δ
u
a
∗
-
z
,which, in view of the lower semicontinuity of R
(
·
)
=
·
H
2
(
Ω
)
2, implies that(41)
J
ɛ
a
∗
≤
liminf
n
→
∞
⁡
J
ɛ
n
a
n
∗
.
Now consider an arbitrarya
∈
K. By (31a), (31b), and (31c), there exists a sequence {
a
n
} such that a
n
∈
K
n, for all n, a
n
→
a in L
∞, and R
(
a
n
)
→
R
(
a
). This, however, implies that(42)
lim
n
→
∞
⁡
J
ρ
n
a
n
=
J
a
,and hence, by the optimality of a
n
∗, we have(43)
J
ɛ
a
∗
≤
liminf
n
→
∞
⁡
J
ɛ
n
a
n
∗
≤
liminf
n
→
∞
⁡
J
ɛ
n
a
n
=
lim
n
→
∞
⁡
J
ɛ
n
a
n
=
J
ɛ
a
.
Sincea
∈
K was arbitrary, this shows that a
∗ is a solution of (8), which completes the proof.


## 3. Discrete Formulas

To represent the discrete weak form in a computable way, we proceed as follows. We represent bases forB
n and V
n by ψ
1
,
ψ
2
,
…
,
ψ
m and φ
1
,
φ
2
,
…
,
φ
n, respectively. The space B
n is then isomorphic to R
m, and for any a
∈
B
n, we define A
∈
R
m by A
i
=
a
(
x
i
), i
=
1,2
,
…
,
m, where {
ψ
1
,
ψ
2
,
…
,
ψ
m
} is a nodal basis corresponding to the nodes {
x
1
,
x
2
,
…
,
x
m
}. Conversely, each A
∈
R
m corresponds to a
∈
B
h defined by a
=
∑
i
=
1
m
A
i
ψ
i. Similarly, u
∈
V
n will correspond to U
∈
R
n, where U
i
=
u
(
y
i
), i
=
1,2
,
…
,
n and u
=
∑
i
=
1
n
U
i
φ
i. Here, y
1
,
y
2
,
…
,
y
n are the nodes of the mesh defining V
n.Then, the stiffness matrixK
(
A
) and the load vector F are given by(44)
K
A
i
j
=
∫
Ω
a
Δ
φ
i
Δ
φ
j
,
i
,
j
=
1,2
,
…
,
n
,
F
i
=
∫
Ω
f
φ
i
,
i
=
1,2
,
…
,
n
.We also define the matrixL
(
V
) by the following condition:(45)
L
V
A
=
K
A
V
∀
A
∈
R
m
,
V
∈
R
n
.Using the above notation, the following formulas can be verified:(46)
J
A
=
1
2
U
-
Z
K
A
U
-
Z
;
D
J
A
=
-
1
2
L
U
+
Z
T
U
-
Z
;
D
2
J
A
=
L
U
T
K
A
-
1
L
U
.


## 4. Continuous Newton Method

After the discretization, the finite-dimensional minimization problem can be solved by any known optimization solver. In this paper, we will employ the continuous Newton method recently proposed by Zhang et al. [22]. Continuous methods for optimization problems solve an optimization problem by following the solution trajectory of a system of ordinary differential equations. Although numerous authors have focused on continuous optimization methods for solving unconstrained and constrained optimization problems, most of these works have been of a theoretical nature. However, recently, many researchers have advocated the usefulness of continuous methods for solving large scale optimization problems and have given significant numerical evidence to show the competitiveness of these methods with the conventional optimization methods (see [23–26]).In this work, we will use the continuous Newton method [22] to solve the finite-dimensional optimization problems given by the regularized MOLS and OLS. For simplicity, the side constraints will be ignored. Zhang et al. [22] propose following the trajectory given by the following initial value problem:(47)
d
A
d
t
=
G
A
,
A
t
0
=
A
0
,where G
(
A
) is defined by(48)
G
A
=
-
∇
2
J
A
t
-
1
∇
J
A
t
if
λ
min
⁡
A
>
δ
2
-
α
t
∇
2
J
A
t
-
1
∇
J
A
t
-
β
t
∇
J
A
t
if
δ
1
≤
λ
min
⁡
A
≤
δ
2
-
∇
J
A
t
if
λ
min
⁡
A
<
δ
1
.Here,J corresponds to either the regularized MOLS or the regularized OLS, the minimum eigenvalue of the Hessian of J
(
A
) is given by λ
min
⁡
(
A
), and δ
2
>
δ
1
>
0. Moreover,(49)
α
a
=
λ
min
⁡
A
-
δ
1
δ
2
-
δ
1
,
β
a
=
1
-
α
A
=
δ
2
-
λ
min
⁡
A
δ
2
-
δ
1
.The idea of the proposed scheme is to combine the steepest descent and Newton’s direction through a convex combination. It involves a simple yet efficient strategy to get around the potential singularities of the Hessian. To be specific, whenλ
min
⁡
(
a
) is between δ
1 and δ
2, if its value is closer to δ
2, then α
(
a
) will be larger, in turn putting more emphasis on the Newton direction. On the other hand, if λ
min
⁡
(
a
) is closer to δ
1, then β
(
a
) will be larger, putting more weight on the gradient direction. The approach is particularly appealing as one of the main challenges in the nonlinear inverse problems is the ill-posedness of the Hessian. Moreover, since the continuous Newton method takes into account the behaviour of the Hessian by examining its minimum eigenvalue, we also gain insight into the relative convexity (or nonconvexity) of the OLS or MOLS functionals. Finally, this approach places less emphasis on computing an optimal regularization parameter, a difficult task in inverse problems.


## 5. Numerical Experiments

In this section, we present the results of numerical experiments designed to gauge the preliminary effectiveness of the MOLS approach when coupled with the continuous Newton method and applied to the flexural rigidity inverse problem.


### 5.1. Experiment Setup

We consider first an example boundary value problem (Example 1) derived from (1a), (1b), and (1c):(50)
Δ
a
x
,
y
Δ
u
x
,
y
=
f
x
,
y
in
Ω
,
u
x
,
y
=
∂
u
∂
n
=
0
on
Γ
,where the solution u and the exact parameter a are given by(51)
u
x
,
y
=
100
x
2
1
-
x
2
y
2
1
-
y
2
,
a
x
,
y
=
3
2
+
sin
⁡
π
x
sin
⁡
π
yand where f
(
x
,
y
) is subsequently defined by (50). The domain Ω is taken as the unit square, Ω
=
(
0,1
)
×
(
0,1
), with the boundary Γ as the square’s outside edges.Additionally, we considered a similar example (Example 2) with the recovery of the same parameter as Example 1, but on a simply supported circular plate (with appropriate changes to the forward problem applied).Finally, we considered an example (Example 3) again on the unit square, but with the recovery of a different parameter and with more complicated mixed boundary conditions (i.e., clamped along the top and bottom edges, simply supported along the left edge, and with the right edge free). The parameter recovered in this case was given by(52)
a
x
,
y
=
3
2
x
2
+
y
2
+
x
y
.Discretization of the solution was performed using cubic Hermite elements for all examples. For Examples 1 and 3, the mesh used was a regular 25 × 25 triangular grid consisting of 1,352 triangles resulting in 3,529 total degrees of freedom for the solution. For the circular domain in Example 2, a mesh was generated using a Delaunay triangulation yielding 848 triangles with corresponding 2,219 degrees of freedom.


### 5.2. Algorithm Details

The optimization was performed using the continuous Newton method outlined in the previous section applied to both the discretized MOLS and OLS functionals. Additionally, for means of comparison, the optimization was also performed using a conjugate gradient trust-region-reflective method (MATLAB’sfminunc) using similar stopping criteria. In all examples, H
1 seminorm regularization was used for computational simplicity with a fixed regularization parameter used per example (see Table 1 for values).Table 1
Performance results for the MOLS and OLS methods.
Method
Algorithm
Iterations
L
2-error
Min.λ
min
⁡
Max.λ
min
⁡
Example  1:ε
=
1
0
-
6, δ
1
=
1
0
-
8, and δ
2
=
3
×
1
0
-
7
MOLS
Continuous Newton
95
5.53e − 05
1.2981e − 06
1.6089e − 06
f
minunc
71
2.71e − 04
—
—
OLS
Continuous Newton
184
3.32e − 02
−5.6394e − 06
4.9100e − 07
fminunc
77
5.39e − 02
—
—
Example  2:ε
=
1
0
-
5, δ
1
=
1
0
-
8, and δ
2
=
1
0
-
6
MOLS
Continuous Newton
26
1.67e − 04
3.9702e − 03
7.2222e − 02
fminunc
63
2.68e − 06
—
—
OLS
Continuous Newton
662
2.44e − 02
2.9112e − 06
4.3008e − 03
fminunc
97
2.44e − 02
—
—
Example  3:ε
=
1
0
-
8, δ
1
=
1
0
-
9, and δ
2
=
3
×
1
0
-
8
MOLS
Continuous Newton
112
2.98e − 03
1.2992e − 08
1.3792e − 08
fminunc
79
2.65e − 03
—
—
OLS
Continuous Newton
442
5.13e − 01
5.5219e − 09
1.1703e − 08
fminunc
69
1.22e − 01
—
—In the continuous Newton algorithm, the Lanczos iteration was used to compute/estimate the minimum eigenvalue of∇
2
J and in those cases where the Lanczos method failed to converge in the given (limited) number of iterations, the method defaulted to the steepest descent trajectory. The resulting ordinary differential equation was solved using MATLAB’sode113 solver with default parameters. The ODE was solved over the time domain [
0
,
T
] with T
=
100 for MOLS and T
=
1000 for OLS. The additional continuous Newton algorithmic parameters were fixed for each example (again, see Table 1 for values).


### 5.3. Numerical Results

Figures1–3, respectively, show the computed parameter at several algorithm steps for the OLS and MOLS approaches using both the continuous Newton andfminunc algorithms. These figures also show the final output of the algorithm and the error between the computed and exact parameters.Figure 1
Example 1: MOLS, continuous Newton.Figure 2
Example 2: MOLS, continuous Newton.Figure 3
Example 3: MOLS, continuous Newton.In Figure4, we compare the minimum eigenvalues of the Hessian of both the MOLS and OLS functional over a number of continuous Newton algorithm iterations applied to Example 1. Here, we see that, for the OLS method, the algorithm transitions to both the “mixed” step direction and, when the minimum eigenvalue becomes exceedingly small or negative, the steepest descent direction. We also note that, for the same regularization parameter, the MOLS method uses the full Newton direction. The minimum and maximum of λ
min
⁡ along with the L
2 error and the total number of algorithm iterations for all examples are provided in Table 1.Figure 4
Comparison ofλ
min
⁡ for Example 1 using MOLS and OLS methods.


## 5.1. Experiment Setup

We consider first an example boundary value problem (Example 1) derived from (1a), (1b), and (1c):(50)
Δ
a
x
,
y
Δ
u
x
,
y
=
f
x
,
y
in
Ω
,
u
x
,
y
=
∂
u
∂
n
=
0
on
Γ
,where the solution u and the exact parameter a are given by(51)
u
x
,
y
=
100
x
2
1
-
x
2
y
2
1
-
y
2
,
a
x
,
y
=
3
2
+
sin
⁡
π
x
sin
⁡
π
yand where f
(
x
,
y
) is subsequently defined by (50). The domain Ω is taken as the unit square, Ω
=
(
0,1
)
×
(
0,1
), with the boundary Γ as the square’s outside edges.Additionally, we considered a similar example (Example 2) with the recovery of the same parameter as Example 1, but on a simply supported circular plate (with appropriate changes to the forward problem applied).Finally, we considered an example (Example 3) again on the unit square, but with the recovery of a different parameter and with more complicated mixed boundary conditions (i.e., clamped along the top and bottom edges, simply supported along the left edge, and with the right edge free). The parameter recovered in this case was given by(52)
a
x
,
y
=
3
2
x
2
+
y
2
+
x
y
.Discretization of the solution was performed using cubic Hermite elements for all examples. For Examples 1 and 3, the mesh used was a regular 25 × 25 triangular grid consisting of 1,352 triangles resulting in 3,529 total degrees of freedom for the solution. For the circular domain in Example 2, a mesh was generated using a Delaunay triangulation yielding 848 triangles with corresponding 2,219 degrees of freedom.


## 5.2. Algorithm Details

The optimization was performed using the continuous Newton method outlined in the previous section applied to both the discretized MOLS and OLS functionals. Additionally, for means of comparison, the optimization was also performed using a conjugate gradient trust-region-reflective method (MATLAB’sfminunc) using similar stopping criteria. In all examples, H
1 seminorm regularization was used for computational simplicity with a fixed regularization parameter used per example (see Table 1 for values).Table 1
Performance results for the MOLS and OLS methods.
Method
Algorithm
Iterations
L
2-error
Min.λ
min
⁡
Max.λ
min
⁡
Example  1:ε
=
1
0
-
6, δ
1
=
1
0
-
8, and δ
2
=
3
×
1
0
-
7
MOLS
Continuous Newton
95
5.53e − 05
1.2981e − 06
1.6089e − 06
f
minunc
71
2.71e − 04
—
—
OLS
Continuous Newton
184
3.32e − 02
−5.6394e − 06
4.9100e − 07
fminunc
77
5.39e − 02
—
—
Example  2:ε
=
1
0
-
5, δ
1
=
1
0
-
8, and δ
2
=
1
0
-
6
MOLS
Continuous Newton
26
1.67e − 04
3.9702e − 03
7.2222e − 02
fminunc
63
2.68e − 06
—
—
OLS
Continuous Newton
662
2.44e − 02
2.9112e − 06
4.3008e − 03
fminunc
97
2.44e − 02
—
—
Example  3:ε
=
1
0
-
8, δ
1
=
1
0
-
9, and δ
2
=
3
×
1
0
-
8
MOLS
Continuous Newton
112
2.98e − 03
1.2992e − 08
1.3792e − 08
fminunc
79
2.65e − 03
—
—
OLS
Continuous Newton
442
5.13e − 01
5.5219e − 09
1.1703e − 08
fminunc
69
1.22e − 01
—
—In the continuous Newton algorithm, the Lanczos iteration was used to compute/estimate the minimum eigenvalue of∇
2
J and in those cases where the Lanczos method failed to converge in the given (limited) number of iterations, the method defaulted to the steepest descent trajectory. The resulting ordinary differential equation was solved using MATLAB’sode113 solver with default parameters. The ODE was solved over the time domain [
0
,
T
] with T
=
100 for MOLS and T
=
1000 for OLS. The additional continuous Newton algorithmic parameters were fixed for each example (again, see Table 1 for values).


## 5.3. Numerical Results

Figures1–3, respectively, show the computed parameter at several algorithm steps for the OLS and MOLS approaches using both the continuous Newton andfminunc algorithms. These figures also show the final output of the algorithm and the error between the computed and exact parameters.Figure 1
Example 1: MOLS, continuous Newton.Figure 2
Example 2: MOLS, continuous Newton.Figure 3
Example 3: MOLS, continuous Newton.In Figure4, we compare the minimum eigenvalues of the Hessian of both the MOLS and OLS functional over a number of continuous Newton algorithm iterations applied to Example 1. Here, we see that, for the OLS method, the algorithm transitions to both the “mixed” step direction and, when the minimum eigenvalue becomes exceedingly small or negative, the steepest descent direction. We also note that, for the same regularization parameter, the MOLS method uses the full Newton direction. The minimum and maximum of λ
min
⁡ along with the L
2 error and the total number of algorithm iterations for all examples are provided in Table 1.Figure 4
Comparison ofλ
min
⁡ for Example 1 using MOLS and OLS methods.


## 6. Discussion

The results of the numerical experiments suggest the following:(i)
The MOLS functional converges faster than the OLS functional (Table1) for the flexural rigidity inverse problem under similar conditions (i.e., similar regularization parameters, optimization algorithm).
(ii)
In at least one instance for the MOLS approach (see Example 2 in Table1), the continuous Newton method shows accelerated convergence over thefminunc algorithm.
(iii)
The eigenvalues of the Hessian of the MOLS functional remain strictly positive compared with the Hessian of the OLS functional, whose minimum eigenvalues can become negative after a number of algorithm iterations (see both Table1 and Figure 4).
(iv)
The more complicated domains and boundary conditions in Examples 2 and 3 have a negative impact in terms of algorithm performance and accuracy in recovery of the flexural rigidity parameter. As can be seen in Table1, this is particularly true in the case of the OLS approach.
(v)
Overall, for either the MOLS or OLS approach, the error in the recovery of the coefficient is most prominent along the boundaries with significantly better reconstruction in the domain’s interior. Again, in particular, the OLS approach seems to have the most difficulty in reconstructing the flexural rigidity along the boundaries.


## 7. Concluding Remarks

Further directions for research are as follows:(i)
Thorough comparison of the performance and stability differences between the OLS and MOLS approaches.
(ii)
Thorough comparison of differing approaches to solve the ODE system arising from the continuous Newton method including a comparison of the relative performance and stability of calculating/estimating the eigenvalues of the Hessian ofJ.
(iii)
Comparison/exploration of differing regularization approaches.
(iv)
The use of mixed methods for solving the variational problem.



---

*Source: 290301-2015-09-08.xml*