# On the Iterative Methods of Linearization, Decrease of Order and Dimension of the Karman-Type PDEs

Authors: A. V. Krysko; J. Awrejcewicz; S. P. Pavlov; M. V. Zhigalov; V. A. Krysko
Journal: The Scientific World Journal
 (2014)

---

## Abstract

Iterative methods to achieve a suitable linearization as well as a decrease of the order and dimension of nonlinear partial differential equations of the eighth order into the biharmonic and Poisson-type differential equations with their simultaneous linearization are proposed in this work. Validity and reliability of the obtained results are discussed using computer programs developed by the authors.


---

## Body

## 1. Introduction

Mathematical models of continuous mechanical structures are described by nonlinear partial differential equations which may be solved analytically only in a few rare cases. However, a direct application of the numerical methods is associated also with big difficulties regarding a high order of both dimension and differential operator, as well as nonlinearity of the PDEs studied.This is why it is tempting to develop approaches that offer a reduction of the input differential equations. The mentioned methods can be divided into three groups:(
1
) linearization; (
2
) order decrease of the PDEs; (
3
) order decrease of a differential operator.The so far existing methods of solutions of nonlinear problems, depending on the introduced linearization level, can be divided into two groups. The first one deals with the linearization of PDEs, whereas the second one is dedicated to the linearization of algebraic equations obtained through the discretization procedures applied to the input PDEs. Below, we consider the methods associated with the first group. This group contains the Newton and Newton-Kantorovich methods [1].One of the linearization methods is the method of quasilinearization, widely illustrated in monograph [2]. It presents a further development of Newton’s method, and it generalizes the method proposed by Kantorovich.On the other hand, there is a seminal approach known as the Agmon-Douglis-Nirenberg (ADN) theory for elliptic PDEs still attracting a large number of imitators [3, 4]. In particular, the abstract least squares theory is developed satisfying the ADN elliptic theory assumptions [5–7].Furthermore, in the case of corners in plane domains the ADN system exhibits singularities, which imply a need for construction of singular exponents and angular functions [8]. Our approach does not have this disadvantage and it is simple in direct applications to the real world systems.The so far briefly addressed approaches linearize the input problem; that is, they reduce it to the solution of linear problems. However, there is one more important question to be solved, that is, a reduction of the space dimension of the initial problem.One of the methods to solve the stated problem is focused on averaging (integration) along such a coordinate on which the object dimension is lesser in comparison to the two remaining coordinates. On the other hand, it is well known that mathematical problems related to the theory of material strength can be formulated as variational problems, that is, the problems of finding extrema of a certain functional. Variational statements create a foundation for the construction of direct difference and variational methods, as it is widely described in monograph [9].We mention only a few works [10–12] devoted to the third group, that is, aiming at a decrease of the PDE order.Note that the so far presented state of the art of the proposed and applied methods allows us to solve each of the mentioned problems separately: either a decrease of the system order or its linearization. However, we show how all these problems can be solved simultaneously.Our paper is focused mainly on the method of dimension decrease and linearization of the Karman-type PDEs. However, the presented approach can be successfully applied to other nonlinear PDEs. In particular, in the modified version two variants of the proposed method are presented:(i)
the first iterative method consists of the reduction of the eighth order linear PDEs to a successive solution of linear PDEs of the fourth order biharmonic equations; that is, the system dimension is reduced twice with simultaneous linearization of the problem;
(ii)
the applied second iterative procedure includes a further order decrease of the earlier obtained (first iterative method) linear system of biharmonic PDEs of the fourth order to the successive solution to the system of the second order Poisson-type equations.In other words, the application of these two iterative procedures implies a fourfold reduction of the PDEs order with the linearization procedure carried out simultaneously.The proposed iterative procedures regarding the nonlinear PDEs order decrease and linearization can also be applied to PDEs with a curvilinear boundary. The application of FDM (finite difference method) to solve biharmonic equations and PDEs of the Poisson-type requires a solution to the so-called Sapondzhyan-Babuška problem. The paradox of Sapondzhyan-Babuška (see [13–15]) was discovered when studying the asymptotic behavior of solutions to an elasticity system in a thin polygonal plate (inscribed in the plate with smooth boundary) as the length of the side of the polygon tends to zero and the number of sides goes to infinity. In Section 2 of our work we prove the proposed iterative procedure to remove this paradox (this problem concerns smoothness of the curvilinear boundary).In Section3 of our work the reliability and validity of the method of variational iterative procedure to solve PDEs described by positively defined operators are illustrated and discussed. Namely, the convergence of the method of variational iterations generalizes the Kantorovich-Vlasov method [16] aimed at the reduction of PDEs to ODEs. On the other hand, as it was pointed out by Vorovich [17], the Kantorovich-Vlasov method generalizes the Galerkin method. It should be emphasized that the choice of approximating functions referring either to two variables in the Galerkin method or to one variable in the Kantorovich-Vlasov method cannot be applied in the method of variational iterations. The system of functions being sought is provided by a solution of the PDEs with regard to two variables assuming that we deal with the 2D problem. Furthermore, the proposed approach can be applied to 3D elliptic equations.Section6 of the paper deals with a comparison of the solutions to the Karman equations obtained via our proposed iterative procedures with those offered by FEM and FDM, as well as with experimental results. Good coincidence of the results is achieved.


## 2. Mathematical Model of a Flexible Karman-Type Plate (Hypotheses, Differential Equations, and Boundary Conditions)

The objects of our investigation are plates of different shapes (in particular, rectangular ones), representing a closed 3D part of spaceR
3 (Figure 1). The following hypotheses are introduced: (i) plate material is elastic and isotropic; (ii) the following Karman relations between deformations and displacements are introduced:
(1)
ε
x
x
=
∂
u
∂
x
+
1
2
(
∂
w
∂
x
)
2
,
ε
x
y
=
∂
u
∂
y
+
∂
v
∂
x
+
∂
w
∂
x
∂
w
∂
y
,
(
u
,
v
,
w
)
,
(
x
,
y
)
.Figure 1
A studied plate.Equations governing the deflectionw
(
x
,
y
,
t
) and stress function F
(
x
,
y
) have the following form [15]:
(2)
Δ
2
w
-
L
(
w
,
F
)
-
q
=
0
,
Δ
2
F
+
1
2
L
(
w
,
w
)
=
0
.The following operators are introduced:(3)
Δ
2
(
·
)
=
1
λ
2
∂
4
(
·
)
∂
x
4
+
λ
2
∂
4
(
·
)
∂
y
4
+
2
∂
4
(
·
)
∂
x
2
∂
y
2
,
L
(
·
,
·
)
=
∂
2
(
·
)
∂
x
2
∂
2
(
·
)
∂
y
2
-
2
∂
2
(
·
)
∂
x
∂
y
∂
2
(
·
)
∂
x
∂
y
+
∂
2
(
·
)
∂
x
2
∂
2
(
·
)
∂
y
2
.Here and further on the following nondimensional quantities are introduced:w
=
h
w
-; F
=
E
h
2
F
¯; t
=
t
0
t
¯; ε
=
ε
¯
/
τ; x
=
a
x
¯; y
=
b
y
-; q
=
q
¯
(
12
(
1
-
ν
2
)
E
h
4
/
a
2
b
2
); τ
=
(
a
b
/
h
)
ρ
/
E
g; λ
=
a
/
b, where a
,
b are the maximal plate dimensions regarding x and y, respectively; h is thickness; g is acceleration due to gravity; ρ
=
γ
h; γ is specific gravity of volume plate material; ν is Poisson’s coefficient; E is the Young modulus; w
,
F are deflection and stress functions, respectively.Let us add boundary conditions of the support on flexible nonstretched (noncompressed) ribs to the system of plates [18, 19]:
(4)
w
|
Γ
=
∂
2
w
∂
n
2
|
Γ
=
F
|
Γ
=
∂
2
F
∂
n
2
|
Γ
=
0
,
where Γ stands for the space boundary occupied by the plate. The following initial conditions are attached to (2):
(5)
w
|
t
=
0
=
w
0
,
w
˙
|
t
=
0
=
w
˙
0
.System of (2) is composed of nonlinear PDEs of the eighth order. Finding a reliable solution to this problem is still a serious problem in spite of achievements of the numerical methods. It should be emphasized that the solution to the mentioned problem was found earlier via either FDM (finite difference method) or FEM (finite element method), or by the Bubnov-Galerkin method. Below, we propose a novel method of order reduction and linearization of PDEs (2).


## 3. Methods of Order Decrease and Linearization of the Karman Equation

There are two ways for construction of the fundamental iterative procedure to solve system (2): (i) system reduction to a successive solution to the Germain-Lagrange type equations, in this case the system order is decreased twice; (ii) system reduction to a Poisson-type equation (in this case the system order is reduced four times). In both mentioned cases the linearization procedure of the input PDEs is carried out.


### 3.1. Iterative Linearization Procedure and Reduction of the Karman Equation into Germain-Lagrange Equations

We keep the biharmonic operator in each of (2), and we shift nonlinear terms into their right-hand sides. Assuming that functions on the right-hand sides are computed with respect to their previous step and that the equations are solved successively, the following iterative procedure is proposed:
(6)
Δ
2
w
(
k
)
=
L
(
w
(
k
-
1
)
,
F
(
k
-
1
)
)
+
q
,
Δ
2
F
(
k
)
=
-
1
2
L
(
w
(
k
)
,
w
(
k
)
)
,
{
x
,
y
}
∈
Ω
.On the first step of the iterative procedure the following biharmonic equation for a given loadq
(
x
,
y
) is solved:
(7)
Δ
2
w
(
1
)
(
x
,
y
)
=
q
(
x
,
y
)
.The value ofw
(
1
)
(
x
,
y
) is substituted into the right-hand side of equation system (6), and as a result a biharmonic equation for F
(
1
)
(
x
,
y
) with the known right-hand side is obtained. The value of the stress function found so far is substituted to the first system equation. The process of finding solutions is continued to achieve the required accuracy.Let us note that as a result of the application of the iterative procedure, the Germain-Lagrange type system of equations are obtained.Let us prove convergence of the constructed iterative procedure. LetH
2
(
Ω
) refer to a Sobolev space of functions ξ
=
{
w
,
F
} such that
(8)
ξ
∈
L
2
(
Ω
)
,
∂
ξ
∂
x
i
∈
L
2
(
Ω
)
,
∂
2
ξ
∂
x
i
∂
x
j
∈
L
2
(
Ω
)
;
i
,
j
=
1,2
,
where L
2
(
Ω
) denotes the space of functions being summed with a square in Ω.LetH
0
2
(
Ω
) denote the closure of functions from D
(
Ω
) (space of functions of class C
∞ in Ω, having compact carrier in Ω) in norm H
2
(
Ω
):
(9)
H
0
2
(
Ω
)
=
D
(
Ω
)
¯
H
2
(
Ω
)
=
{
ξ
∈
H
2
(
Ω
)
∣
ξ
|
Γ
=
∂
ξ
∂
n
|
Γ
=
0
}
.Since spaceΩ is bounded and its boundary Γ is efficiently regular, then map ξ
→
∥
Δ
ξ
∥
0
,
Ω defines the norm in H
0
2
(
Ω
) being equivalent to the norm generated by spaces H
2
(
Ω
).Assume thatq
∈
H
-
2
(
Ω
) (H
-
2
(
Ω
) denotes a conjugation to H
2
(
Ω
)). It is known [17] that in this case problems (2) and (4) have a solution (it may happen that it shall be nonunique).Novel Variation Formulation of the Problem. Let us denote by (
·
,
·
) a scalar product in L
2
(
Ω
): (
ξ
,
η
)
=
∫
Ω
ξ
η
d
Ω, and by β
(
w
,
F
,
μ
) a three linear form defined on (
H
0
2
(
Ω
)
)
3:
(10)
β
(
w
,
F
,
μ
)
=
(
Δ
F
,
Δ
μ
)
+
1
2
(
L
(
w
,
w
)
,
μ
)
.Let us define the set(11)
M
=
{
w
,
F
∈
H
0
2
(
Ω
)
∣
∀
μ
∈
H
0
2
(
Ω
)
,
β
(
w
,
F
,
μ
)
=
0
}
,
and square function J
(
w
,
F
)
:
M
→
R
(12)
J
(
w
,
F
)
=
1
2
∥
Δ
w
∥
0
,
Ω
2
+
1
2
∥
Δ
F
∥
0
,
Ω
2
-
(
q
,
w
)
.Theorem 1.
The problem of minimizing (12) on set (11) has, at least, one solution.Proof.
Let{
w
n
,
F
n
}
∈
M be the minimizing series; that is, we have
(13)
J
(
w
n
,
F
n
)
⟶
inf
⁡
{
w
,
F
}
∈
M
J
(
w
,
F
)
,
which exists, since J is a square functional.
For arbitraryw
,
F
∈
H
0
2
(
Ω
) the following inequality holds:
(14)
J
(
w
,
F
)
≥
c
1
∥
w
∥
2
,
Ω
2
+
c
2
∥
F
∥
2
,
Ω
2
-
c
3
·
∥
w
∥
2
,
Ω
,
where ∥
·
∥
2
,
Ω denotes the norm in H
2
(
Ω
) and c
i are certain positive constants. Then, (13) yields c
1
∥
w
n
∥
2
,
Ω
2
+
c
2
∥
F
n
∥
2
,
Ω
2
-
c
3
∥
w
n
∥
2
,
Ω
≤
J
(
w
n
,
F
n
)
≤
J
(
w
0
,
F
0
)
=
A, where w
0
,
F
0 are the arbitrary functions (initial approximation).
Then, the following estimation holds:(15)
c
1
(
∥
w
n
∥
2
,
Ω
2
-
c
3
2
c
1
)
2
+
c
2
∥
F
n
∥
2
,
Ω
2
≤
A
+
c
3
2
4
c
1
,
∥
w
n
∥
2
,
Ω
≤
c
4
,
∥
F
n
∥
2
,
Ω
≤
c
5
.
Therefore, series{
w
n
,
F
n
} is bounded in (
H
0
2
(
Ω
)
)
2. Consequently, one may choose a series {
w
k
,
F
k
} that w
k
→
w
~, F
k
→
F
~ is weak in H
0
2
(
Ω
). Since H
0
2
(
Ω
)
→
L
2
(
Ω
) is compact, then w
k
→
w
~, F
k
→
F
~ is strong in L
2
(
Ω
).
We show that interval{
w
~
,
F
~
} of the minimized series belongs to M; that is, β
(
w
~
,
F
~
,
μ
)
=
0, for all μ
∈
H
0
2
(
Ω
).
Since(
L
(
w
k
,
w
k
)
,
μ
)
=
(
L
(
w
k
,
μ
)
,
w
k
)
∀
μ
∈
H
0
2
(
Ω
), L
(
w
k
,
μ
)
→
L
(
w
~
,
μ
) is weak in H
0
2
(
Ω
), and w
k
→
w
~ is strong in L
2
(
Ω
), we get (
L
(
w
k
,
w
k
)
,
μ
)
=
(
L
(
w
~
,
w
~
)
,
μ
) and consequently, β
(
w
~
,
F
~
,
μ
)
=
0 for all μ
∈
H
0
2
(
Ω
). This means that
(16)
{
w
~
,
F
~
}
∈
M
.
However,J
(
w
,
F
) is half-continuous from below in weak topology on (
H
2
(
Ω
)
)
2, and therefore the following inequality holds: lim
⁡
k
→
∞
¯
J
(
w
k
,
F
k
)
≥
J
(
w
~
,
F
~
).
Then (13) and (16) imply that J
(
w
~
,
F
~
)
≤
inf
⁡
(
w
,
F
)
∈
M
J
(
w
,
F
). Therefore, the following equation holds: J
(
w
~
,
F
~
)
=
inf
⁡
(
w
,
F
)
∈
M
J
(
w
,
F
), which means that {
w
~
,
F
~
}
∈
M is a solution to the minimization problem.Let us explain how points of the minimum of functional (12) are linked with solutions to problems (6) and (4). For this purpose a notation of weak solution shall be introduced.A weak solution to problems (6) and (4) is defined by the pair of functions {
w
,
F
}
∈
M, satisfying the following:
(17)
(
Δ
w
,
Δ
μ
)
-
(
L
(
w
,
F
)
,
μ
)
=
(
q
,
μ
)
∀
μ
∈
H
0
2
(
Ω
)
.Theorem 2.
Points of the functional minimum (12) are weak solutions to problems (6) and (4).Proof.
Let{
w
,
F
}
∈
M be one of the functional (12) minimum points. Let us take η
=
w
+
t
δ
w
δ
w
∈
H
0
2
(
Ω
) and let us choose ξ
=
F
+
δ
F
δ
F
∈
H
0
2
(
Ω
) such that {
η
,
ξ
}
∈
M, that is, in the way that β
(
w
,
F
,
μ
)
=
0 for all μ
∈
H
0
2
(
Ω
). Then J
(
w
,
F
)
≤
J
(
η
,
ξ
). This yields
(18)
t
(
Δ
w
,
Δ
δ
w
)
+
(
Δ
F
,
Δ
δ
F
)
-
t
(
q
,
δ
w
)
+
t
2
2
∥
Δ
δ
w
∥
2
,
Ω
2
+
1
2
∥
Δ
δ
F
∥
2
,
Ω
2
≥
0
,
∀
t
∈
R
,
δ
w
∈
H
0
2
(
Ω
)
,
and by taking μ
=
F, condition β
(
η
,
ξ
,
μ
)
=
0 yields
(19)
(
Δ
F
,
Δ
δ
F
)
=
-
t
(
L
(
w
,
δ
w
)
,
F
)
-
t
2
2
(
L
(
δ
w
,
δ
w
)
,
F
)
.
Substituting (19) into (18), dividing the obtained expression by t and going to the limit for t
→
0, the following inequality is obtained:
(20)
(
Δ
w
,
Δ
δ
w
)
-
(
L
(
w
,
F
)
,
δ
w
)
≥
(
q
,
μ
)
.
Substitutingδ
w by -
δ
w in (20), one obtains equality, that is, (17).
Let us use the following notation(
Φ
(
w
,
F
)
,
μ
)
=
a
1
(
Δ
w
,
Δ
μ
)
-
(
L
(
w
,
F
)
,
μ
)
-
(
q
,
μ
).
Equation (17) can be given in the following form:
(21)
(
Φ
(
w
,
F
)
,
μ
)
=
0
,
and it is clear that Φ
(
w
,
F
)
∈
H
-
2
(
Ω
).
Therefore, each point of the minimum of functional (17) on M satisfies (21), and hence it is a weak solution to problems (6) and (4).Therefore, it has been shown that finding a solution to problems (6) and (4) is equivalent to finding a solution to the problem of minimization (13) with the occurrence of constraints {
w
,
F
}
∈
M. The reduced problem can be solved by various methods to find a minimum taking into account the mentioned constarints. Once a solution to the problem of finding an extreme is chosen, various algorithms to solve problems (6) and (4) can be applied.Below, we focus on the method of gradient projection with a restoring constraint [18], which for linear constraints allows for essential simplification of finding a solution to the stated problem.Let us construct an iteration process of minimizingJ
(
w
,
F
) on M using the following scheme:(a)
elementw
0
∈
H
0
2
(
Ω
) is taken arbitrarily;
(b)
after computation ofw
n, F
n
∈
H
0
2
(
Ω
) and w
n
+
1
∈
H
0
2
(
Ω
) is defined successively by solutions to the following problems:
(22)
β
(
w
n
,
F
n
,
μ
)
=
0
,
F
n
∈
H
0
2
(
Ω
)
∀
μ
∈
H
0
2
(
Ω
)
,
(23)
(
Δ
w
n
+
1
,
Δ
μ
)
=
(
Δ
w
n
,
Δ
μ
)
-
ρ
n
(
Φ
(
w
n
,
F
n
)
,
μ
)
hhhhhhhhhhhhhhhhhhhhhhhh
∀
μ
∈
H
0
2
(
Ω
)
;
(c)
coefficientρ
n on step (b) is defined by the condition
(24)
J
(
w
n
+
1
,
F
n
+
1
)
-
J
(
w
n
,
F
n
)
≤
ε
(
Φ
(
w
n
,
F
n
)
,
w
n
+
1
-
w
n
)
,
0
<
ε
<
1
,where ε stands for a parameter of the method.Theorem 3.
For the iteration process (22) to (25) (
Φ
(
w
n
,
F
n
)
,
μ
)
→
0 for n
→
0 an arbitrary initial point {
w
0
,
F
0
}
∈
M, obtained through this procedure series {
w
n
,
F
n
} includes a subseries convergent to the weak solution to the problem ((6) and (4)).Proof.
A possibility of constructing the series{
w
n
,
F
n
} is yielded by an observation that for all ρ
n
w
n
+
1
∈
H
0
2
(
Ω
) and, consequently, L
(
w
n
+
1
,
w
n
+
1
)
∈
H
-
2
(
Ω
), ∇
k
2
w
n
+
1
∈
H
-
2
(
Ω
) [19]. It means that the coupling equation β
(
w
n
+
1
,
F
n
+
1
,
μ
)
=
0 is solvable. Consider the following difference:
(25)
Δ
J
n
=
J
(
w
n
+
1
,
F
n
+
1
)
-
J
(
w
n
,
F
n
)
=
1
2
(
Δ
(
w
n
+
1
-
w
n
)
,
Δ
(
w
n
+
1
+
w
n
)
)
+
1
2
(
Δ
(
F
n
+
1
-
F
n
)
,
Δ
(
F
n
+
1
+
F
n
)
)
-
(
q
,
w
n
+
1
-
w
n
)
.
Owing to{
w
n
,
F
n
}
∈
M, {
w
n
+
1
,
F
n
+
1
}
∈
M, (25) gives
(26)
Δ
J
n
=
(
Φ
(
w
n
,
F
n
)
,
δ
w
)
+
1
2
∥
Δ
δ
w
∥
0
,
Ω
2
+
1
2
∥
Δ
δ
F
∥
0
,
Ω
2
,
where δ
w
=
w
n
+
1
-
w
n, δ
F
=
F
n
+
1
-
F
n. Taking (23) into account, one observes that δ
w serves as a general solution to the boundary value problem:
(27)
Δ
2
δ
w
=
-
ρ
n
Φ
(
w
n
,
F
n
)
,
δ
w
∈
H
0
2
(
Ω
)
.
Further on it means that(28)
δ
w
=
-
ρ
n
G
[
Φ
(
w
n
,
F
n
)
]
,
where G
[
•
]
:
H
-
2
(
Ω
)
→
H
0
2
(
Ω
) stands for the linear bounded operator being inversed to operator Δ
2
(
•
). Therefore
(29)
Δ
J
n
=
-
ρ
n
(
Φ
(
w
n
,
F
n
)
,
G
[
Φ
(
w
n
,
F
n
)
]
)
+
1
2
∥
Δ
δ
w
∥
0
,
Ω
2
+
1
2
∥
Δ
δ
F
∥
0
,
Ω
2
.
Let us proceed to the second order terms. Taking in (11) μ
=
δ
w and applying (28), one gets
(30)
∥
Δ
δ
w
∥
0
,
Ω
2
=
-
ρ
n
(
Φ
(
w
n
,
F
n
)
,
δ
w
)
=
ρ
n
2
(
Φ
(
w
n
,
F
n
)
,
G
[
Φ
(
w
n
,
F
n
)
]
)
.
Let us estimate the last term in (29). Since {
w
n
,
F
n
}
∈
M and {
w
n
+
1
,
F
n
+
1
}
∈
M, then for δ
F the following equation should be satisfied:
(31)
(
Δ
δ
F
,
Δ
μ
)
+
(
L
(
w
n
,
δ
w
)
,
μ
)
+
(
∇
k
2
δ
w
,
μ
)
+
1
2
(
L
(
δ
w
,
δ
w
)
,
μ
)
=
0
δ
F
∈
H
0
2
(
Ω
)
,
∀
μ
∈
H
0
2
(
Ω
)
.
This, in particular, yields(32)
∥
Δ
δ
F
∥
0
,
Ω
≤
c
7
(
∥
∇
k
2
δ
w
∥
L
1
(
Ω
)
∥
L
(
w
n
,
δ
w
)
∥
L
1
(
Ω
)
F
i
f
+
∥
L
(
δ
w
,
δ
w
)
∥
L
1
(
Ω
)
+
∥
∇
k
2
δ
w
∥
L
1
(
Ω
)
)
.
However,w
n belongs to the bounded set in H
0
2
(
Ω
) for arbitrary n. It implies that ∥
Δ
δ
F
∥
0
,
Ω
≤
c
8
∥
Δ
δ
w
∥
0
,
Ω
2 or equivalently
(33)
∥
Δ
δ
F
∥
0
,
Ω
2
≤
c
9
ρ
n
4
(
Φ
(
w
n
,
F
n
)
,
G
[
Φ
(
w
n
,
F
n
)
]
)
2
.
Substituting (30), (33) into (29), and taking into account both positive defiantness (in the sense of (
Φ
(
w
n
,
F
n
)
,
G
[
Φ
(
w
n
,
F
n
)
]
)
≥
α
∥
Φ
(
w
n
,
F
n
)
∥
2) and the constraints of operator G
[
•
], one gets
(34)
Δ
J
n
≤
-
ρ
n
c
10
∥
Φ
(
w
n
,
F
n
)
∥
2
×
(
-
1
+
ρ
n
2
+
c
1
ρ
n
3
2
∥
Φ
(
w
n
,
F
n
)
∥
2
)
.
The latter estimation shows that the valuesρ
n
≠
0 are responsible for the satisfaction of inequality (24). For this purpose ρ
n should be chosen in the following way:
(35)
ρ
n
2
+
c
1
ρ
n
3
2
∥
Φ
(
w
n
,
F
n
)
∥
2
≤
1
-
ε
.
It can always be done, since0
<
ε
<
1.
Takingρ
n in accordance with the algorithm applied so far, the following estimations are obtained on each step:
(36)
Δ
J
n
≤
-
ρ
n
ε
∥
Φ
(
w
n
,
F
n
)
∥
2
,
which means that for the arbitrarily taken n we have J
n
+
1
-
J
n
≤
0. Since functional J is bounded from below, the last inequality yields for n
→
∞
Δ
J
n
→
0. Besides, (36) gives
(37)
∥
Φ
(
w
n
,
F
n
)
∥
2
≤
-
Δ
J
n
ε
ρ
n
.
Let us emphasize that the so far introduced algorithm of the choice ofρ
n guarantees that for arbitrary n we have ρ
n
≥
ρ
0
>
0. In fact, because Δ
J
n
≤
0, then
(38)
J
(
w
n
,
F
n
)
≤
J
(
w
0
,
F
0
)
=
A
.
Owing to (38), norms ∥
w
n
∥
2
,
Ω
,
∥
F
n
∥
2
,
Ω are bounded. Therefore, also the norm ∥
Φ
(
w
n
,
F
n
)
∥ is bounded. In addition, taking (37) into account, we have ∥
Φ
(
w
n
,
F
n
)
∥
→
0 for n
→
∞, and consequently, also (
Φ
(
w
n
,
F
n
)
,
μ
)
→
0 for n
→
∞for all μ
∈
H
0
2
(
Ω
). The occurrence of a convergent subseries follows now from a bound of norms ∥
w
n
∥
2
,
Ω
,
∥
F
n
∥
2
,
Ω (see proof of Theorem 1).We have shown in the above the convergence of the reduction procedure of system (2) to the successive solution to the biharmonic Germain-Lagrange type equation. The applied procedure linearizes and decreases the order of the input equations. We have proposed a further development of this approach on the basis of reduction of the biharmonic equation to that of the Poisson-type. The latter approach allows us to decrease four times the order of system (2).


### 3.2. Iterative Procedure of Reduction of the Germain-Lagrange Equations Type to the Poisson Equations Type

The following original iterative procedure is proposed.We consider a biharmonic equation given in the bounded convex spaceΩ
∈
R
2:
(39)
Δ
2
w
(
x
,
y
)
=
g
(
x
,
y
)
.On the space boundary the following boundary conditions are given:(40)
w
|
Γ
=
0
,
Δ
w
-
χ
∂
w
∂
n
|
Γ
=
0
,
where χ denotes the curvature of boundary Γ.Let us introduce the following new functionM
w
=
Δ
w. Substituting this function into (39) the following system of two Poisson-type equations is obtained:
(41)
Δ
M
w
(
x
,
y
)
=
g
(
x
,
y
)
,
Δ
w
(
x
,
y
)
=
M
w
(
x
,
y
)
,
x
,
y
∈
Ω
.Boundary conditions have the following form:(42)
w
(
x
,
y
)
|
x
,
y
∈
Γ
=
0
,
M
w
(
x
,
y
)
|
x
,
y
∈
Γ
=
0
.Therefore, a solution of the biharmonic equation is divided into a solution of two Poisson-type equations. Below, we prove convergence of the proposed procedure.Let us define the following set for the functionξ
=
{
w
,
M
w
}:
(43)
E
=
{
ξ
∈
ρ
∞
(
Ω
)
∣
ξ
|
Γ
=
0
}
,
where ρ
∞
(
Ω
) is the set of functions infinitely many times differentiable on Ω
¯
∈
R
2. Closure of set (43) in norm H
2
(
Ω
) is a subspace in H
2
(
Ω
) which is denoted by V
(
Ω
). It is clear that V
(
Ω
)
=
H
2
(
Ω
)
∩
H
0
1
(
Ω
).It is known (see [13]) that a solution to problems (39) and (40) is equivalent to minimization on V
(
Ω
) of the following functional:
(44)
J
(
v
)
=
1
2
∫
Ω
|
Δ
ξ
|
2
d
Ω
-
∫
Ω
g
ξ
d
Ω
-
1
2
∫
Γ
χ
|
∂
ξ
∂
n
|
2
d
s
.Hybrid Variation Problem Formulation. We assume that instead of functional (44) the following one is minimized:
(45)
Φ
(
ξ
,
ψ
,
α
)
=
1
2
∫
Ω
|
ψ
|
2
d
Ω
-
∫
Ω
g
ξ
d
Ω
-
1
2
∫
Γ
χ
|
α
|
2
d
s
,
on such triads (
ξ
,
ψ
,
α
)
∈
V
(
Ω
)
×
L
2
(
Ω
)
×
L
2
(
Γ
) that their elements are coupled through equalities -
Δ
ξ
=
ψ
,
(
∂
w
/
∂
n
)
|
Γ
=
α.Let us define the space of the following functions:(46)
P
(
Ω
)
=
{
(
ξ
,
ψ
,
α
)
∈
H
0
1
(
Ω
)
×
L
2
(
Ω
)
×
L
2
(
Γ
)
∣
g
g
∀
μ
∈
H
1
(
Ω
)
,
β
[
(
ξ
,
ψ
,
α
)
,
μ
]
=
0
(
ξ
,
ψ
,
α
)
}
,
where the bilinear form β
[
·
,
·
] is defined as
(47)
β
[
(
ξ
,
ψ
,
α
)
,
μ
]
=
1
2
∫
Ω
∇
ξ
∇
μ
d
Ω
-
∫
Ω
ψ
μ
d
Ω
-
∫
Γ
α
μ
d
s
.Theorem 4.
If the spaceΩ is convex and has a Lipschitz continuous boundary Γ, then, first, the map (
ξ
,
ψ
,
α
)
∈
P
(
Ω
)
→
|
ψ
|
0
,
Ω
2, now and later |
ξ
|
m
,
Ω
=
(
∑
|
k
|
=
m
∫
Ω
|
∂
k
ξ
|
2
d
Ω
)
1
/
2, ∥
ξ
∥
m
,
Ω
=
(
∑
|
k
|
≤
m
∫
Ω
|
∂
k
ξ
|
2
d
Ω
)
1
/
2
) is the norm on space P
(
Ω
) equivalent to the real dot product form (
ξ
,
ψ
,
α
)
∈
P
(
Ω
)
→
(
|
ξ
|
1
,
Ω
2
+
|
ψ
|
0
,
Ω
2
+
|
α
|
0
,
Γ
2
)
1
/
2 transforming P
(
Ω
) into a Hilbert space; second, if (
ξ
,
ψ
,
α
)
∈
P
(
Ω
), then
(48)
(
ξ
,
ψ
,
α
)
∈
V
×
L
2
(
Ω
)
×
L
2
(
Γ
)
,
-
Δ
ξ
=
ψ
,
∂
w
∂
n
|
Γ
=
α
,
and if (48) is satisfied, then (
ξ
,
ψ
,
α
)
∈
P
(
Ω
).Proof.
In the beginning we show the second statement. SinceΩ has a continuous boundary, then the following Green formula holds:
(49)
∫
Ω
∇
ξ
∇
μ
d
Ω
=
-
∫
Ω
Δ
ξ
μ
d
Ω
+
∫
Γ
∂
ξ
∂
n
μ
d
s
,
hhhhhhhhh
∀
ξ
∈
H
2
(
Ω
)
,
∀
μ
∈
H
1
(
Ω
)
.
Let (
ξ
,
ψ
,
α
)
∈
P
(
Ω
). Then ξ
∈
H
0
1
(
Ω
), ψ
∈
L
2
(
Ω
), α
∈
L
2
(
Γ
), and β
[
(
ξ
,
ψ
,
α
)
,
μ
]
=
0 for all μ
∈
H
1
(
Ω
).
The last condition, in particular for∀
μ
∈
H
0
1
(
Ω
) yields
(50)
∫
Ω
∇
ξ
∇
μ
d
Ω
=
∫
Ω
ψ
μ
d
Ω
.
It follows from (50) that v appears as a solution to the Dirichlet problem for the operator -
Δ for ξ
|
Γ
=
0. Since space Ω is convex, therefore ξ
∈
H
2
(
Ω
) ([13], Section 7.1, page 373), and consequently ξ
∈
H
2
(
Ω
)
∩
H
0
1
(
Ω
). Using now (50) for μ
∈
H
0
1
(
Ω
), we get -
Δ
ξ
=
ψ. Using the Green formula for μ
∈
H
1
(
Ω
), we find that (
∂
ξ
/
∂
n
)
|
Γ
=
α. Assume that (48) holds. We show that (
ξ
,
ψ
,
α
)
∈
P
(
Ω
). Since ξ
∈
V
(
Ω
)
⊂
H
2
(
Ω
) and -
Δ
ξ
=
ψ, (
∂
ξ
/
∂
n
)
|
Γ
=
α, then (49) yields ∫
Ω
∇
ξ
∇
μ
d
Ω
=
∫
Ω
ψ
μ
d
Ω
+
∫
Γ
α
μ
d
s for all μ
∈
H
1
(
Ω
); that is, β
[
(
ξ
,
ψ
,
α
)
,
μ
]
=
0 for all μ
∈
H
1
(
Ω
). Besides ξ
∈
V
(
Ω
)
⊂
H
0
1
(
Ω
) and the second statement is proved.
Let us prove the first statement. Endowed with the multiplication normP
(
Ω
) is a Hilbert space. Let (
ξ
,
ψ
,
α
)
∈
P
(
Ω
). Then, as it has been shown, ξ
∈
H
2
(
Ω
)
∩
H
0
1
(
Ω
). For μ
=
ξ condition β
[
(
ξ
,
ψ
,
α
)
,
μ
]
=
0 yields
(51)
|
ξ
|
1
,
Ω
2
≤
C
1
|
ψ
|
0
,
Ω
|
ξ
|
0
,
Ω
.
Let us introduce spaceM
⊂
H
1
(
Ω
) such that H
1
(
Ω
)
=
H
0
1
(
Ω
)
⊕
M. Besides, let us introduce the following operator B
:
H
1
(
Ω
)
→
L
2
(
Ω
) defined in the following way: for all ψ
∈
L
2
(
Ω
)
α
=
B
ψ
∈
L
2
(
Γ
) we have a unique solution to the following equation:
(52)
∫
Γ
α
μ
d
s
=
∫
Ω
∇
ξ
∇
μ
d
Ω
-
∫
Ω
ψ
μ
d
Ω
∀
μ
∈
M
.
Under the condition thatξ
∈
H
0
1
(
Ω
) satisfies the following:
(53)
∫
Ω
∇
ξ
∇
μ
d
Ω
=
∫
Ω
ψ
μ
d
Ω
μ
∈
H
0
1
(
Ω
)
.
It is not difficult to verify that under the theorem conditions,B
ψ
=
-
B
Δ
ξ
=
(
∂
ξ
/
∂
n
)
|
Γ; that is, B stands for the operator of external normal derivative for ξ
∈
H
2
(
Ω
)
∩
H
0
1
(
Ω
). This operator is bounded; that is, ξ
∈
H
2
(
Ω
)
∩
H
0
1
(
Ω
). We denote its norm by ∥
B
∥. Then ∥
B
∥
=
sup
⁡
v
∈
H
2
(
Ω
)
∩
H
0
1
(
Ω
)
(
∥
∂
ξ
/
∂
n
∥
0
,
Γ
/
|
Δ
ξ
|
0
,
Ω
), where ∥
•
∥
0
,
Γ denotes the norm associated with the scalar product (
α
,
β
)
L
2
(
Γ
)
=
∫
Γ
α
β
d
s. Therefore, for all ψ
∈
L
2
(
Ω
)
∥
α
∥
0
,
Γ
≤
∥
B
∥
|
ψ
|
0
,
Ω. Hence, taking (52) into account we get (
∥
ξ
∥
1
,
Ω
+
∥
ψ
∥
0
,
Ω
+
∥
α
∥
0
,
Γ
)
≤
C
2
∥
ψ
∥
0
,
Ω and the theorem is proved.Results of this theorem allow us to transit from minimization of functional (45) on space V
(
Ω
) to minimization of functional (45) on space P
(
Ω
).Theorem 5.
Letw
∈
V
(
Ω
) be a solution to problem (44), then
(54)
Φ
(
w
,
-
Δ
w
,
∂
w
∂
n
)
⟶
inf
⁡
(
v
,
ψ
,
α
)
∈
P
_
(
Ω
)
Φ
(
ξ
,
ψ
,
α
)
.
In this case the triad (
w
,
-
Δ
w
,
∂
w
/
∂
n
)
∈
P
(
Ω
) is the unique solution to the problem of minimization of (54).Proof.
We prove that a symmetric bilinear forma
(
[
(
ξ
,
ψ
,
α
)
,
(
η
,
φ
,
β
)
]
)
=
∫
Ω
ψ
φ
d
Ω
-
∫
Γ
χ
α
β
d
s, (
ξ
,
ψ
,
α
), (
η
,
φ
,
β
)
∈
P
(
Ω
) is continuous and elliptic on P
(
Ω
).
Owing to Theorem4, if (
ξ
,
ψ
,
α
)
,
(
η
,
φ
,
β
)
∈
P
(
Ω
), then -
Δ
ξ
=
ψ; (
∂
ξ
/
∂
n
)
|
Γ
=
α; -
Δ
η
=
φ; (
∂
η
/
∂
n
)
|
Γ
=
β. Then we have
(55)
a
(
[
(
ξ
,
ψ
,
α
)
,
(
η
,
ϕ
,
β
)
]
)
=
∫
Ω
Δ
ξ
Δ
η
d
Ω
-
∫
Γ
∂
ξ
∂
n
∂
η
∂
n
d
s
.
Forη
=
ξ, φ
=
ψ, β
=
α from (55) we obtain ([13], Section 1.2, page 38)
(56)
a
(
[
(
ξ
,
ψ
,
α
)
,
(
η
,
ψ
,
α
)
]
)
=
∫
Ω
[
(
∂
2
ξ
∂
x
2
)
2
+
2
∂
2
ξ
∂
x
2
∂
2
ξ
∂
y
2
+
(
∂
2
ξ
∂
x
2
)
2
+
2
(
∂
2
ξ
∂
x
∂
y
)
2
]
d
Ω
≥
|
ψ
|
0
,
Ω
2
.
Then, elliptic property has been proved byP
.Continuity of the bilinear form is evident. It means that the problem of minimization
(57)
Φ
(
ξ
*
,
ψ
*
,
α
*
)
⟶
inf
⁡
(
u
,
φ
,
β
)
∈
P
(
Ω
)
Φ
(
η
,
φ
,
β
)
has a solution which is unique. Let us find a link between a solution to problem (37) as well as problems (41) and (42). If (
ξ
*
,
ψ
*
,
α
*
)
∈
P
(
Ω
) is a solution to problem (57), then the following relations should hold:
(58)
∫
Ω
ψ
*
φ
d
Ω
-
∫
Ω
η
g
d
Ω
-
∫
Γ
χ
α
*
β
d
s
=
0
hhhhhhhhhhhhhhhh
∀
(
η
,
φ
,
β
)
∈
P
(
Ω
)
.
Since(
ξ
*
,
ψ
*
,
α
*
)
∈
P
(
Ω
), then -
Δ
ξ
*
=
ψ
*, (
∂
ξ
*
/
∂
n
)
|
Γ
=
α
*, and ξ
*
∈
P
(
Ω
). Therefore, taking (58) into account we get ∫
Ω
Δ
ξ
*
Δ
η
d
Ω
-
∫
Γ
χ
(
∂
ξ
*
/
∂
n
)
(
∂
η
/
∂
n
)
d
s
=
∫
Γ
η
g
d
s. Therefore, ξ
* coincides with the solution w to problems (39) and (40), and ψ
*
=
-
Δ
w, α
*
=
(
∂
w
/
∂
n
)
|
Γ.Remark 6.
Since the space is convex and its boundary is regular, then forg
∈
H
-
1
(
Ω
) a solution to problems (39) and (40),
(59)
w
∈
H
3
(
Ω
)
∩
H
0
1
(
Ω
)
,
Δ
w
∈
H
1
(
Ω
)
.Solution to the Minimization Problem (44). We show that a solution to problem (44) can be reduced to a solution of successive Dirichlet problems for the operator -
Δ.For further analysis it is suitable to introduce a linear transformationA
:
L
2
(
Ω
)
→
H
0
1
(
Ω
) in the following way: if ψ
∈
L
2
(
Ω
) is a given function, then the function ξ
=
A
ψ
∈
H
0
1
(
Ω
) is a unique solution to the equation ∫
Ω
∇
ξ
∇
μ
d
Ω
=
∫
Ω
ψ
μ
d
Ω
∀
μ
∈
H
0
1
(
Ω
), for ξ
∈
H
0
1
(
Ω
). This means that space P
(
Ω
), defined by (43), can be presented in the following form:
(60)
P
(
Ω
)
=
{
(
ξ
,
ψ
,
α
)
∈
H
0
1
(
Ω
)
×
L
2
(
Ω
)
×
L
2
(
Γ
)
∣
f
ξ
=
A
ψ
,
α
=
B
ψ
(
ξ
,
ψ
,
α
)
∈
H
0
1
(
Ω
)
×
L
2
(
Ω
)
}
.Problem (44) is equivalent to the following problem of optimal control:
(61)
min
⁡
ψ
∈
L
2
(
Ω
)
[
1
2
∫
Ω
|
ψ
|
2
d
Ω
-
∫
Ω
g
ξ
d
Ω
-
1
2
∫
Γ
χ
|
α
|
2
d
s
]
,
where the state function ξ and α are coupled via control ψ
∈
L
2
(
Ω
) through the following state equations:
(62)
ξ
∈
H
0
1
(
Ω
)
,
ξ
=
A
ψ
,
α
∈
L
2
(
Γ
)
,
α
=
B
ψ
.As it follows from Remark6, although the optimal control ψ is sought on L
2
(
Ω
), its regularity is higher for g
∈
H
-
1
(
Ω
)
ψ
∈
H
1
(
Ω
)
. In this case the following trace is defined: ψ
|
Γ
=
λ
,
λ
∈
M. Furthermore, besides (62), we require that ∫
Ω
∇
ψ
λ
∇
μ
d
Ω
=
∫
Ω
g
μ
d
Ω for all μ
∈
H
0
1
(
Ω
) and ψ
λ
-
λ
∈
H
0
1
(
Ω
). Then, if ξ
λ
=
A
ψ
λ and α
λ
=
B
ψ
λ, (61) implies that min
⁡
ψ
∈
L
2
(
Ω
)
Φ
(
ξ
,
ψ
,
α
)
=
min
⁡
λ
∈
M
D
(
λ
), where
(63)
D
(
λ
)
=
-
1
2
∫
Ω
|
ψ
λ
|
2
d
Ω
-
∫
Γ
λ
α
λ
d
s
-
1
2
∫
Γ
χ
|
α
λ
|
2
d
s
∀
λ
∈
M
.The fundamental idea consists now in the application of a gradient method to the problem of minimization (63).Let us takeM
′ as a dual space for space M, and let 〈
·
,
·
〉 denote the relation of duality between spaces M and M
′. We denote by D
′
∈
M
′ a derivative of the functional D
(
λ
). Let us introduce a map S
:
M
→
H
1
(
Ω
) in the following way: for λ
∈
M
φ
λ
∘
=
S
(
λ
) is a unique function from H
1
(
Ω
) satisfying the condition
(64)
∫
Ω
∇
ϕ
∘
λ
∇
μ
d
Ω
=
0
,
∀
μ
∈
H
0
1
(
Ω
)
,
ϕ
∘
λ
-
λ
∈
H
0
1
(
Ω
)
.Theorem 7.
For an arbitraryλ
∈
M defined in (63), the functional is differentiable and its derivative is defined by the relation
(65)
∀
μ
∈
M
〈
D
′
(
λ
)
,
μ
〉
=
∫
Ω
φ
∘
θ
φ
∘
μ
d
Ω
,
w
h
e
r
e
φ
∘
μ
=
S
(
μ
)
,
φ
∘
θ
=
S
(
θ
λ
)
,
θ
λ
=
λ
+
χ
α
λ
,
θ
λ
∈
M
.Proof.
Differentiating (63) yields
(66)
〈
D
′
(
λ
)
,
μ
〉
=
-
∫
Ω
ψ
λ
ϕ
∘
μ
d
Ω
-
∫
Γ
α
λ
μ
d
Ω
-
∫
Γ
(
λ
+
χ
α
λ
)
β
∘
μ
d
s
,
where φ
∘
μ
=
S
(
μ
)
,
β
∘
μ
=
B
φ
∘
μ.
From (66) and taking (52) into account we get
(67)
〈
D
′
(
λ
)
,
μ
〉
=
-
∫
Ω
∇
v
λ
∇
ϕ
∘
μ
d
Ω
-
∫
Γ
(
λ
+
χ
α
λ
)
β
∘
μ
d
s
.
First term in (67) is equal to zero due to (64). Let us introduce the function φ
∘
θ
=
S
(
θ
λ
), where θ
λ is defined by (66) and let u
∘
μ
=
A
φ
∘
μ. Then, (66) implies 〈
D
′
(
λ
)
,
μ
〉
=
∫
Ω
ϕ
∘
θ
ϕ
∘
μ
d
Ω
-
∫
Γ
∇
ϕ
∘
θ
∇
u
∘
μ
d
Ω.
However, the last term in this equality is equal to zero due to (64), which ends the proof. The gradient method applied to minimization (63) consists now in the determination of a series {
λ
n
}
n
=
0
∞ of functions λ
n
∈
M via the following iteration scheme:
(68)
μ
∈
M
(
λ
n
+
1
-
λ
n
,
μ
)
M
=
-
ρ
n
〈
D
′
(
λ
n
)
,
μ
〉
,
where (
·
,
·
)
M is the arbitrary scalar product in space M, ρ is the positive parameter, and λ
0 is the arbitrary function of M.
Therefore, one iteration (68) corresponds to successive solutions of the following problems: (a)
find for a given functionλ
n
∈
M a unique function ψ
n
∈
H
1
(
Ω
), satisfying the following relations:
(69)
ψ
n
-
λ
n
∈
H
0
1
(
Ω
)
,
(70)
∀
μ
∈
H
0
1
(
Ω
)
∫
Ω
∇
ψ
n
∇
μ
d
Ω
=
∫
Ω
g
μ
d
Ω
;
(b)
find a functionξ
n
∈
H
0
1
(
Ω
), satisfying the following relation:
(71)
∀
μ
∈
H
0
1
(
Ω
)
∫
Ω
∇
ξ
n
∇
μ
d
Ω
=
∫
Ω
ψ
n
μ
d
Ω
;
(c)
find a functionα
n
∈
L
2
(
Γ
), satisfying the following relation:
(72)
∀
μ
∈
M
∫
Γ
α
n
μ
d
Ω
=
∫
Ω
∇
ξ
n
∇
μ
d
Ω
-
∫
Ω
ψ
n
μ
d
Ω
;
(d)
find a functionλ
n
+
1
∈
M, satisfying the following relation:
(73)
∀
μ
∈
M
(
λ
n
+
1
-
λ
n
,
μ
)
M
=
-
ρ
∫
Ω
ϕ
θ
n
∘
ϕ
μ
∘
d
Ω
;where φ
θ
n
∘
=
S
θ
λ
n, θ
λ
n
=
λ
n
-
χ
B
ψ
n, and φ
μ
∘
=
S
μ.We show that by a proper choice of parameterρ
>
0 the iteration process ((69)–(73)) is convergent for arbitrarily taken initial approximation.Let us first define the mapC
:
H
1
(
Ω
)
→
M in the following way: for any function ψ
∈
H
1
(
Ω
) function C
ψ
∈
M is unique satisfying the following condition:
(74)
∀
μ
∈
M
(
C
ψ
,
μ
)
M
=
∫
Ω
ψ
ϕ
μ
∘
d
Ω
.Let us take∥
C
∥
=
sup
⁡
ψ
∈
H
1
(
Ω
)
(
|
C
ψ
|
μ
/
|
ψ
|
0
,
Ω
), where |
·
|
M denotes the norm associated with the scalar product (
·
,
·
)
M. It is clear that this norm exists, since the map μ
∈
M
→
φ
∘
μ
∈
H
1
(
Ω
) is bounded.Theorem 8.
If parameterρ satisfies the following condition:
(75)
0
<
ρ
<
2
∥
C
∥
2
(
1
-
|
χ
|
L
∞
(
Ω
)
∥
S
∥
∥
B
∥
)
2
,
then the iteration process (69)–(73) is convergent in the sense that
(76)
lim
⁡
n
→
∞
ψ
n
=
ψ
в
L
2
(
Ω
)
,
lim
⁡
n
→
∞
ξ
n
=
ξ
в
H
0
1
(
Ω
)
,
lim
⁡
n
→
∞
α
n
=
α
в
L
2
(
Γ
)
.Proof.
It is sufficient to show thatlim
⁡
n
→
∞
ψ
n
=
0 in L
2
(
Ω
) in the particular case when g
=
0. If we use the definition (74) of map C, then the recurrent formula (73) gives
(77)
λ
n
+
1
=
λ
n
-
ρ
C
φ
θ
n
∘
,
where
φ
θ
n
∘
=
S
θ
λ
n
,
θ
λ
n
=
λ
n
+
χ
B
ψ
n
,
and therefore
(78)
|
λ
n
+
1
|
M
2
=
|
λ
n
|
M
2
-
2
ρ
(
C
ϕ
θ
n
∘
,
λ
n
)
M
+
ρ
2
|
C
ϕ
θ
n
∘
|
M
2
.
Consider the term(
C
φ
θ
n
∘
,
λ
n
)
M:
(79)
(
C
φ
θ
n
∘
,
λ
n
)
M
=
∫
Ω
φ
θ
n
∘
ψ
n
d
Ω
=
-
∫
Γ
(
λ
n
+
χ
B
ψ
n
)
B
ψ
n
d
s
=
∫
Ω
|
ψ
n
|
2
d
Ω
-
∫
Γ
χ
|
B
ψ
n
|
2
d
s
≥
|
ψ
n
|
0
,
Ω
2
.
Let us estimate the norm|
φ
θ
n
∘
|:
(80)
|
φ
θ
n
∘
|
0
,
Ω
=
|
S
(
λ
n
+
χ
B
ψ
n
)
|
0
,
Ω
≤
|
ψ
n
|
0
,
Ω
+
∥
S
∥
|
χ
|
L
∞
(
Ω
)
∥
B
∥
|
ψ
n
|
0
,
Ω
≤
(
1
+
|
χ
|
L
∞
(
Ω
)
∥
S
∥
∥
B
∥
)
|
ψ
n
|
0
,
Ω
.
The latter inequalities and (78) imply the following estimation:
(81)
|
λ
n
+
1
|
M
2
-
|
λ
n
|
M
2
≤
-
ρ
[
2
η
-
∥
C
∥
2
(
1
+
|
χ
|
L
∞
(
Ω
)
∥
S
∥
∥
B
∥
)
2
ρ
]
|
ψ
n
|
0
,
Ω
2
.
Hence, in particular, we get(82)
lim
⁡
n
→
∞
|
ψ
|
0
,
Ω
n
=
0
,
if ρ satisfies inequalities (75). Besides, we have
(83)
lim
⁡
n
→
∞
ξ
n
=
lim
⁡
n
→
∞
A
ψ
n
=
0
в
H
0
1
(
Ω
)
,
lim
⁡
n
→
∞
α
n
=
lim
⁡
n
→
∞
B
ψ
n
=
0
в
L
2
(
Γ
)
,
which finishes the proof.Since convergence of the considered method is guaranteed, any choice of subspaceM, satisfying the condition H
1
(
Ω
)
=
H
0
1
(
Ω
)
⊕
M and any choice of scalar product (
·
,
·
)
M on space is allowed. However, the choice influences parameter ρ, as well as the computation time on each iterations. Finally, we point out a few remarks regarding practical computations of λ
n
+
1
∈
M.If the scalar product inM is defined via the following formula:
(84)
(
λ
,
μ
)
M
=
∫
Ω
S
λ
·
S
μ
d
Ω
,
then as λ
n
+
1 any function from M can be taken, assuming that the following condition is satisfied:
(85)
λ
n
+
1
|
Γ
=
(
λ
n
-
ρ
θ
λ
n
)
|
Γ
.Equality (85) can be understood in the sense of trace equality on a boundary. In fact, if (85) is satisfied, then
(86)
(
λ
n
+
1
-
λ
n
,
μ
)
M
=
-
ρ
∫
Ω
S
θ
n
·
S
μ
d
Ω
=
-
ρ
∫
Ω
φ
θ
n
°
φ
μ
∘
d
Ω
,
which means that conditions of Theorem 7 are satisfied.We may choose also the following scalar product:(87)
either
(
λ
,
μ
)
M
=
∫
Ω
∇
λ
∇
μ
d
Ω
.
or
(
λ
,
μ
)
M
=
∫
Ω
λ
μ
d
Ω
.However, in the latter case one needs to compute gradient〈
D
′
(
λ
)
,
μ
〉 on each step, which extends the computational time.Final Remarks. (1) The proof has been carried out for equations in the hybrid form (2). It can be relatively easily extended into equations regarding displacements. (2) Results can be extended on other types of the differential equations, including nonlinear ones, consisting of a biharmonic operator.


### 3.3. Iterative Procedure for the Reduction of the Karman Equation into the Poisson Equation

In the preceding sections we have proved convergence of the iterative procedures for linearization of (2) by reducing the solution of the eighth order system of nonlinear differential equations into that of the solution to a biharmonic equation, as well as the reduction of the biharmonic equation to the Poisson-type equation in the case of a curvilinear boundary using the finite element method (FEM).While considering a space with the rectangular boundary, we may extend the procedure reported in Section3.1 by introduction of new variables into the iterative procedure of solution to the Poisson-type equations without difficulties.In the case of spaces with the curvilinear boundary, the procedure described in Section3.2 can be applied to solve (2) using the iterative procedure, whose convergence has been proved in Section 3.1.For this purpose new variablesM
w
(
x
,
y
) and M
F
(
x
,
y
) are introduced
(88)
M
w
(
x
,
y
)
=
Δ
w
(
x
,
y
)
,
M
F
(
x
,
y
)
=
Δ
F
(
x
,
y
)
.Then each of differential equations (2) is divided into two Poisson-type equations. The iterative procedure of solution of the obtained system of four Poisson-type equations has the following form:
(89)
Δ
M
w
(
k
)
=
q
+
L
(
w
(
k
-
1
)
,
F
(
k
-
1
)
)
,
Δ
w
(
k
)
=
M
w
(
k
)
,
Δ
M
F
(
k
)
=
L
(
w
(
k
)
,
w
(
k
)
)
,
Δ
F
(
k
)
=
M
F
(
k
)
,
{
x
,
y
}
∈
Ω
.Boundary conditions (4) are transformed to the following form:
(90)
w
|
Γ
=
M
w
|
Γ
=
F
|
Γ
=
M
F
|
Γ
=
0
.The given procedure (89) has advantages over procedure (6), while solving each equation since instead of the fourth order equation that of the second order is solved. Because equations are solved by numerical methods (FDM, FEM) and the approximation of the biharmonic operator has high requirements on the approximating functions, then for the Poisson-type equation one may simplify the procedure (89) of finding a solution by choosing simple approximating functions.In the FDM case, an order of algebraic equations system, after a discretization of the biharmonic equation, is higher for the second order equation, and hence higher expectations are required from computer abilities while solving the problem numerically.


## 3.1. Iterative Linearization Procedure and Reduction of the Karman Equation into Germain-Lagrange Equations

We keep the biharmonic operator in each of (2), and we shift nonlinear terms into their right-hand sides. Assuming that functions on the right-hand sides are computed with respect to their previous step and that the equations are solved successively, the following iterative procedure is proposed:
(6)
Δ
2
w
(
k
)
=
L
(
w
(
k
-
1
)
,
F
(
k
-
1
)
)
+
q
,
Δ
2
F
(
k
)
=
-
1
2
L
(
w
(
k
)
,
w
(
k
)
)
,
{
x
,
y
}
∈
Ω
.On the first step of the iterative procedure the following biharmonic equation for a given loadq
(
x
,
y
) is solved:
(7)
Δ
2
w
(
1
)
(
x
,
y
)
=
q
(
x
,
y
)
.The value ofw
(
1
)
(
x
,
y
) is substituted into the right-hand side of equation system (6), and as a result a biharmonic equation for F
(
1
)
(
x
,
y
) with the known right-hand side is obtained. The value of the stress function found so far is substituted to the first system equation. The process of finding solutions is continued to achieve the required accuracy.Let us note that as a result of the application of the iterative procedure, the Germain-Lagrange type system of equations are obtained.Let us prove convergence of the constructed iterative procedure. LetH
2
(
Ω
) refer to a Sobolev space of functions ξ
=
{
w
,
F
} such that
(8)
ξ
∈
L
2
(
Ω
)
,
∂
ξ
∂
x
i
∈
L
2
(
Ω
)
,
∂
2
ξ
∂
x
i
∂
x
j
∈
L
2
(
Ω
)
;
i
,
j
=
1,2
,
where L
2
(
Ω
) denotes the space of functions being summed with a square in Ω.LetH
0
2
(
Ω
) denote the closure of functions from D
(
Ω
) (space of functions of class C
∞ in Ω, having compact carrier in Ω) in norm H
2
(
Ω
):
(9)
H
0
2
(
Ω
)
=
D
(
Ω
)
¯
H
2
(
Ω
)
=
{
ξ
∈
H
2
(
Ω
)
∣
ξ
|
Γ
=
∂
ξ
∂
n
|
Γ
=
0
}
.Since spaceΩ is bounded and its boundary Γ is efficiently regular, then map ξ
→
∥
Δ
ξ
∥
0
,
Ω defines the norm in H
0
2
(
Ω
) being equivalent to the norm generated by spaces H
2
(
Ω
).Assume thatq
∈
H
-
2
(
Ω
) (H
-
2
(
Ω
) denotes a conjugation to H
2
(
Ω
)). It is known [17] that in this case problems (2) and (4) have a solution (it may happen that it shall be nonunique).Novel Variation Formulation of the Problem. Let us denote by (
·
,
·
) a scalar product in L
2
(
Ω
): (
ξ
,
η
)
=
∫
Ω
ξ
η
d
Ω, and by β
(
w
,
F
,
μ
) a three linear form defined on (
H
0
2
(
Ω
)
)
3:
(10)
β
(
w
,
F
,
μ
)
=
(
Δ
F
,
Δ
μ
)
+
1
2
(
L
(
w
,
w
)
,
μ
)
.Let us define the set(11)
M
=
{
w
,
F
∈
H
0
2
(
Ω
)
∣
∀
μ
∈
H
0
2
(
Ω
)
,
β
(
w
,
F
,
μ
)
=
0
}
,
and square function J
(
w
,
F
)
:
M
→
R
(12)
J
(
w
,
F
)
=
1
2
∥
Δ
w
∥
0
,
Ω
2
+
1
2
∥
Δ
F
∥
0
,
Ω
2
-
(
q
,
w
)
.Theorem 1.
The problem of minimizing (12) on set (11) has, at least, one solution.Proof.
Let{
w
n
,
F
n
}
∈
M be the minimizing series; that is, we have
(13)
J
(
w
n
,
F
n
)
⟶
inf
⁡
{
w
,
F
}
∈
M
J
(
w
,
F
)
,
which exists, since J is a square functional.
For arbitraryw
,
F
∈
H
0
2
(
Ω
) the following inequality holds:
(14)
J
(
w
,
F
)
≥
c
1
∥
w
∥
2
,
Ω
2
+
c
2
∥
F
∥
2
,
Ω
2
-
c
3
·
∥
w
∥
2
,
Ω
,
where ∥
·
∥
2
,
Ω denotes the norm in H
2
(
Ω
) and c
i are certain positive constants. Then, (13) yields c
1
∥
w
n
∥
2
,
Ω
2
+
c
2
∥
F
n
∥
2
,
Ω
2
-
c
3
∥
w
n
∥
2
,
Ω
≤
J
(
w
n
,
F
n
)
≤
J
(
w
0
,
F
0
)
=
A, where w
0
,
F
0 are the arbitrary functions (initial approximation).
Then, the following estimation holds:(15)
c
1
(
∥
w
n
∥
2
,
Ω
2
-
c
3
2
c
1
)
2
+
c
2
∥
F
n
∥
2
,
Ω
2
≤
A
+
c
3
2
4
c
1
,
∥
w
n
∥
2
,
Ω
≤
c
4
,
∥
F
n
∥
2
,
Ω
≤
c
5
.
Therefore, series{
w
n
,
F
n
} is bounded in (
H
0
2
(
Ω
)
)
2. Consequently, one may choose a series {
w
k
,
F
k
} that w
k
→
w
~, F
k
→
F
~ is weak in H
0
2
(
Ω
). Since H
0
2
(
Ω
)
→
L
2
(
Ω
) is compact, then w
k
→
w
~, F
k
→
F
~ is strong in L
2
(
Ω
).
We show that interval{
w
~
,
F
~
} of the minimized series belongs to M; that is, β
(
w
~
,
F
~
,
μ
)
=
0, for all μ
∈
H
0
2
(
Ω
).
Since(
L
(
w
k
,
w
k
)
,
μ
)
=
(
L
(
w
k
,
μ
)
,
w
k
)
∀
μ
∈
H
0
2
(
Ω
), L
(
w
k
,
μ
)
→
L
(
w
~
,
μ
) is weak in H
0
2
(
Ω
), and w
k
→
w
~ is strong in L
2
(
Ω
), we get (
L
(
w
k
,
w
k
)
,
μ
)
=
(
L
(
w
~
,
w
~
)
,
μ
) and consequently, β
(
w
~
,
F
~
,
μ
)
=
0 for all μ
∈
H
0
2
(
Ω
). This means that
(16)
{
w
~
,
F
~
}
∈
M
.
However,J
(
w
,
F
) is half-continuous from below in weak topology on (
H
2
(
Ω
)
)
2, and therefore the following inequality holds: lim
⁡
k
→
∞
¯
J
(
w
k
,
F
k
)
≥
J
(
w
~
,
F
~
).
Then (13) and (16) imply that J
(
w
~
,
F
~
)
≤
inf
⁡
(
w
,
F
)
∈
M
J
(
w
,
F
). Therefore, the following equation holds: J
(
w
~
,
F
~
)
=
inf
⁡
(
w
,
F
)
∈
M
J
(
w
,
F
), which means that {
w
~
,
F
~
}
∈
M is a solution to the minimization problem.Let us explain how points of the minimum of functional (12) are linked with solutions to problems (6) and (4). For this purpose a notation of weak solution shall be introduced.A weak solution to problems (6) and (4) is defined by the pair of functions {
w
,
F
}
∈
M, satisfying the following:
(17)
(
Δ
w
,
Δ
μ
)
-
(
L
(
w
,
F
)
,
μ
)
=
(
q
,
μ
)
∀
μ
∈
H
0
2
(
Ω
)
.Theorem 2.
Points of the functional minimum (12) are weak solutions to problems (6) and (4).Proof.
Let{
w
,
F
}
∈
M be one of the functional (12) minimum points. Let us take η
=
w
+
t
δ
w
δ
w
∈
H
0
2
(
Ω
) and let us choose ξ
=
F
+
δ
F
δ
F
∈
H
0
2
(
Ω
) such that {
η
,
ξ
}
∈
M, that is, in the way that β
(
w
,
F
,
μ
)
=
0 for all μ
∈
H
0
2
(
Ω
). Then J
(
w
,
F
)
≤
J
(
η
,
ξ
). This yields
(18)
t
(
Δ
w
,
Δ
δ
w
)
+
(
Δ
F
,
Δ
δ
F
)
-
t
(
q
,
δ
w
)
+
t
2
2
∥
Δ
δ
w
∥
2
,
Ω
2
+
1
2
∥
Δ
δ
F
∥
2
,
Ω
2
≥
0
,
∀
t
∈
R
,
δ
w
∈
H
0
2
(
Ω
)
,
and by taking μ
=
F, condition β
(
η
,
ξ
,
μ
)
=
0 yields
(19)
(
Δ
F
,
Δ
δ
F
)
=
-
t
(
L
(
w
,
δ
w
)
,
F
)
-
t
2
2
(
L
(
δ
w
,
δ
w
)
,
F
)
.
Substituting (19) into (18), dividing the obtained expression by t and going to the limit for t
→
0, the following inequality is obtained:
(20)
(
Δ
w
,
Δ
δ
w
)
-
(
L
(
w
,
F
)
,
δ
w
)
≥
(
q
,
μ
)
.
Substitutingδ
w by -
δ
w in (20), one obtains equality, that is, (17).
Let us use the following notation(
Φ
(
w
,
F
)
,
μ
)
=
a
1
(
Δ
w
,
Δ
μ
)
-
(
L
(
w
,
F
)
,
μ
)
-
(
q
,
μ
).
Equation (17) can be given in the following form:
(21)
(
Φ
(
w
,
F
)
,
μ
)
=
0
,
and it is clear that Φ
(
w
,
F
)
∈
H
-
2
(
Ω
).
Therefore, each point of the minimum of functional (17) on M satisfies (21), and hence it is a weak solution to problems (6) and (4).Therefore, it has been shown that finding a solution to problems (6) and (4) is equivalent to finding a solution to the problem of minimization (13) with the occurrence of constraints {
w
,
F
}
∈
M. The reduced problem can be solved by various methods to find a minimum taking into account the mentioned constarints. Once a solution to the problem of finding an extreme is chosen, various algorithms to solve problems (6) and (4) can be applied.Below, we focus on the method of gradient projection with a restoring constraint [18], which for linear constraints allows for essential simplification of finding a solution to the stated problem.Let us construct an iteration process of minimizingJ
(
w
,
F
) on M using the following scheme:(a)
elementw
0
∈
H
0
2
(
Ω
) is taken arbitrarily;
(b)
after computation ofw
n, F
n
∈
H
0
2
(
Ω
) and w
n
+
1
∈
H
0
2
(
Ω
) is defined successively by solutions to the following problems:
(22)
β
(
w
n
,
F
n
,
μ
)
=
0
,
F
n
∈
H
0
2
(
Ω
)
∀
μ
∈
H
0
2
(
Ω
)
,
(23)
(
Δ
w
n
+
1
,
Δ
μ
)
=
(
Δ
w
n
,
Δ
μ
)
-
ρ
n
(
Φ
(
w
n
,
F
n
)
,
μ
)
hhhhhhhhhhhhhhhhhhhhhhhh
∀
μ
∈
H
0
2
(
Ω
)
;
(c)
coefficientρ
n on step (b) is defined by the condition
(24)
J
(
w
n
+
1
,
F
n
+
1
)
-
J
(
w
n
,
F
n
)
≤
ε
(
Φ
(
w
n
,
F
n
)
,
w
n
+
1
-
w
n
)
,
0
<
ε
<
1
,where ε stands for a parameter of the method.Theorem 3.
For the iteration process (22) to (25) (
Φ
(
w
n
,
F
n
)
,
μ
)
→
0 for n
→
0 an arbitrary initial point {
w
0
,
F
0
}
∈
M, obtained through this procedure series {
w
n
,
F
n
} includes a subseries convergent to the weak solution to the problem ((6) and (4)).Proof.
A possibility of constructing the series{
w
n
,
F
n
} is yielded by an observation that for all ρ
n
w
n
+
1
∈
H
0
2
(
Ω
) and, consequently, L
(
w
n
+
1
,
w
n
+
1
)
∈
H
-
2
(
Ω
), ∇
k
2
w
n
+
1
∈
H
-
2
(
Ω
) [19]. It means that the coupling equation β
(
w
n
+
1
,
F
n
+
1
,
μ
)
=
0 is solvable. Consider the following difference:
(25)
Δ
J
n
=
J
(
w
n
+
1
,
F
n
+
1
)
-
J
(
w
n
,
F
n
)
=
1
2
(
Δ
(
w
n
+
1
-
w
n
)
,
Δ
(
w
n
+
1
+
w
n
)
)
+
1
2
(
Δ
(
F
n
+
1
-
F
n
)
,
Δ
(
F
n
+
1
+
F
n
)
)
-
(
q
,
w
n
+
1
-
w
n
)
.
Owing to{
w
n
,
F
n
}
∈
M, {
w
n
+
1
,
F
n
+
1
}
∈
M, (25) gives
(26)
Δ
J
n
=
(
Φ
(
w
n
,
F
n
)
,
δ
w
)
+
1
2
∥
Δ
δ
w
∥
0
,
Ω
2
+
1
2
∥
Δ
δ
F
∥
0
,
Ω
2
,
where δ
w
=
w
n
+
1
-
w
n, δ
F
=
F
n
+
1
-
F
n. Taking (23) into account, one observes that δ
w serves as a general solution to the boundary value problem:
(27)
Δ
2
δ
w
=
-
ρ
n
Φ
(
w
n
,
F
n
)
,
δ
w
∈
H
0
2
(
Ω
)
.
Further on it means that(28)
δ
w
=
-
ρ
n
G
[
Φ
(
w
n
,
F
n
)
]
,
where G
[
•
]
:
H
-
2
(
Ω
)
→
H
0
2
(
Ω
) stands for the linear bounded operator being inversed to operator Δ
2
(
•
). Therefore
(29)
Δ
J
n
=
-
ρ
n
(
Φ
(
w
n
,
F
n
)
,
G
[
Φ
(
w
n
,
F
n
)
]
)
+
1
2
∥
Δ
δ
w
∥
0
,
Ω
2
+
1
2
∥
Δ
δ
F
∥
0
,
Ω
2
.
Let us proceed to the second order terms. Taking in (11) μ
=
δ
w and applying (28), one gets
(30)
∥
Δ
δ
w
∥
0
,
Ω
2
=
-
ρ
n
(
Φ
(
w
n
,
F
n
)
,
δ
w
)
=
ρ
n
2
(
Φ
(
w
n
,
F
n
)
,
G
[
Φ
(
w
n
,
F
n
)
]
)
.
Let us estimate the last term in (29). Since {
w
n
,
F
n
}
∈
M and {
w
n
+
1
,
F
n
+
1
}
∈
M, then for δ
F the following equation should be satisfied:
(31)
(
Δ
δ
F
,
Δ
μ
)
+
(
L
(
w
n
,
δ
w
)
,
μ
)
+
(
∇
k
2
δ
w
,
μ
)
+
1
2
(
L
(
δ
w
,
δ
w
)
,
μ
)
=
0
δ
F
∈
H
0
2
(
Ω
)
,
∀
μ
∈
H
0
2
(
Ω
)
.
This, in particular, yields(32)
∥
Δ
δ
F
∥
0
,
Ω
≤
c
7
(
∥
∇
k
2
δ
w
∥
L
1
(
Ω
)
∥
L
(
w
n
,
δ
w
)
∥
L
1
(
Ω
)
F
i
f
+
∥
L
(
δ
w
,
δ
w
)
∥
L
1
(
Ω
)
+
∥
∇
k
2
δ
w
∥
L
1
(
Ω
)
)
.
However,w
n belongs to the bounded set in H
0
2
(
Ω
) for arbitrary n. It implies that ∥
Δ
δ
F
∥
0
,
Ω
≤
c
8
∥
Δ
δ
w
∥
0
,
Ω
2 or equivalently
(33)
∥
Δ
δ
F
∥
0
,
Ω
2
≤
c
9
ρ
n
4
(
Φ
(
w
n
,
F
n
)
,
G
[
Φ
(
w
n
,
F
n
)
]
)
2
.
Substituting (30), (33) into (29), and taking into account both positive defiantness (in the sense of (
Φ
(
w
n
,
F
n
)
,
G
[
Φ
(
w
n
,
F
n
)
]
)
≥
α
∥
Φ
(
w
n
,
F
n
)
∥
2) and the constraints of operator G
[
•
], one gets
(34)
Δ
J
n
≤
-
ρ
n
c
10
∥
Φ
(
w
n
,
F
n
)
∥
2
×
(
-
1
+
ρ
n
2
+
c
1
ρ
n
3
2
∥
Φ
(
w
n
,
F
n
)
∥
2
)
.
The latter estimation shows that the valuesρ
n
≠
0 are responsible for the satisfaction of inequality (24). For this purpose ρ
n should be chosen in the following way:
(35)
ρ
n
2
+
c
1
ρ
n
3
2
∥
Φ
(
w
n
,
F
n
)
∥
2
≤
1
-
ε
.
It can always be done, since0
<
ε
<
1.
Takingρ
n in accordance with the algorithm applied so far, the following estimations are obtained on each step:
(36)
Δ
J
n
≤
-
ρ
n
ε
∥
Φ
(
w
n
,
F
n
)
∥
2
,
which means that for the arbitrarily taken n we have J
n
+
1
-
J
n
≤
0. Since functional J is bounded from below, the last inequality yields for n
→
∞
Δ
J
n
→
0. Besides, (36) gives
(37)
∥
Φ
(
w
n
,
F
n
)
∥
2
≤
-
Δ
J
n
ε
ρ
n
.
Let us emphasize that the so far introduced algorithm of the choice ofρ
n guarantees that for arbitrary n we have ρ
n
≥
ρ
0
>
0. In fact, because Δ
J
n
≤
0, then
(38)
J
(
w
n
,
F
n
)
≤
J
(
w
0
,
F
0
)
=
A
.
Owing to (38), norms ∥
w
n
∥
2
,
Ω
,
∥
F
n
∥
2
,
Ω are bounded. Therefore, also the norm ∥
Φ
(
w
n
,
F
n
)
∥ is bounded. In addition, taking (37) into account, we have ∥
Φ
(
w
n
,
F
n
)
∥
→
0 for n
→
∞, and consequently, also (
Φ
(
w
n
,
F
n
)
,
μ
)
→
0 for n
→
∞for all μ
∈
H
0
2
(
Ω
). The occurrence of a convergent subseries follows now from a bound of norms ∥
w
n
∥
2
,
Ω
,
∥
F
n
∥
2
,
Ω (see proof of Theorem 1).We have shown in the above the convergence of the reduction procedure of system (2) to the successive solution to the biharmonic Germain-Lagrange type equation. The applied procedure linearizes and decreases the order of the input equations. We have proposed a further development of this approach on the basis of reduction of the biharmonic equation to that of the Poisson-type. The latter approach allows us to decrease four times the order of system (2).


## 3.2. Iterative Procedure of Reduction of the Germain-Lagrange Equations Type to the Poisson Equations Type

The following original iterative procedure is proposed.We consider a biharmonic equation given in the bounded convex spaceΩ
∈
R
2:
(39)
Δ
2
w
(
x
,
y
)
=
g
(
x
,
y
)
.On the space boundary the following boundary conditions are given:(40)
w
|
Γ
=
0
,
Δ
w
-
χ
∂
w
∂
n
|
Γ
=
0
,
where χ denotes the curvature of boundary Γ.Let us introduce the following new functionM
w
=
Δ
w. Substituting this function into (39) the following system of two Poisson-type equations is obtained:
(41)
Δ
M
w
(
x
,
y
)
=
g
(
x
,
y
)
,
Δ
w
(
x
,
y
)
=
M
w
(
x
,
y
)
,
x
,
y
∈
Ω
.Boundary conditions have the following form:(42)
w
(
x
,
y
)
|
x
,
y
∈
Γ
=
0
,
M
w
(
x
,
y
)
|
x
,
y
∈
Γ
=
0
.Therefore, a solution of the biharmonic equation is divided into a solution of two Poisson-type equations. Below, we prove convergence of the proposed procedure.Let us define the following set for the functionξ
=
{
w
,
M
w
}:
(43)
E
=
{
ξ
∈
ρ
∞
(
Ω
)
∣
ξ
|
Γ
=
0
}
,
where ρ
∞
(
Ω
) is the set of functions infinitely many times differentiable on Ω
¯
∈
R
2. Closure of set (43) in norm H
2
(
Ω
) is a subspace in H
2
(
Ω
) which is denoted by V
(
Ω
). It is clear that V
(
Ω
)
=
H
2
(
Ω
)
∩
H
0
1
(
Ω
).It is known (see [13]) that a solution to problems (39) and (40) is equivalent to minimization on V
(
Ω
) of the following functional:
(44)
J
(
v
)
=
1
2
∫
Ω
|
Δ
ξ
|
2
d
Ω
-
∫
Ω
g
ξ
d
Ω
-
1
2
∫
Γ
χ
|
∂
ξ
∂
n
|
2
d
s
.Hybrid Variation Problem Formulation. We assume that instead of functional (44) the following one is minimized:
(45)
Φ
(
ξ
,
ψ
,
α
)
=
1
2
∫
Ω
|
ψ
|
2
d
Ω
-
∫
Ω
g
ξ
d
Ω
-
1
2
∫
Γ
χ
|
α
|
2
d
s
,
on such triads (
ξ
,
ψ
,
α
)
∈
V
(
Ω
)
×
L
2
(
Ω
)
×
L
2
(
Γ
) that their elements are coupled through equalities -
Δ
ξ
=
ψ
,
(
∂
w
/
∂
n
)
|
Γ
=
α.Let us define the space of the following functions:(46)
P
(
Ω
)
=
{
(
ξ
,
ψ
,
α
)
∈
H
0
1
(
Ω
)
×
L
2
(
Ω
)
×
L
2
(
Γ
)
∣
g
g
∀
μ
∈
H
1
(
Ω
)
,
β
[
(
ξ
,
ψ
,
α
)
,
μ
]
=
0
(
ξ
,
ψ
,
α
)
}
,
where the bilinear form β
[
·
,
·
] is defined as
(47)
β
[
(
ξ
,
ψ
,
α
)
,
μ
]
=
1
2
∫
Ω
∇
ξ
∇
μ
d
Ω
-
∫
Ω
ψ
μ
d
Ω
-
∫
Γ
α
μ
d
s
.Theorem 4.
If the spaceΩ is convex and has a Lipschitz continuous boundary Γ, then, first, the map (
ξ
,
ψ
,
α
)
∈
P
(
Ω
)
→
|
ψ
|
0
,
Ω
2, now and later |
ξ
|
m
,
Ω
=
(
∑
|
k
|
=
m
∫
Ω
|
∂
k
ξ
|
2
d
Ω
)
1
/
2, ∥
ξ
∥
m
,
Ω
=
(
∑
|
k
|
≤
m
∫
Ω
|
∂
k
ξ
|
2
d
Ω
)
1
/
2
) is the norm on space P
(
Ω
) equivalent to the real dot product form (
ξ
,
ψ
,
α
)
∈
P
(
Ω
)
→
(
|
ξ
|
1
,
Ω
2
+
|
ψ
|
0
,
Ω
2
+
|
α
|
0
,
Γ
2
)
1
/
2 transforming P
(
Ω
) into a Hilbert space; second, if (
ξ
,
ψ
,
α
)
∈
P
(
Ω
), then
(48)
(
ξ
,
ψ
,
α
)
∈
V
×
L
2
(
Ω
)
×
L
2
(
Γ
)
,
-
Δ
ξ
=
ψ
,
∂
w
∂
n
|
Γ
=
α
,
and if (48) is satisfied, then (
ξ
,
ψ
,
α
)
∈
P
(
Ω
).Proof.
In the beginning we show the second statement. SinceΩ has a continuous boundary, then the following Green formula holds:
(49)
∫
Ω
∇
ξ
∇
μ
d
Ω
=
-
∫
Ω
Δ
ξ
μ
d
Ω
+
∫
Γ
∂
ξ
∂
n
μ
d
s
,
hhhhhhhhh
∀
ξ
∈
H
2
(
Ω
)
,
∀
μ
∈
H
1
(
Ω
)
.
Let (
ξ
,
ψ
,
α
)
∈
P
(
Ω
). Then ξ
∈
H
0
1
(
Ω
), ψ
∈
L
2
(
Ω
), α
∈
L
2
(
Γ
), and β
[
(
ξ
,
ψ
,
α
)
,
μ
]
=
0 for all μ
∈
H
1
(
Ω
).
The last condition, in particular for∀
μ
∈
H
0
1
(
Ω
) yields
(50)
∫
Ω
∇
ξ
∇
μ
d
Ω
=
∫
Ω
ψ
μ
d
Ω
.
It follows from (50) that v appears as a solution to the Dirichlet problem for the operator -
Δ for ξ
|
Γ
=
0. Since space Ω is convex, therefore ξ
∈
H
2
(
Ω
) ([13], Section 7.1, page 373), and consequently ξ
∈
H
2
(
Ω
)
∩
H
0
1
(
Ω
). Using now (50) for μ
∈
H
0
1
(
Ω
), we get -
Δ
ξ
=
ψ. Using the Green formula for μ
∈
H
1
(
Ω
), we find that (
∂
ξ
/
∂
n
)
|
Γ
=
α. Assume that (48) holds. We show that (
ξ
,
ψ
,
α
)
∈
P
(
Ω
). Since ξ
∈
V
(
Ω
)
⊂
H
2
(
Ω
) and -
Δ
ξ
=
ψ, (
∂
ξ
/
∂
n
)
|
Γ
=
α, then (49) yields ∫
Ω
∇
ξ
∇
μ
d
Ω
=
∫
Ω
ψ
μ
d
Ω
+
∫
Γ
α
μ
d
s for all μ
∈
H
1
(
Ω
); that is, β
[
(
ξ
,
ψ
,
α
)
,
μ
]
=
0 for all μ
∈
H
1
(
Ω
). Besides ξ
∈
V
(
Ω
)
⊂
H
0
1
(
Ω
) and the second statement is proved.
Let us prove the first statement. Endowed with the multiplication normP
(
Ω
) is a Hilbert space. Let (
ξ
,
ψ
,
α
)
∈
P
(
Ω
). Then, as it has been shown, ξ
∈
H
2
(
Ω
)
∩
H
0
1
(
Ω
). For μ
=
ξ condition β
[
(
ξ
,
ψ
,
α
)
,
μ
]
=
0 yields
(51)
|
ξ
|
1
,
Ω
2
≤
C
1
|
ψ
|
0
,
Ω
|
ξ
|
0
,
Ω
.
Let us introduce spaceM
⊂
H
1
(
Ω
) such that H
1
(
Ω
)
=
H
0
1
(
Ω
)
⊕
M. Besides, let us introduce the following operator B
:
H
1
(
Ω
)
→
L
2
(
Ω
) defined in the following way: for all ψ
∈
L
2
(
Ω
)
α
=
B
ψ
∈
L
2
(
Γ
) we have a unique solution to the following equation:
(52)
∫
Γ
α
μ
d
s
=
∫
Ω
∇
ξ
∇
μ
d
Ω
-
∫
Ω
ψ
μ
d
Ω
∀
μ
∈
M
.
Under the condition thatξ
∈
H
0
1
(
Ω
) satisfies the following:
(53)
∫
Ω
∇
ξ
∇
μ
d
Ω
=
∫
Ω
ψ
μ
d
Ω
μ
∈
H
0
1
(
Ω
)
.
It is not difficult to verify that under the theorem conditions,B
ψ
=
-
B
Δ
ξ
=
(
∂
ξ
/
∂
n
)
|
Γ; that is, B stands for the operator of external normal derivative for ξ
∈
H
2
(
Ω
)
∩
H
0
1
(
Ω
). This operator is bounded; that is, ξ
∈
H
2
(
Ω
)
∩
H
0
1
(
Ω
). We denote its norm by ∥
B
∥. Then ∥
B
∥
=
sup
⁡
v
∈
H
2
(
Ω
)
∩
H
0
1
(
Ω
)
(
∥
∂
ξ
/
∂
n
∥
0
,
Γ
/
|
Δ
ξ
|
0
,
Ω
), where ∥
•
∥
0
,
Γ denotes the norm associated with the scalar product (
α
,
β
)
L
2
(
Γ
)
=
∫
Γ
α
β
d
s. Therefore, for all ψ
∈
L
2
(
Ω
)
∥
α
∥
0
,
Γ
≤
∥
B
∥
|
ψ
|
0
,
Ω. Hence, taking (52) into account we get (
∥
ξ
∥
1
,
Ω
+
∥
ψ
∥
0
,
Ω
+
∥
α
∥
0
,
Γ
)
≤
C
2
∥
ψ
∥
0
,
Ω and the theorem is proved.Results of this theorem allow us to transit from minimization of functional (45) on space V
(
Ω
) to minimization of functional (45) on space P
(
Ω
).Theorem 5.
Letw
∈
V
(
Ω
) be a solution to problem (44), then
(54)
Φ
(
w
,
-
Δ
w
,
∂
w
∂
n
)
⟶
inf
⁡
(
v
,
ψ
,
α
)
∈
P
_
(
Ω
)
Φ
(
ξ
,
ψ
,
α
)
.
In this case the triad (
w
,
-
Δ
w
,
∂
w
/
∂
n
)
∈
P
(
Ω
) is the unique solution to the problem of minimization of (54).Proof.
We prove that a symmetric bilinear forma
(
[
(
ξ
,
ψ
,
α
)
,
(
η
,
φ
,
β
)
]
)
=
∫
Ω
ψ
φ
d
Ω
-
∫
Γ
χ
α
β
d
s, (
ξ
,
ψ
,
α
), (
η
,
φ
,
β
)
∈
P
(
Ω
) is continuous and elliptic on P
(
Ω
).
Owing to Theorem4, if (
ξ
,
ψ
,
α
)
,
(
η
,
φ
,
β
)
∈
P
(
Ω
), then -
Δ
ξ
=
ψ; (
∂
ξ
/
∂
n
)
|
Γ
=
α; -
Δ
η
=
φ; (
∂
η
/
∂
n
)
|
Γ
=
β. Then we have
(55)
a
(
[
(
ξ
,
ψ
,
α
)
,
(
η
,
ϕ
,
β
)
]
)
=
∫
Ω
Δ
ξ
Δ
η
d
Ω
-
∫
Γ
∂
ξ
∂
n
∂
η
∂
n
d
s
.
Forη
=
ξ, φ
=
ψ, β
=
α from (55) we obtain ([13], Section 1.2, page 38)
(56)
a
(
[
(
ξ
,
ψ
,
α
)
,
(
η
,
ψ
,
α
)
]
)
=
∫
Ω
[
(
∂
2
ξ
∂
x
2
)
2
+
2
∂
2
ξ
∂
x
2
∂
2
ξ
∂
y
2
+
(
∂
2
ξ
∂
x
2
)
2
+
2
(
∂
2
ξ
∂
x
∂
y
)
2
]
d
Ω
≥
|
ψ
|
0
,
Ω
2
.
Then, elliptic property has been proved byP
.Continuity of the bilinear form is evident. It means that the problem of minimization
(57)
Φ
(
ξ
*
,
ψ
*
,
α
*
)
⟶
inf
⁡
(
u
,
φ
,
β
)
∈
P
(
Ω
)
Φ
(
η
,
φ
,
β
)
has a solution which is unique. Let us find a link between a solution to problem (37) as well as problems (41) and (42). If (
ξ
*
,
ψ
*
,
α
*
)
∈
P
(
Ω
) is a solution to problem (57), then the following relations should hold:
(58)
∫
Ω
ψ
*
φ
d
Ω
-
∫
Ω
η
g
d
Ω
-
∫
Γ
χ
α
*
β
d
s
=
0
hhhhhhhhhhhhhhhh
∀
(
η
,
φ
,
β
)
∈
P
(
Ω
)
.
Since(
ξ
*
,
ψ
*
,
α
*
)
∈
P
(
Ω
), then -
Δ
ξ
*
=
ψ
*, (
∂
ξ
*
/
∂
n
)
|
Γ
=
α
*, and ξ
*
∈
P
(
Ω
). Therefore, taking (58) into account we get ∫
Ω
Δ
ξ
*
Δ
η
d
Ω
-
∫
Γ
χ
(
∂
ξ
*
/
∂
n
)
(
∂
η
/
∂
n
)
d
s
=
∫
Γ
η
g
d
s. Therefore, ξ
* coincides with the solution w to problems (39) and (40), and ψ
*
=
-
Δ
w, α
*
=
(
∂
w
/
∂
n
)
|
Γ.Remark 6.
Since the space is convex and its boundary is regular, then forg
∈
H
-
1
(
Ω
) a solution to problems (39) and (40),
(59)
w
∈
H
3
(
Ω
)
∩
H
0
1
(
Ω
)
,
Δ
w
∈
H
1
(
Ω
)
.Solution to the Minimization Problem (44). We show that a solution to problem (44) can be reduced to a solution of successive Dirichlet problems for the operator -
Δ.For further analysis it is suitable to introduce a linear transformationA
:
L
2
(
Ω
)
→
H
0
1
(
Ω
) in the following way: if ψ
∈
L
2
(
Ω
) is a given function, then the function ξ
=
A
ψ
∈
H
0
1
(
Ω
) is a unique solution to the equation ∫
Ω
∇
ξ
∇
μ
d
Ω
=
∫
Ω
ψ
μ
d
Ω
∀
μ
∈
H
0
1
(
Ω
), for ξ
∈
H
0
1
(
Ω
). This means that space P
(
Ω
), defined by (43), can be presented in the following form:
(60)
P
(
Ω
)
=
{
(
ξ
,
ψ
,
α
)
∈
H
0
1
(
Ω
)
×
L
2
(
Ω
)
×
L
2
(
Γ
)
∣
f
ξ
=
A
ψ
,
α
=
B
ψ
(
ξ
,
ψ
,
α
)
∈
H
0
1
(
Ω
)
×
L
2
(
Ω
)
}
.Problem (44) is equivalent to the following problem of optimal control:
(61)
min
⁡
ψ
∈
L
2
(
Ω
)
[
1
2
∫
Ω
|
ψ
|
2
d
Ω
-
∫
Ω
g
ξ
d
Ω
-
1
2
∫
Γ
χ
|
α
|
2
d
s
]
,
where the state function ξ and α are coupled via control ψ
∈
L
2
(
Ω
) through the following state equations:
(62)
ξ
∈
H
0
1
(
Ω
)
,
ξ
=
A
ψ
,
α
∈
L
2
(
Γ
)
,
α
=
B
ψ
.As it follows from Remark6, although the optimal control ψ is sought on L
2
(
Ω
), its regularity is higher for g
∈
H
-
1
(
Ω
)
ψ
∈
H
1
(
Ω
)
. In this case the following trace is defined: ψ
|
Γ
=
λ
,
λ
∈
M. Furthermore, besides (62), we require that ∫
Ω
∇
ψ
λ
∇
μ
d
Ω
=
∫
Ω
g
μ
d
Ω for all μ
∈
H
0
1
(
Ω
) and ψ
λ
-
λ
∈
H
0
1
(
Ω
). Then, if ξ
λ
=
A
ψ
λ and α
λ
=
B
ψ
λ, (61) implies that min
⁡
ψ
∈
L
2
(
Ω
)
Φ
(
ξ
,
ψ
,
α
)
=
min
⁡
λ
∈
M
D
(
λ
), where
(63)
D
(
λ
)
=
-
1
2
∫
Ω
|
ψ
λ
|
2
d
Ω
-
∫
Γ
λ
α
λ
d
s
-
1
2
∫
Γ
χ
|
α
λ
|
2
d
s
∀
λ
∈
M
.The fundamental idea consists now in the application of a gradient method to the problem of minimization (63).Let us takeM
′ as a dual space for space M, and let 〈
·
,
·
〉 denote the relation of duality between spaces M and M
′. We denote by D
′
∈
M
′ a derivative of the functional D
(
λ
). Let us introduce a map S
:
M
→
H
1
(
Ω
) in the following way: for λ
∈
M
φ
λ
∘
=
S
(
λ
) is a unique function from H
1
(
Ω
) satisfying the condition
(64)
∫
Ω
∇
ϕ
∘
λ
∇
μ
d
Ω
=
0
,
∀
μ
∈
H
0
1
(
Ω
)
,
ϕ
∘
λ
-
λ
∈
H
0
1
(
Ω
)
.Theorem 7.
For an arbitraryλ
∈
M defined in (63), the functional is differentiable and its derivative is defined by the relation
(65)
∀
μ
∈
M
〈
D
′
(
λ
)
,
μ
〉
=
∫
Ω
φ
∘
θ
φ
∘
μ
d
Ω
,
w
h
e
r
e
φ
∘
μ
=
S
(
μ
)
,
φ
∘
θ
=
S
(
θ
λ
)
,
θ
λ
=
λ
+
χ
α
λ
,
θ
λ
∈
M
.Proof.
Differentiating (63) yields
(66)
〈
D
′
(
λ
)
,
μ
〉
=
-
∫
Ω
ψ
λ
ϕ
∘
μ
d
Ω
-
∫
Γ
α
λ
μ
d
Ω
-
∫
Γ
(
λ
+
χ
α
λ
)
β
∘
μ
d
s
,
where φ
∘
μ
=
S
(
μ
)
,
β
∘
μ
=
B
φ
∘
μ.
From (66) and taking (52) into account we get
(67)
〈
D
′
(
λ
)
,
μ
〉
=
-
∫
Ω
∇
v
λ
∇
ϕ
∘
μ
d
Ω
-
∫
Γ
(
λ
+
χ
α
λ
)
β
∘
μ
d
s
.
First term in (67) is equal to zero due to (64). Let us introduce the function φ
∘
θ
=
S
(
θ
λ
), where θ
λ is defined by (66) and let u
∘
μ
=
A
φ
∘
μ. Then, (66) implies 〈
D
′
(
λ
)
,
μ
〉
=
∫
Ω
ϕ
∘
θ
ϕ
∘
μ
d
Ω
-
∫
Γ
∇
ϕ
∘
θ
∇
u
∘
μ
d
Ω.
However, the last term in this equality is equal to zero due to (64), which ends the proof. The gradient method applied to minimization (63) consists now in the determination of a series {
λ
n
}
n
=
0
∞ of functions λ
n
∈
M via the following iteration scheme:
(68)
μ
∈
M
(
λ
n
+
1
-
λ
n
,
μ
)
M
=
-
ρ
n
〈
D
′
(
λ
n
)
,
μ
〉
,
where (
·
,
·
)
M is the arbitrary scalar product in space M, ρ is the positive parameter, and λ
0 is the arbitrary function of M.
Therefore, one iteration (68) corresponds to successive solutions of the following problems: (a)
find for a given functionλ
n
∈
M a unique function ψ
n
∈
H
1
(
Ω
), satisfying the following relations:
(69)
ψ
n
-
λ
n
∈
H
0
1
(
Ω
)
,
(70)
∀
μ
∈
H
0
1
(
Ω
)
∫
Ω
∇
ψ
n
∇
μ
d
Ω
=
∫
Ω
g
μ
d
Ω
;
(b)
find a functionξ
n
∈
H
0
1
(
Ω
), satisfying the following relation:
(71)
∀
μ
∈
H
0
1
(
Ω
)
∫
Ω
∇
ξ
n
∇
μ
d
Ω
=
∫
Ω
ψ
n
μ
d
Ω
;
(c)
find a functionα
n
∈
L
2
(
Γ
), satisfying the following relation:
(72)
∀
μ
∈
M
∫
Γ
α
n
μ
d
Ω
=
∫
Ω
∇
ξ
n
∇
μ
d
Ω
-
∫
Ω
ψ
n
μ
d
Ω
;
(d)
find a functionλ
n
+
1
∈
M, satisfying the following relation:
(73)
∀
μ
∈
M
(
λ
n
+
1
-
λ
n
,
μ
)
M
=
-
ρ
∫
Ω
ϕ
θ
n
∘
ϕ
μ
∘
d
Ω
;where φ
θ
n
∘
=
S
θ
λ
n, θ
λ
n
=
λ
n
-
χ
B
ψ
n, and φ
μ
∘
=
S
μ.We show that by a proper choice of parameterρ
>
0 the iteration process ((69)–(73)) is convergent for arbitrarily taken initial approximation.Let us first define the mapC
:
H
1
(
Ω
)
→
M in the following way: for any function ψ
∈
H
1
(
Ω
) function C
ψ
∈
M is unique satisfying the following condition:
(74)
∀
μ
∈
M
(
C
ψ
,
μ
)
M
=
∫
Ω
ψ
ϕ
μ
∘
d
Ω
.Let us take∥
C
∥
=
sup
⁡
ψ
∈
H
1
(
Ω
)
(
|
C
ψ
|
μ
/
|
ψ
|
0
,
Ω
), where |
·
|
M denotes the norm associated with the scalar product (
·
,
·
)
M. It is clear that this norm exists, since the map μ
∈
M
→
φ
∘
μ
∈
H
1
(
Ω
) is bounded.Theorem 8.
If parameterρ satisfies the following condition:
(75)
0
<
ρ
<
2
∥
C
∥
2
(
1
-
|
χ
|
L
∞
(
Ω
)
∥
S
∥
∥
B
∥
)
2
,
then the iteration process (69)–(73) is convergent in the sense that
(76)
lim
⁡
n
→
∞
ψ
n
=
ψ
в
L
2
(
Ω
)
,
lim
⁡
n
→
∞
ξ
n
=
ξ
в
H
0
1
(
Ω
)
,
lim
⁡
n
→
∞
α
n
=
α
в
L
2
(
Γ
)
.Proof.
It is sufficient to show thatlim
⁡
n
→
∞
ψ
n
=
0 in L
2
(
Ω
) in the particular case when g
=
0. If we use the definition (74) of map C, then the recurrent formula (73) gives
(77)
λ
n
+
1
=
λ
n
-
ρ
C
φ
θ
n
∘
,
where
φ
θ
n
∘
=
S
θ
λ
n
,
θ
λ
n
=
λ
n
+
χ
B
ψ
n
,
and therefore
(78)
|
λ
n
+
1
|
M
2
=
|
λ
n
|
M
2
-
2
ρ
(
C
ϕ
θ
n
∘
,
λ
n
)
M
+
ρ
2
|
C
ϕ
θ
n
∘
|
M
2
.
Consider the term(
C
φ
θ
n
∘
,
λ
n
)
M:
(79)
(
C
φ
θ
n
∘
,
λ
n
)
M
=
∫
Ω
φ
θ
n
∘
ψ
n
d
Ω
=
-
∫
Γ
(
λ
n
+
χ
B
ψ
n
)
B
ψ
n
d
s
=
∫
Ω
|
ψ
n
|
2
d
Ω
-
∫
Γ
χ
|
B
ψ
n
|
2
d
s
≥
|
ψ
n
|
0
,
Ω
2
.
Let us estimate the norm|
φ
θ
n
∘
|:
(80)
|
φ
θ
n
∘
|
0
,
Ω
=
|
S
(
λ
n
+
χ
B
ψ
n
)
|
0
,
Ω
≤
|
ψ
n
|
0
,
Ω
+
∥
S
∥
|
χ
|
L
∞
(
Ω
)
∥
B
∥
|
ψ
n
|
0
,
Ω
≤
(
1
+
|
χ
|
L
∞
(
Ω
)
∥
S
∥
∥
B
∥
)
|
ψ
n
|
0
,
Ω
.
The latter inequalities and (78) imply the following estimation:
(81)
|
λ
n
+
1
|
M
2
-
|
λ
n
|
M
2
≤
-
ρ
[
2
η
-
∥
C
∥
2
(
1
+
|
χ
|
L
∞
(
Ω
)
∥
S
∥
∥
B
∥
)
2
ρ
]
|
ψ
n
|
0
,
Ω
2
.
Hence, in particular, we get(82)
lim
⁡
n
→
∞
|
ψ
|
0
,
Ω
n
=
0
,
if ρ satisfies inequalities (75). Besides, we have
(83)
lim
⁡
n
→
∞
ξ
n
=
lim
⁡
n
→
∞
A
ψ
n
=
0
в
H
0
1
(
Ω
)
,
lim
⁡
n
→
∞
α
n
=
lim
⁡
n
→
∞
B
ψ
n
=
0
в
L
2
(
Γ
)
,
which finishes the proof.Since convergence of the considered method is guaranteed, any choice of subspaceM, satisfying the condition H
1
(
Ω
)
=
H
0
1
(
Ω
)
⊕
M and any choice of scalar product (
·
,
·
)
M on space is allowed. However, the choice influences parameter ρ, as well as the computation time on each iterations. Finally, we point out a few remarks regarding practical computations of λ
n
+
1
∈
M.If the scalar product inM is defined via the following formula:
(84)
(
λ
,
μ
)
M
=
∫
Ω
S
λ
·
S
μ
d
Ω
,
then as λ
n
+
1 any function from M can be taken, assuming that the following condition is satisfied:
(85)
λ
n
+
1
|
Γ
=
(
λ
n
-
ρ
θ
λ
n
)
|
Γ
.Equality (85) can be understood in the sense of trace equality on a boundary. In fact, if (85) is satisfied, then
(86)
(
λ
n
+
1
-
λ
n
,
μ
)
M
=
-
ρ
∫
Ω
S
θ
n
·
S
μ
d
Ω
=
-
ρ
∫
Ω
φ
θ
n
°
φ
μ
∘
d
Ω
,
which means that conditions of Theorem 7 are satisfied.We may choose also the following scalar product:(87)
either
(
λ
,
μ
)
M
=
∫
Ω
∇
λ
∇
μ
d
Ω
.
or
(
λ
,
μ
)
M
=
∫
Ω
λ
μ
d
Ω
.However, in the latter case one needs to compute gradient〈
D
′
(
λ
)
,
μ
〉 on each step, which extends the computational time.Final Remarks. (1) The proof has been carried out for equations in the hybrid form (2). It can be relatively easily extended into equations regarding displacements. (2) Results can be extended on other types of the differential equations, including nonlinear ones, consisting of a biharmonic operator.


## 3.3. Iterative Procedure for the Reduction of the Karman Equation into the Poisson Equation

In the preceding sections we have proved convergence of the iterative procedures for linearization of (2) by reducing the solution of the eighth order system of nonlinear differential equations into that of the solution to a biharmonic equation, as well as the reduction of the biharmonic equation to the Poisson-type equation in the case of a curvilinear boundary using the finite element method (FEM).While considering a space with the rectangular boundary, we may extend the procedure reported in Section3.1 by introduction of new variables into the iterative procedure of solution to the Poisson-type equations without difficulties.In the case of spaces with the curvilinear boundary, the procedure described in Section3.2 can be applied to solve (2) using the iterative procedure, whose convergence has been proved in Section 3.1.For this purpose new variablesM
w
(
x
,
y
) and M
F
(
x
,
y
) are introduced
(88)
M
w
(
x
,
y
)
=
Δ
w
(
x
,
y
)
,
M
F
(
x
,
y
)
=
Δ
F
(
x
,
y
)
.Then each of differential equations (2) is divided into two Poisson-type equations. The iterative procedure of solution of the obtained system of four Poisson-type equations has the following form:
(89)
Δ
M
w
(
k
)
=
q
+
L
(
w
(
k
-
1
)
,
F
(
k
-
1
)
)
,
Δ
w
(
k
)
=
M
w
(
k
)
,
Δ
M
F
(
k
)
=
L
(
w
(
k
)
,
w
(
k
)
)
,
Δ
F
(
k
)
=
M
F
(
k
)
,
{
x
,
y
}
∈
Ω
.Boundary conditions (4) are transformed to the following form:
(90)
w
|
Γ
=
M
w
|
Γ
=
F
|
Γ
=
M
F
|
Γ
=
0
.The given procedure (89) has advantages over procedure (6), while solving each equation since instead of the fourth order equation that of the second order is solved. Because equations are solved by numerical methods (FDM, FEM) and the approximation of the biharmonic operator has high requirements on the approximating functions, then for the Poisson-type equation one may simplify the procedure (89) of finding a solution by choosing simple approximating functions.In the FDM case, an order of algebraic equations system, after a discretization of the biharmonic equation, is higher for the second order equation, and hence higher expectations are required from computer abilities while solving the problem numerically.


## 4. The Method of Variational Iterations (MVI) of PDEs Solutions

### 4.1. Validation of Convergence

The method of variational iterations (MVI) was applied first in 1933 by Shunok who considered a deflection of cylindrical panels. However, this work did not meet with the response of others, and then it was rediscovered in the sixties of the previous century by Kantorovich and Krylov [20], who applied it in his investigation of rectangular plates. Then the MVI found wide application in solving various problems of plates and shells (see the list of references reported in [21]).Here we prove validity and reliability of the mentioned method for a class of equations with positively defined operators, that is, biharmonic and harmonic ones. In other words, we prove a theorem on convergence of the MVI for iterative procedures (6) and (89).Formally, the MVI scheme is as follows. Assume that our aim is to find a solution to the following:(91)
T
ω
(
x
,
y
)
=
g
(
x
,
y
)
,
x
,
y
∈
Ω
(
x
,
y
)
,
where T stands for a certain operator defined on set D
(
T
) of the Hilbert space L
2
(
Ω
); g
(
x
,
y
) is the function given for two variables x and y; ω
(
x
,
y
) is the function of these two variables being sought; Ω
(
x
,
y
) is the space of changes of variables x and y.IfΩ
(
x
,
y
)
=
X
×
Y (X is the certain bounded set of variables x, Y is the bounded set of y), then a solution to (91) can be given in the following form:
(92)
ω
N
(
x
,
y
)
=
∑
i
=
1
N
u
i
(
x
)
v
i
(
y
)
,
where functions u
i
(
x
) and v
i
(
y
) are defined by the following system of equations:
(93)
∫
X
(
T
ω
N
-
g
)
u
1
(
x
)
d
x
=
0
,
⋮
∫
X
(
T
ω
N
-
g
)
u
N
(
x
)
d
x
=
0
,
∫
Y
(
T
ω
N
-
g
)
v
1
(
y
)
d
y
=
0
,
⋮
∫
Y
(
T
ω
N
-
g
)
v
N
(
y
)
d
y
=
0
.It is found in the following way: we have a certain system composed ofN functions regarding one variable, for instance, u
1
0
(
x
)
,
u
2
0
(
x
)
,
…
,
u
N
0
(
x
), and then from the first N equations of system (93) the system of N functions v
1
1
(
x
)
,
v
2
1
(
x
)
,
…
,
v
N
1
(
x
) is defined. Then, the so far obtained functions represent a new choice of the functions regarding the variable x
-
u
1
2
(
x
)
,
u
2
2
(
x
)
,
…
,
u
N
2
(
x
), and the latter serves to get a new set of functions regarding variable y
-
v
1
3
(
x
)
,
v
2
3
(
x
)
,
…
,
v
N
3
(
x
), and so forth.Definition 9.
We say that a process of computation, when one given system of functions is replaced by the second system, is the MVI step. The number of steps needed to define a certain choice of functions corresponds to the superscript (number) of functions being considered. Truncating the process of finding functionsu
i
(
x
) and v
i
(
y
) on the kth step, which, for example, corresponds to the choice of functions v
1
k
(
y
)
,
v
2
k
(
y
)
,
…
,
v
N
k
(
y
), we define the function
(94)
ω
N
k
=
∑
i
=
1
N
u
i
k
-
1
(
x
)
v
i
k
(
y
)
,
taken as the approximating solution of (91) obtained by MVI.Remark 10.
Here and further on, we shall take as operatorT a certain differential operator defined on set D
(
T
) of the Hilbert space L
2
(
Ω
). Then, on each step system (93) shall be transformed to a system of ODEs which can be solved further.Remark 11.
We call functionω
N
(
x
,
y
) the Nth approximation to (91) if the number of series terms in (92) is equal to N.Let us study the case of first approximation; that is, the following solution of (91) is sought:
(95)
ω
1
(
x
,
y
)
=
u
(
x
)
v
(
y
)
,
where functions u
(
x
) and v
(
y
) are defined through the illustrated way from the following system of equations:
(96)
∫
X
(
T
u
(
x
)
·
v
(
y
)
-
g
)
u
(
x
)
d
x
=
0
,
∫
Y
(
T
u
(
x
)
·
v
(
y
)
-
g
)
u
(
y
)
d
y
=
0
.Let the operatorT in (91) be positive definite. Let us introduce the following notation: H
T
(
X
×
Y
) is the energy space of the operator T; [
·
,
·
] is the scalar product of elements in H
T; ω
0 is the exact solution to (91).Theorem 12.
IfT is a positive definite operator with the space of action D
(
T
)
⊂
H
T, then the sequence of elements α
k
=
∥
ω
1
k
(
x
,
y
)
-
ω
0
∥
H
T is monotonously decreasing; that is, for arbitrary i and j if i
≥
j, then
(97)
∥
ω
1
i
-
ω
0
∥
H
T
≤
∥
ω
1
j
-
ω
0
∥
H
T
.Proof.
We consider a subsetM
1
1 of the space H
T which has the following form:
(98)
M
1
1
=
{
ω
(
x
,
y
)
∣
ω
(
x
,
y
)
=
u
0
(
x
)
v
(
y
)
,
f
u
0
(
x
)
∈
H
T
(
x
)
,
v
(
y
)
∈
H
T
(
Y
)
}
.
It is clear, that setM
1
1 represents a subspace of space H
T
(
X
×
Y
) (generally, of infinite dimension). Therefore, one may define ω
0 projection onto space M
1
1. As it is known that element u
0
(
x
)
v
*
(
y
)
∈
M
1
1 stands for the projection of ω
0 onto M
1
1 if the following condition is satisfied:
(99)
[
u
0
(
x
)
v
*
(
y
)
-
ω
0
,
u
0
(
x
)
v
(
y
)
]
H
T
=
0
for arbitrary elements u
0
(
x
)
v
(
y
)
∈
M
1
1. It is clear that if u
0
(
x
)
v
*
(
y
)
∈
M
1
1, then (99) coincides with the first equation of system (97).
Since the elementu
0
(
x
)
v
1
(
y
) obtained through the first step of MVI is a projection of element ω
0 onto the subspace M
1
1, hence the following inequality holds:
(100)
∥
ω
1
1
(
x
,
y
)
-
ω
0
∥
H
T
≤
∥
u
0
(
x
)
v
(
y
)
-
ω
0
∥
H
T
for arbitrary elements u
0
(
x
)
v
(
y
)
∈
M
1
1. An analogous construction allows us to get a similar inequality for the subspaces; that is, we have
(101)
M
2
1
=
{
ω
(
x
,
y
)
∣
ω
(
x
,
y
)
=
u
0
(
x
)
v
1
(
y
)
,
f
f
u
(
x
)
∈
H
T
(
x
)
,
v
1
(
y
)
∈
H
T
(
Y
)
}
.
In the case corresponding to the second MVI step,(102)
∥
u
2
(
x
)
v
1
(
y
)
-
ω
0
∥
H
T
≤
∥
u
(
x
)
v
1
(
y
)
-
ω
0
∥
H
T
for arbitrary elements u
(
x
)
v
1
(
y
)
∈
M
2
1. It follows from (100) and (102) that ∥
u
2
(
x
)
v
1
(
y
)
-
ω
0
∥
H
T
≤
∥
u
0
(
x
)
v
1
(
y
)
-
ω
0
∥
H
T. Considerations similar to those so far provided and obtained for the kth MVI step prove the theorem as well as inequality (97) with the help of induction.Remark 13.
Results of Theorem12 are extended into the case of Nth approximation, and therefore inequality (97) is given in the following form:
(103)
∥
ω
N
n
(
x
,
y
)
-
ω
0
∥
H
T
≤
∥
ω
N
n
(
x
,
y
)
-
ω
0
∥
H
T
,
m
≥
n
.In order to prove the theorem we introduce the following lemma.Lemma 14.
Let each of elements of the basis system of spaceH
T have the following form:
(104)
θ
i
(
x
,
y
)

=

φ
i
(
x
)
ψ
i
(
y
)
,

∀
i

φ
i
(
x
)

∈

H
T
(
X
)
,
φ
i
(
y
)

∈

H
T
(
Y
)
.If for the initial MVI approximation one takes any component of a certain basis functionθ
i, that is, u
0
(
x
)
≡
φ
i
(
x
), then for an arbitrary number k of the MVI steps the following inequality holds:
(105)
∥
ω
1
k
(
x
,
y
)
-
ω
0
∥
H
T
≤
∥
c
φ
i
(
x
)
ψ
i
(
y
)
-
ω
0
∥
H
T
,
where c is the arbitrarily taken real number.Proof.
Sinceu
0
(
x
)
v
1
(
y
)
≡
φ
i
(
x
)
v
1
(
y
), then Theorem 12 yields
(106)
∥
ω
1
k
(
x
,
y
)
-
ω
0
∥
H
T
≤
∥
φ
i
(
x
)
v
1
(
y
)
-
ω
0
∥
H
T
≤
∥
c
φ
i
(
x
)
ψ
i
(
y
)
-
ω
0
∥
H
T
.
On the basis of the given lemma we formulate one of the MVI convergence criterions. Initially, we identify spaceH
T with space W
2
m
0
(
Ω
) which is generated through a closure regarding the norm
(107)
∥
ω
∥
W
2
m
=
{
∫
Ω
∑
k
=
0
m
∑
(
k
)
|
∂
k
ω
∂
k
1
x
∂
k
2
y
|
d
x
d
y
}
2
of a set of infinitely differentiable functions C
∞
0
(
Ω
) with a compact carrier inΩ.Theorem 15.
Let each of elements of the basis system of spaceW
2
m
0
(
X
×
Y
) have the following form:
(108)
θ
i
(
x
,
y
)
=
φ
i
(
x
)
ψ
i
(
y
)
,
where {
φ
i
(
x
)
} stands for the basis system in space W
2
m
0
(
X
) and {
ψ
i
(
y
)
} in space W
2
m
0
(
Y
), and in order to get an arbitrary Nth order approximation of the MVI we take components of the elements of the basis system {
θ
i
(
x
,
y
)
} as the initial functions. Then, for sufficiently large N the MVI gives a unique approximate solution ω
N, and the sequence {
ω
N
} is convergent with respect to the norm of space W
2
m
0
(
X
×
Y
) to the exact solution ω
0 irrespectively of the number of steps k. This construction can be carried out for each Nth approximation; that is,
(109)
∥
ω
N
k
-
ω
0
∥
W
2
m
0
⟶
0
,
N
⟶
∞
.Proof.
If we prove the theorem regarding approximations obtained via the first step of MVI, then owing to the Lemma the obtained results shall be valid for the arbitrarykth step. Therefore, let us consider the Nth approximation of problem (91), obtained on the first step. In a way similar to considerations regarding Theorem 12, one may show that each ω
N
1 stands for a projection of element ω
0 onto the following subspace:
(110)
M
N
1
=
{
ω
(
x
,
y
)
∣
ω
(
x
,
y
)
=
∑
i
=
1
N
u
0
(
x
)
v
1
(
y
)
}
,
where u
i
0
(
x
) denotes N the fixed elements from system {
φ
i
(
x
)
}, and for arbitrary i and v
i
(
y
) they cover the whole space W
2
m
0
(
Y
). Therefore, ω
N
1
=
P
N
ω
0, where P
N stands for the operator of an orthogonal projection onto subspace M
N
1 which is bounded. Since elements of the basis system θ
i
(
x
,
y
) have form (108), it is limiting dense ([20], page 191) in W
2
m
0
(
X
×
Y
). Then the proof is carried out in a way similar to that of Theorem 16.2 (see [20], page 216), since all conditions of its application are satisfied.Remark 16.
Results of Theorem12 and the Lemma indicate that using the MVI one may get an approximate solution to (91) in a way not worse than that obtained via the Ritz method in accordance with the corresponding subspace.
The MVI can be extended also on the case of a large number of variables. For instance, a sought solution to (91) is the function of three variables x
,
y
,
z, and hence the approximate solution of MVI can be sought in the following form:
(111)
ω
(
x
,
y
,
z
)
=
u
(
x
)
v
(
y
,
z
)
.Remark 17.
Note that during application of the MVI there is no need to construct the initial condition, satisfying, say, boundary conditions of the stated problem. Let us assume that operatorT defines a certain boundary value problem. Let us introduce an arbitrary function from a space of the definition of the differential operator of the studied problem. Then on the first (second) step we get a system of functions satisfying boundary conditions regarding one (two) of the variables.Observe that the MVI, on each step, defines only one of the functions appearing in the representation of solution (108). The following method develops MVI and allows us, using only a first step, to estimate at once two functions regarding two directions of the coordinates.


### 4.2. Numerical Results

Let us consider the following rectangular plate(112)
Δ
2
ω
=
g
(
x
,
y
)
D
,
where ω
(
x
,
y
) is the normal plate deflection in point x
,
y; g
(
x
,
y
) is the intensity of the normal load; D
=
E
h
3
[
12
(
1
-
ν
2
)
]
-
1; E
,
ν are the Young modulus and Poisson constant, respectively; 2
h is the plate thickness, y
-
=
0; 1
,
x
-
=
0
;
1, where y
-
=
y
/
b, x
-
=
x
/
a; a and b are the plate dimensions; and g
-
(
x
,
y
)
=
g
a
2
b
2
[
E
h
4
]
-
1. Plate space in the plane x
,
y is denoted by Ω, and its contour is Γ. Below, we study two types of boundary conditions.In the case of a simple supportw
=
∂
2
w
/
∂
n
2
|
Γ
=
0, a solution obtained via variational iterations is compared with that obtained through the first order approximation of the Bubnov method and with the solution represented by double trigonometric series. The following deflection function has been assumed: W
(
x
,
y
)
=
A
sin
(
π
x
)
sin
(
π
y
), which satisfies the boundary conditions. Substituting this into (112) and applying the Bubnov procedure we obtain A
=
(
4
/
π
6
)
q, and for q
=
1 we have A
=
0.00416. The deflection function w
(
x
,
y
)
=
∑
m
,
n
A
m
,
n
sin
⁡
(
m
π
x
/
a
)
sin
⁡
(
n
π
y
/
b
) is substituted into (20) then multiplied by sin
⁡
(
m
π
x
/
a
), sin
⁡
(
n
π
y
/
b
), and integrated regarding the plate surface. The following deflection value is obtained:
(113)
w
(
x
,
y
)
=
q
24
D
(
x
4
-
2
a
x
3
+
a
3
x
)
-
4
q
a
4
π
5
D
sin
m
π
x
a
×
∑
m
,
n
1
m
5
[
a
m
t
h
(
α
m
)
+
2
2
c
h
(
α
m
)
c
h
(
2
α
m
y
b
)
ddddddddddd
-
α
m
2
c
h
(
α
m
)
2
y
b
s
h
(
2
α
m
y
b
)
]
.It should be emphasized that the obtained series converges fast, and in practice it is sufficient to keep only its first term. For the square plate, the deflection measured in its center is 0.00406 (Table1).Table 1
Solution methods.
Reduction to the Poisson-type equation (MVI)
Bubnov’s method
Solution in series
0.004054
0.00416
0.00406In the second case the plate contours are clamped; that is,w
=
(
∂
w
/
∂
n
)
|
Γ
=
0. Here computations are carried out in the first approximation, where owing to Remark 16, we take sin
(
π
x
) as the input function; that is, this function does not satisfy damping conditions on the plate’s contour. A solution to the obtained ODEs was carried out via the difference method (FDM) with the plate partition 60
×
60 and successive solution to the obtained algebraic equations by the Gauss method.Results of the quarter plate deflection function obtained on the liney
-
=
0.5 are given in Table 2. These results coincide with the conclusion of Theorem 12 regarding monotonous series {
a
k
} behavior. The exact value of deflection equal to 0.0138 is taken from monograph [20].Table 2
Deflection of the plate quarter.
Step
y
-
=
0.5
,
x
-
=
0.1
y
-
=
0.5
,
x
-
=
0.3
y
-
=
0.5
,
x
-
=
0.5
1
0.55543336
·
1
0
-
3
0.19404364
·
1
0
-
2
0.22896260
·
1
0
-
2
2
0.16923198
·
1
0
-
2
0.79178517
·
1
0
-
2
0.10638649
·
1
0
-
1
3
0.21179814
·
1
0
-
2
0.10124784
·
1
0
-
1
0.13723563
·
1
0
-
1
4
0.21349905
·
1
0
-
2
0.102096
·
1
0
-
1
0.13840985
·
1
0
-
1
5
0.21349905
·
1
0
-
2
0.10211107
·
1
0
-
1
0.13842917
·
1
0
-
1
6
0.21352879
·
1
0
-
2
0.10211157
·
1
0
-
1
0.13843016
·
1
0
-
1


## 4.1. Validation of Convergence

The method of variational iterations (MVI) was applied first in 1933 by Shunok who considered a deflection of cylindrical panels. However, this work did not meet with the response of others, and then it was rediscovered in the sixties of the previous century by Kantorovich and Krylov [20], who applied it in his investigation of rectangular plates. Then the MVI found wide application in solving various problems of plates and shells (see the list of references reported in [21]).Here we prove validity and reliability of the mentioned method for a class of equations with positively defined operators, that is, biharmonic and harmonic ones. In other words, we prove a theorem on convergence of the MVI for iterative procedures (6) and (89).Formally, the MVI scheme is as follows. Assume that our aim is to find a solution to the following:(91)
T
ω
(
x
,
y
)
=
g
(
x
,
y
)
,
x
,
y
∈
Ω
(
x
,
y
)
,
where T stands for a certain operator defined on set D
(
T
) of the Hilbert space L
2
(
Ω
); g
(
x
,
y
) is the function given for two variables x and y; ω
(
x
,
y
) is the function of these two variables being sought; Ω
(
x
,
y
) is the space of changes of variables x and y.IfΩ
(
x
,
y
)
=
X
×
Y (X is the certain bounded set of variables x, Y is the bounded set of y), then a solution to (91) can be given in the following form:
(92)
ω
N
(
x
,
y
)
=
∑
i
=
1
N
u
i
(
x
)
v
i
(
y
)
,
where functions u
i
(
x
) and v
i
(
y
) are defined by the following system of equations:
(93)
∫
X
(
T
ω
N
-
g
)
u
1
(
x
)
d
x
=
0
,
⋮
∫
X
(
T
ω
N
-
g
)
u
N
(
x
)
d
x
=
0
,
∫
Y
(
T
ω
N
-
g
)
v
1
(
y
)
d
y
=
0
,
⋮
∫
Y
(
T
ω
N
-
g
)
v
N
(
y
)
d
y
=
0
.It is found in the following way: we have a certain system composed ofN functions regarding one variable, for instance, u
1
0
(
x
)
,
u
2
0
(
x
)
,
…
,
u
N
0
(
x
), and then from the first N equations of system (93) the system of N functions v
1
1
(
x
)
,
v
2
1
(
x
)
,
…
,
v
N
1
(
x
) is defined. Then, the so far obtained functions represent a new choice of the functions regarding the variable x
-
u
1
2
(
x
)
,
u
2
2
(
x
)
,
…
,
u
N
2
(
x
), and the latter serves to get a new set of functions regarding variable y
-
v
1
3
(
x
)
,
v
2
3
(
x
)
,
…
,
v
N
3
(
x
), and so forth.Definition 9.
We say that a process of computation, when one given system of functions is replaced by the second system, is the MVI step. The number of steps needed to define a certain choice of functions corresponds to the superscript (number) of functions being considered. Truncating the process of finding functionsu
i
(
x
) and v
i
(
y
) on the kth step, which, for example, corresponds to the choice of functions v
1
k
(
y
)
,
v
2
k
(
y
)
,
…
,
v
N
k
(
y
), we define the function
(94)
ω
N
k
=
∑
i
=
1
N
u
i
k
-
1
(
x
)
v
i
k
(
y
)
,
taken as the approximating solution of (91) obtained by MVI.Remark 10.
Here and further on, we shall take as operatorT a certain differential operator defined on set D
(
T
) of the Hilbert space L
2
(
Ω
). Then, on each step system (93) shall be transformed to a system of ODEs which can be solved further.Remark 11.
We call functionω
N
(
x
,
y
) the Nth approximation to (91) if the number of series terms in (92) is equal to N.Let us study the case of first approximation; that is, the following solution of (91) is sought:
(95)
ω
1
(
x
,
y
)
=
u
(
x
)
v
(
y
)
,
where functions u
(
x
) and v
(
y
) are defined through the illustrated way from the following system of equations:
(96)
∫
X
(
T
u
(
x
)
·
v
(
y
)
-
g
)
u
(
x
)
d
x
=
0
,
∫
Y
(
T
u
(
x
)
·
v
(
y
)
-
g
)
u
(
y
)
d
y
=
0
.Let the operatorT in (91) be positive definite. Let us introduce the following notation: H
T
(
X
×
Y
) is the energy space of the operator T; [
·
,
·
] is the scalar product of elements in H
T; ω
0 is the exact solution to (91).Theorem 12.
IfT is a positive definite operator with the space of action D
(
T
)
⊂
H
T, then the sequence of elements α
k
=
∥
ω
1
k
(
x
,
y
)
-
ω
0
∥
H
T is monotonously decreasing; that is, for arbitrary i and j if i
≥
j, then
(97)
∥
ω
1
i
-
ω
0
∥
H
T
≤
∥
ω
1
j
-
ω
0
∥
H
T
.Proof.
We consider a subsetM
1
1 of the space H
T which has the following form:
(98)
M
1
1
=
{
ω
(
x
,
y
)
∣
ω
(
x
,
y
)
=
u
0
(
x
)
v
(
y
)
,
f
u
0
(
x
)
∈
H
T
(
x
)
,
v
(
y
)
∈
H
T
(
Y
)
}
.
It is clear, that setM
1
1 represents a subspace of space H
T
(
X
×
Y
) (generally, of infinite dimension). Therefore, one may define ω
0 projection onto space M
1
1. As it is known that element u
0
(
x
)
v
*
(
y
)
∈
M
1
1 stands for the projection of ω
0 onto M
1
1 if the following condition is satisfied:
(99)
[
u
0
(
x
)
v
*
(
y
)
-
ω
0
,
u
0
(
x
)
v
(
y
)
]
H
T
=
0
for arbitrary elements u
0
(
x
)
v
(
y
)
∈
M
1
1. It is clear that if u
0
(
x
)
v
*
(
y
)
∈
M
1
1, then (99) coincides with the first equation of system (97).
Since the elementu
0
(
x
)
v
1
(
y
) obtained through the first step of MVI is a projection of element ω
0 onto the subspace M
1
1, hence the following inequality holds:
(100)
∥
ω
1
1
(
x
,
y
)
-
ω
0
∥
H
T
≤
∥
u
0
(
x
)
v
(
y
)
-
ω
0
∥
H
T
for arbitrary elements u
0
(
x
)
v
(
y
)
∈
M
1
1. An analogous construction allows us to get a similar inequality for the subspaces; that is, we have
(101)
M
2
1
=
{
ω
(
x
,
y
)
∣
ω
(
x
,
y
)
=
u
0
(
x
)
v
1
(
y
)
,
f
f
u
(
x
)
∈
H
T
(
x
)
,
v
1
(
y
)
∈
H
T
(
Y
)
}
.
In the case corresponding to the second MVI step,(102)
∥
u
2
(
x
)
v
1
(
y
)
-
ω
0
∥
H
T
≤
∥
u
(
x
)
v
1
(
y
)
-
ω
0
∥
H
T
for arbitrary elements u
(
x
)
v
1
(
y
)
∈
M
2
1. It follows from (100) and (102) that ∥
u
2
(
x
)
v
1
(
y
)
-
ω
0
∥
H
T
≤
∥
u
0
(
x
)
v
1
(
y
)
-
ω
0
∥
H
T. Considerations similar to those so far provided and obtained for the kth MVI step prove the theorem as well as inequality (97) with the help of induction.Remark 13.
Results of Theorem12 are extended into the case of Nth approximation, and therefore inequality (97) is given in the following form:
(103)
∥
ω
N
n
(
x
,
y
)
-
ω
0
∥
H
T
≤
∥
ω
N
n
(
x
,
y
)
-
ω
0
∥
H
T
,
m
≥
n
.In order to prove the theorem we introduce the following lemma.Lemma 14.
Let each of elements of the basis system of spaceH
T have the following form:
(104)
θ
i
(
x
,
y
)

=

φ
i
(
x
)
ψ
i
(
y
)
,

∀
i

φ
i
(
x
)

∈

H
T
(
X
)
,
φ
i
(
y
)

∈

H
T
(
Y
)
.If for the initial MVI approximation one takes any component of a certain basis functionθ
i, that is, u
0
(
x
)
≡
φ
i
(
x
), then for an arbitrary number k of the MVI steps the following inequality holds:
(105)
∥
ω
1
k
(
x
,
y
)
-
ω
0
∥
H
T
≤
∥
c
φ
i
(
x
)
ψ
i
(
y
)
-
ω
0
∥
H
T
,
where c is the arbitrarily taken real number.Proof.
Sinceu
0
(
x
)
v
1
(
y
)
≡
φ
i
(
x
)
v
1
(
y
), then Theorem 12 yields
(106)
∥
ω
1
k
(
x
,
y
)
-
ω
0
∥
H
T
≤
∥
φ
i
(
x
)
v
1
(
y
)
-
ω
0
∥
H
T
≤
∥
c
φ
i
(
x
)
ψ
i
(
y
)
-
ω
0
∥
H
T
.
On the basis of the given lemma we formulate one of the MVI convergence criterions. Initially, we identify spaceH
T with space W
2
m
0
(
Ω
) which is generated through a closure regarding the norm
(107)
∥
ω
∥
W
2
m
=
{
∫
Ω
∑
k
=
0
m
∑
(
k
)
|
∂
k
ω
∂
k
1
x
∂
k
2
y
|
d
x
d
y
}
2
of a set of infinitely differentiable functions C
∞
0
(
Ω
) with a compact carrier inΩ.Theorem 15.
Let each of elements of the basis system of spaceW
2
m
0
(
X
×
Y
) have the following form:
(108)
θ
i
(
x
,
y
)
=
φ
i
(
x
)
ψ
i
(
y
)
,
where {
φ
i
(
x
)
} stands for the basis system in space W
2
m
0
(
X
) and {
ψ
i
(
y
)
} in space W
2
m
0
(
Y
), and in order to get an arbitrary Nth order approximation of the MVI we take components of the elements of the basis system {
θ
i
(
x
,
y
)
} as the initial functions. Then, for sufficiently large N the MVI gives a unique approximate solution ω
N, and the sequence {
ω
N
} is convergent with respect to the norm of space W
2
m
0
(
X
×
Y
) to the exact solution ω
0 irrespectively of the number of steps k. This construction can be carried out for each Nth approximation; that is,
(109)
∥
ω
N
k
-
ω
0
∥
W
2
m
0
⟶
0
,
N
⟶
∞
.Proof.
If we prove the theorem regarding approximations obtained via the first step of MVI, then owing to the Lemma the obtained results shall be valid for the arbitrarykth step. Therefore, let us consider the Nth approximation of problem (91), obtained on the first step. In a way similar to considerations regarding Theorem 12, one may show that each ω
N
1 stands for a projection of element ω
0 onto the following subspace:
(110)
M
N
1
=
{
ω
(
x
,
y
)
∣
ω
(
x
,
y
)
=
∑
i
=
1
N
u
0
(
x
)
v
1
(
y
)
}
,
where u
i
0
(
x
) denotes N the fixed elements from system {
φ
i
(
x
)
}, and for arbitrary i and v
i
(
y
) they cover the whole space W
2
m
0
(
Y
). Therefore, ω
N
1
=
P
N
ω
0, where P
N stands for the operator of an orthogonal projection onto subspace M
N
1 which is bounded. Since elements of the basis system θ
i
(
x
,
y
) have form (108), it is limiting dense ([20], page 191) in W
2
m
0
(
X
×
Y
). Then the proof is carried out in a way similar to that of Theorem 16.2 (see [20], page 216), since all conditions of its application are satisfied.Remark 16.
Results of Theorem12 and the Lemma indicate that using the MVI one may get an approximate solution to (91) in a way not worse than that obtained via the Ritz method in accordance with the corresponding subspace.
The MVI can be extended also on the case of a large number of variables. For instance, a sought solution to (91) is the function of three variables x
,
y
,
z, and hence the approximate solution of MVI can be sought in the following form:
(111)
ω
(
x
,
y
,
z
)
=
u
(
x
)
v
(
y
,
z
)
.Remark 17.
Note that during application of the MVI there is no need to construct the initial condition, satisfying, say, boundary conditions of the stated problem. Let us assume that operatorT defines a certain boundary value problem. Let us introduce an arbitrary function from a space of the definition of the differential operator of the studied problem. Then on the first (second) step we get a system of functions satisfying boundary conditions regarding one (two) of the variables.Observe that the MVI, on each step, defines only one of the functions appearing in the representation of solution (108). The following method develops MVI and allows us, using only a first step, to estimate at once two functions regarding two directions of the coordinates.


## 4.2. Numerical Results

Let us consider the following rectangular plate(112)
Δ
2
ω
=
g
(
x
,
y
)
D
,
where ω
(
x
,
y
) is the normal plate deflection in point x
,
y; g
(
x
,
y
) is the intensity of the normal load; D
=
E
h
3
[
12
(
1
-
ν
2
)
]
-
1; E
,
ν are the Young modulus and Poisson constant, respectively; 2
h is the plate thickness, y
-
=
0; 1
,
x
-
=
0
;
1, where y
-
=
y
/
b, x
-
=
x
/
a; a and b are the plate dimensions; and g
-
(
x
,
y
)
=
g
a
2
b
2
[
E
h
4
]
-
1. Plate space in the plane x
,
y is denoted by Ω, and its contour is Γ. Below, we study two types of boundary conditions.In the case of a simple supportw
=
∂
2
w
/
∂
n
2
|
Γ
=
0, a solution obtained via variational iterations is compared with that obtained through the first order approximation of the Bubnov method and with the solution represented by double trigonometric series. The following deflection function has been assumed: W
(
x
,
y
)
=
A
sin
(
π
x
)
sin
(
π
y
), which satisfies the boundary conditions. Substituting this into (112) and applying the Bubnov procedure we obtain A
=
(
4
/
π
6
)
q, and for q
=
1 we have A
=
0.00416. The deflection function w
(
x
,
y
)
=
∑
m
,
n
A
m
,
n
sin
⁡
(
m
π
x
/
a
)
sin
⁡
(
n
π
y
/
b
) is substituted into (20) then multiplied by sin
⁡
(
m
π
x
/
a
), sin
⁡
(
n
π
y
/
b
), and integrated regarding the plate surface. The following deflection value is obtained:
(113)
w
(
x
,
y
)
=
q
24
D
(
x
4
-
2
a
x
3
+
a
3
x
)
-
4
q
a
4
π
5
D
sin
m
π
x
a
×
∑
m
,
n
1
m
5
[
a
m
t
h
(
α
m
)
+
2
2
c
h
(
α
m
)
c
h
(
2
α
m
y
b
)
ddddddddddd
-
α
m
2
c
h
(
α
m
)
2
y
b
s
h
(
2
α
m
y
b
)
]
.It should be emphasized that the obtained series converges fast, and in practice it is sufficient to keep only its first term. For the square plate, the deflection measured in its center is 0.00406 (Table1).Table 1
Solution methods.
Reduction to the Poisson-type equation (MVI)
Bubnov’s method
Solution in series
0.004054
0.00416
0.00406In the second case the plate contours are clamped; that is,w
=
(
∂
w
/
∂
n
)
|
Γ
=
0. Here computations are carried out in the first approximation, where owing to Remark 16, we take sin
(
π
x
) as the input function; that is, this function does not satisfy damping conditions on the plate’s contour. A solution to the obtained ODEs was carried out via the difference method (FDM) with the plate partition 60
×
60 and successive solution to the obtained algebraic equations by the Gauss method.Results of the quarter plate deflection function obtained on the liney
-
=
0.5 are given in Table 2. These results coincide with the conclusion of Theorem 12 regarding monotonous series {
a
k
} behavior. The exact value of deflection equal to 0.0138 is taken from monograph [20].Table 2
Deflection of the plate quarter.
Step
y
-
=
0.5
,
x
-
=
0.1
y
-
=
0.5
,
x
-
=
0.3
y
-
=
0.5
,
x
-
=
0.5
1
0.55543336
·
1
0
-
3
0.19404364
·
1
0
-
2
0.22896260
·
1
0
-
2
2
0.16923198
·
1
0
-
2
0.79178517
·
1
0
-
2
0.10638649
·
1
0
-
1
3
0.21179814
·
1
0
-
2
0.10124784
·
1
0
-
1
0.13723563
·
1
0
-
1
4
0.21349905
·
1
0
-
2
0.102096
·
1
0
-
1
0.13840985
·
1
0
-
1
5
0.21349905
·
1
0
-
2
0.10211107
·
1
0
-
1
0.13842917
·
1
0
-
1
6
0.21352879
·
1
0
-
2
0.10211157
·
1
0
-
1
0.13843016
·
1
0
-
1


## 5. Numerical Study of the Karman Equations

### 5.1. Iterative Procedure of Linearization and Variational Iteration for System (6)

A simultaneous use of the iterative procedure is described in Section3 and the method of variational iterations of system (6) allows us to carry out three remarkable procedures:(1)
decrease of the order of the system twice (from the 8th to 4th order);
(2)
linearization of the sought nonlinear systems;
(3)
transition from PDEs to ODEs with constant coefficients.This result is particularly important in the analysis of elliptic type PDEs.Next, we present numerical results of our method using an example of the computation of flexible isotropic square plates of constant thickness for three types of boundary conditions (114)–(116).Consider the following:(114)
w
=
∂
2
w
∂
n
2
=
F
=
∂
2
F
∂
n
2
=
0
,
x
=
y
=
0
,
x
=
y
=
1
,
(115)
w
=
∂
2
w
∂
n
2
=
F
=
∂
F
∂
n
=
0
;
x
=
y
=
0
;
x
=
y
=
1
,
(116)
w
=
∂
w
∂
n
=
F
=
∂
F
∂
n
=
0
;
x
=
y
=
0
;
x
=
y
=
1
.For simplicity, we apply the MVI using its first approximationN
=
1. ODEs are reduced to AE (algebraic equations) through FDM with approximation 0(h
2), which is solved using the Gauss method. The interval of integration [
0
,
1
] was divided into 100 parts. Relation q
[
w
(
0.5,0.5
)
] is illustrated in Figure 2. Curves (
1
), (
2
), and (
3
) refer to boundary conditions (114), (115), and (116), respectively. Curves (
2
) and (
3
) are obtained for the Poisson coefficient ν
=
0.33 and curve (
1
) for ν
=
0.1. Circles refer to experimental results [21]; stars refer to the solution obtained by FDM [22], where FDM was applied directly to (2) and nonlinear AEs was solved by Newton’s method. Plane mesh step is 20 × 20. Computations were carried out with step Δ
q
=
10, where in order to accelerate convergence of the iterative procedure, w and F were taken from the previous step.Figure 2
q versus w.Dependence of the deflection change in the plate center on the number of iterations is shown in Figure3 (curve (
1
)). Boundary conditions correspond to the plate support on flexible noncompressed ribs (114). The remaining parameters are N
=
20, q
=
60, ν
=
0.28, and ε
=
1
0
-
3. We require 16 iterations to achieve a priori given accuracy.Figure 3
w versus n.One may see from Figure3 that the deflection oscillates in the vicinity of a certain averaged value, and it tends to it with an increase of the iteration number. It can be explained in the following way. Since the initial approximation is given by the linear equations, the observed deflection shall be larger than a real one. Substitution of w into the second equation of system (6) shows that the stress function value is also larger than a real one. Therefore, taking into account the obtained values of the stress function F, the deflection estimated through the first equation of system (6) is lesser than the real one. In other words, the deflection obtained via odd (even) iteration is lesser (larger) than the real one, and it can be estimated by the following formula:
(117)
w
=
w
odd
+
w
even
2
.Formula (117) allowed us to reduce the number of iterations up to seven (see Figure 3, curve (
2
)). However, an increase of the load implies the convergence decrease.


### 5.2. Iterative Linearization Procedure (Poisson-Type Equations)

In order to solve system (89) we used the MVI and FDM. Results of a comparison of solutions for systems (6) and (89), using the MVI and FDM applied to system (2) in the case of a square plate, are shown in Figure 4. One may see that the result obtained by the application of procedure (6)—curve (
2
)—and (89)—curve (3)—differs slightly from the result obtained via the FDM—(1). Furthermore, a comparison of the results obtained through iterative procedures (6) (dot curve) and (89) (dashed curve) shows that the compared results for small deflection practically coincide.Figure 4
Load deflection function.The use of formula (117), in order to increase the convergence, implies that the function “load-deflection” practically remains unaffected, but the number of iterations increases. Figure 5 shows results regarding the procedure (89): (i) without the application of (117)—curve (
1
)—and (ii) with the application of formula (117)—curve (
2
). One may observe that the iterative process converges twice as fast.Figure 5
Plate deflection versus number of iterations.


## 5.1. Iterative Procedure of Linearization and Variational Iteration for System (6)

A simultaneous use of the iterative procedure is described in Section3 and the method of variational iterations of system (6) allows us to carry out three remarkable procedures:(1)
decrease of the order of the system twice (from the 8th to 4th order);
(2)
linearization of the sought nonlinear systems;
(3)
transition from PDEs to ODEs with constant coefficients.This result is particularly important in the analysis of elliptic type PDEs.Next, we present numerical results of our method using an example of the computation of flexible isotropic square plates of constant thickness for three types of boundary conditions (114)–(116).Consider the following:(114)
w
=
∂
2
w
∂
n
2
=
F
=
∂
2
F
∂
n
2
=
0
,
x
=
y
=
0
,
x
=
y
=
1
,
(115)
w
=
∂
2
w
∂
n
2
=
F
=
∂
F
∂
n
=
0
;
x
=
y
=
0
;
x
=
y
=
1
,
(116)
w
=
∂
w
∂
n
=
F
=
∂
F
∂
n
=
0
;
x
=
y
=
0
;
x
=
y
=
1
.For simplicity, we apply the MVI using its first approximationN
=
1. ODEs are reduced to AE (algebraic equations) through FDM with approximation 0(h
2), which is solved using the Gauss method. The interval of integration [
0
,
1
] was divided into 100 parts. Relation q
[
w
(
0.5,0.5
)
] is illustrated in Figure 2. Curves (
1
), (
2
), and (
3
) refer to boundary conditions (114), (115), and (116), respectively. Curves (
2
) and (
3
) are obtained for the Poisson coefficient ν
=
0.33 and curve (
1
) for ν
=
0.1. Circles refer to experimental results [21]; stars refer to the solution obtained by FDM [22], where FDM was applied directly to (2) and nonlinear AEs was solved by Newton’s method. Plane mesh step is 20 × 20. Computations were carried out with step Δ
q
=
10, where in order to accelerate convergence of the iterative procedure, w and F were taken from the previous step.Figure 2
q versus w.Dependence of the deflection change in the plate center on the number of iterations is shown in Figure3 (curve (
1
)). Boundary conditions correspond to the plate support on flexible noncompressed ribs (114). The remaining parameters are N
=
20, q
=
60, ν
=
0.28, and ε
=
1
0
-
3. We require 16 iterations to achieve a priori given accuracy.Figure 3
w versus n.One may see from Figure3 that the deflection oscillates in the vicinity of a certain averaged value, and it tends to it with an increase of the iteration number. It can be explained in the following way. Since the initial approximation is given by the linear equations, the observed deflection shall be larger than a real one. Substitution of w into the second equation of system (6) shows that the stress function value is also larger than a real one. Therefore, taking into account the obtained values of the stress function F, the deflection estimated through the first equation of system (6) is lesser than the real one. In other words, the deflection obtained via odd (even) iteration is lesser (larger) than the real one, and it can be estimated by the following formula:
(117)
w
=
w
odd
+
w
even
2
.Formula (117) allowed us to reduce the number of iterations up to seven (see Figure 3, curve (
2
)). However, an increase of the load implies the convergence decrease.


## 5.2. Iterative Linearization Procedure (Poisson-Type Equations)

In order to solve system (89) we used the MVI and FDM. Results of a comparison of solutions for systems (6) and (89), using the MVI and FDM applied to system (2) in the case of a square plate, are shown in Figure 4. One may see that the result obtained by the application of procedure (6)—curve (
2
)—and (89)—curve (3)—differs slightly from the result obtained via the FDM—(1). Furthermore, a comparison of the results obtained through iterative procedures (6) (dot curve) and (89) (dashed curve) shows that the compared results for small deflection practically coincide.Figure 4
Load deflection function.The use of formula (117), in order to increase the convergence, implies that the function “load-deflection” practically remains unaffected, but the number of iterations increases. Figure 5 shows results regarding the procedure (89): (i) without the application of (117)—curve (
1
)—and (ii) with the application of formula (117)—curve (
2
). One may observe that the iterative process converges twice as fast.Figure 5
Plate deflection versus number of iterations.


## 6. Summary

In this work we have proposed and theoretically established (theorems with proofs) some iterative procedures dedicated to a decrease of the order and then linearization of the Karman nonlinear PDEs. It has been shown that the result obtained via the modified Kantorovich-Vlasov method coincides with the exact solution. Furthermore, the theoretical considerations have been supported by the numerical analysis of the Karman equations using two introduced iterative procedures.



---

*Source: 792829-2014-03-16.xml*