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

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.

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 Section 3 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.

Section 6 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 space 𝑅³ (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𝜕𝑥2𝜕𝑤𝜕𝑥²,𝜀_𝑥𝑦=𝜕𝑢+𝜕𝑦𝜕𝑣+𝜕𝑥𝜕𝑤𝜕𝑥𝜕𝑤,(),.𝜕𝑦𝑢,𝑣,𝑤𝑥,𝑦 (1)

Figure 1 

A studied plate.

Equations governing the deflection 𝑤(𝑥,𝑦,𝑡) and stress function 𝐹(𝑥,𝑦) have the following form [15]: Δ²()Δ𝑤−𝐿𝑤,𝐹−𝑞=0,₂1𝐹+2𝐿(𝑤,𝑤)=0. (2)

The following operators are introduced: Δ²(⋅)=1𝜆²𝜕⁴(⋅)𝜕𝑥⁴+𝜆²𝜕⁴(⋅)𝜕𝑦⁴𝜕+2⁴(⋅)𝜕𝑥²𝜕𝑦²,𝐿()=𝜕⋅,⋅²(⋅)𝜕𝑥²𝜕²(⋅)𝜕𝑦²𝜕−2²(⋅)𝜕𝜕𝑥𝜕𝑦²(⋅)+𝜕𝜕𝑥𝜕𝑦²(⋅)𝜕𝑥²𝜕²(⋅)𝜕𝑦². (3)

Here and further on the following nondimensional quantities are introduced: 𝑤=ℎ𝑤; 𝐹=𝐸ℎ²𝐹; 𝑡=𝑡₀𝑡; 𝜀=𝜀/𝜏; 𝑥=𝑎𝑥; 𝑦=𝑏𝑦; 𝑞=𝑞(12(1−𝜈²)𝐸ℎ⁴/𝑎²𝑏²); 𝜏=(𝑎𝑏/ℎ)𝜌/𝐸𝑔; 𝜆=𝑎/𝑏, where 𝑎,𝑏 are the maximal plate dimensions regarding 𝑥 and 𝑦, respectively; ℎ is thickness; 𝑔 is acceleration due to gravity; 𝜌=𝛾ℎ; 𝛾 is specific gravity of volume plate material; 𝜈 is Poisson’s coefficient; 𝐸 is the Young modulus; 𝑤,𝐹 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]: 𝑤|_Γ=𝜕²𝑤𝜕𝑛²_Γ=𝐹|_Γ=𝜕²𝐹𝜕𝑛²_Γ=0, (4) where Γ stands for the space boundary occupied by the plate. The following initial conditions are attached to (2): 𝑤|_𝑡=0=𝑤₀,𝑤̇|_𝑡=0=𝑤̇₀. (5)

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: Δ²𝑤^(𝑘)=𝐿𝑤^(𝑘−1),𝐹^(𝑘−1)Δ+𝑞,₂𝐹^(𝑘)1=−2𝐿𝑤^(𝑘),𝑤^(𝑘),𝑥,𝑦∈Ω. (6)

On the first step of the iterative procedure the following biharmonic equation for a given load 𝑞(𝑥,𝑦) is solved: Δ²𝑤⁽¹⁾𝑥,𝑦=𝑞𝑥,𝑦. (7)

The value of 𝑤⁽¹⁾(𝑥,𝑦) is substituted into the right-hand side of equation system (6), and as a result a biharmonic equation for 𝐹⁽¹⁾(𝑥,𝑦) 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. Let 𝐻²(Ω) refer to a Sobolev space of functions 𝜉={𝑤,𝐹} such that 𝜉∈𝐿²(Ω),𝜕𝜉𝜕𝑥_𝑖∈𝐿²(Ω),𝜕²𝜉𝜕𝑥_𝑖𝜕𝑥_𝑗∈𝐿²(Ω);𝑖,𝑗=1,2, (8) where 𝐿²(Ω) denotes the space of functions being summed with a square in Ω.

Let 𝐻²₀(Ω) denote the closure of functions from 𝐷(Ω) (space of functions of class 𝐶^∞ in Ω, having compact carrier in Ω) in norm 𝐻²(Ω): 𝐻²₀(Ω)=𝐷(Ω)^𝐻2(Ω)=𝜉∈𝐻²(Ω)∣𝜉_Γ=𝜕𝜉𝜕𝑛_Γ=0. (9)

Since space Ω is bounded and its boundary Γ is efficiently regular, then map 𝜉→‖Δ𝜉‖_0,Ω defines the norm in 𝐻²₀(Ω) being equivalent to the norm generated by spaces 𝐻²(Ω).

Assume that 𝑞∈𝐻⁻²(Ω) (𝐻⁻²(Ω) denotes a conjugation to 𝐻²(Ω)). 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 𝐿²(Ω): (𝜉,𝜂)=∫_Ω𝜉𝜂𝑑Ω, and by 𝛽(𝑤,𝐹,𝜇) a three linear form defined on (𝐻²₀(Ω))³: 𝛽=+1𝑤,𝐹,𝜇Δ𝐹,Δ𝜇2𝐿(𝑤,𝑤),𝜇. (10)

Let us define the set 𝑀=𝑤,𝐹∈𝐻²₀(Ω)∣∀𝜇∈𝐻²₀(Ω),𝛽𝑤,𝐹,𝜇=0, (11) and square function 𝐽(𝑤,𝐹)∶𝑀→𝑅𝐽()=1𝑤,𝐹2‖Δ𝑤‖²_0,Ω+12‖Δ𝐹‖²_0,Ω−𝑞,𝑤. (12)

Theorem 1. The problem of minimizing (12) on set (11) has, at least, one solution.

Proof. Let {𝑤_𝑛,𝐹_𝑛}∈𝑀 be the minimizing series; that is, we have 𝐽𝑤_𝑛,𝐹_𝑛⟶inf_{{}𝑤,𝐹∈𝑀}𝐽(𝑤,𝐹), (13) which exists, since 𝐽 is a square functional.
For arbitrary 𝑤,𝐹∈𝐻²₀(Ω) the following inequality holds: 𝐽(𝑤,𝐹)≥𝑐₁‖𝑤‖²_2,Ω+𝑐₂‖𝐹‖²_2,Ω−𝑐₃⋅‖𝑤‖_2,Ω, (14) where ‖⋅‖_2,Ω denotes the norm in 𝐻²(Ω) and 𝑐_𝑖 are certain positive constants. Then, (13) yields 𝑐₁‖𝑤_𝑛‖²_2,Ω+𝑐₂‖𝐹_𝑛‖²_2,Ω−𝑐₃‖𝑤_𝑛‖_2,Ω≤𝐽(𝑤_𝑛,𝐹_𝑛)≤𝐽(𝑤₀,𝐹₀)=𝐴, where 𝑤₀,𝐹₀ are the arbitrary functions (initial approximation).
Then, the following estimation holds: 𝑐₁𝑤_𝑛²_2,Ω−𝑐₃2𝑐₁²+𝑐₂𝐹_𝑛²_2,Ω𝑐≤𝐴+²³4𝑐₁,𝑤_𝑛2,Ω≤𝑐₄,𝐹_𝑛2,Ω≤𝑐₅. (15)
Therefore, series {𝑤_𝑛,𝐹_𝑛} is bounded in (𝐻²₀(Ω))². Consequently, one may choose a series {𝑤_𝑘,𝐹_𝑘} that 𝑤_𝑘→𝑤, 𝐹_𝑘→𝐹 is weak in 𝐻²₀(Ω). Since 𝐻²₀(Ω)→𝐿²(Ω) is compact, then 𝑤_𝑘→𝑤, 𝐹_𝑘→𝐹 is strong in 𝐿²(Ω).
We show that interval {𝑤,𝐹} of the minimized series belongs to 𝑀; that is, 𝛽(𝑤,𝐹,𝜇)=0, for all 𝜇∈𝐻²₀(Ω).
Since (𝐿(𝑤_𝑘,𝑤_𝑘),𝜇)=(𝐿(𝑤_𝑘,𝜇),𝑤_𝑘)∀𝜇∈𝐻²₀(Ω), 𝐿(𝑤_𝑘,𝜇)→𝐿(𝑤,𝜇) is weak in 𝐻²₀(Ω), and 𝑤_𝑘→𝑤 is strong in 𝐿²(Ω), we get (𝐿(𝑤_𝑘,𝑤_𝑘),𝜇)=(𝐿(𝑤,𝑤),𝜇) and consequently, 𝛽(𝑤,𝐹,𝜇)=0 for all 𝜇∈𝐻²₀(Ω). This means that 𝑤,𝐹∈𝑀. (16)
However, 𝐽(𝑤,𝐹) is half-continuous from below in weak topology on (𝐻²(Ω))², and therefore the following inequality holds: lim_𝑘→∞𝐽(𝑤_𝑘,𝐹_𝑘)≥𝐽(𝑤,𝐹).
Then (13) and (16) imply that 𝐽(𝑤,𝐹)≤inf_{(𝑤,𝐹)∈𝑀}𝐽(𝑤,𝐹). Therefore, the following equation holds: 𝐽(𝑤,𝐹)=inf_{(𝑤,𝐹)∈𝑀}𝐽(𝑤,𝐹), which means that {𝑤,𝐹}∈𝑀 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 {𝑤,𝐹}∈𝑀, satisfying the following: Δ𝑤,Δ𝜇−𝐿(𝑤,𝐹),𝜇=𝑞,𝜇∀𝜇∈𝐻²₀(Ω). (17)

Theorem 2. Points of the functional minimum (12) are weak solutions to problems (6) and (4).

Proof. Let {𝑤,𝐹}∈𝑀 be one of the functional (12) minimum points. Let us take 𝜂=𝑤+𝑡𝛿𝑤𝛿𝑤∈𝐻²₀(Ω) and let us choose 𝜉=𝐹+𝛿𝐹𝛿𝐹∈𝐻²₀(Ω) such that {𝜂,𝜉}∈𝑀, that is, in the way that 𝛽(𝑤,𝐹,𝜇)=0 for all 𝜇∈𝐻²₀(Ω). Then 𝐽(𝑤,𝐹)≤𝐽(𝜂,𝜉). This yields 𝑡()+()+𝑡Δ𝑤,Δ𝛿𝑤Δ𝐹,Δ𝛿𝐹−𝑡𝑞,𝛿𝑤²2‖Δ𝛿𝑤‖²_2,Ω+12‖Δ𝛿𝐹‖²_2,Ω≥0,∀𝑡∈𝑅,𝛿𝑤∈𝐻²₀(Ω), (18) and by taking 𝜇=𝐹, condition 𝛽(𝜂,𝜉,𝜇)=0 yields ()(𝐿())−𝑡Δ𝐹,Δ𝛿𝐹=−𝑡𝑤,𝛿𝑤,𝐹²2(𝐿(𝛿𝑤,𝛿𝑤),𝐹). (19)
Substituting (19) into (18), dividing the obtained expression by 𝑡 and going to the limit for 𝑡→0, the following inequality is obtained: (Δ𝑤,Δ𝛿𝑤)−(𝐿(𝑤,𝐹),𝛿𝑤)≥𝑞,𝜇. (20)
Substituting 𝛿𝑤 by −𝛿𝑤 in (20), one obtains equality, that is, (17).
Let us use the following notation (Φ(𝑤,𝐹),𝜇)=𝑎₁(Δ𝑤,Δ𝜇)−(𝐿(𝑤,𝐹),𝜇)−(𝑞,𝜇).
Equation (17) can be given in the following form: Φ(𝑤,𝐹),𝜇=0, (21) and it is clear that Φ(𝑤,𝐹)∈𝐻⁻²(Ω).
Therefore, each point of the minimum of functional (17) on 𝑀 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 {𝑤,𝐹}∈𝑀. 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 minimizing 𝐽(𝑤,𝐹) on 𝑀 using the following scheme:(a)element 𝑤₀∈𝐻²₀(Ω) is taken arbitrarily;(b)after computation of 𝑤_𝑛, 𝐹_𝑛∈𝐻²₀(Ω) and 𝑤_𝑛+1∈𝐻²₀(Ω) is defined successively by solutions to the following problems: 𝛽𝑤_𝑛,𝐹_𝑛,𝜇=0,𝐹_𝑛∈𝐻²₀(Ω)∀𝜇∈𝐻²₀(Ω),Δ𝑤_𝑛+1,Δ𝜇=Δ𝑤_𝑛,Δ𝜇−𝜌_𝑛Φ𝑤_𝑛,𝐹_𝑛,𝜇∀𝜇∈𝐻²₀(Ω); (22)(23)(c)coefficient 𝜌_𝑛 on step (b) is defined by the condition 𝐽𝑤_𝑛+1,𝐹_𝑛+1−𝐽𝑤_𝑛,𝐹_𝑛≤𝜀Φ𝑤_𝑛,𝐹_𝑛,𝑤_𝑛+1−𝑤_𝑛,0<𝜀<1, (24)where 𝜀 stands for a parameter of the method.

Theorem 3. For the iteration process (22) to (25) (Φ(𝑤_𝑛,𝐹_𝑛),𝜇)→0 for 𝑛→0 an arbitrary initial point {𝑤₀,𝐹₀}∈𝑀, obtained through this procedure series {𝑤_𝑛,𝐹_𝑛} includes a subseries convergent to the weak solution to the problem ((6) and (4)).

Proof. A possibility of constructing the series {𝑤_𝑛,𝐹_𝑛} is yielded by an observation that for all 𝜌_𝑛𝑤_𝑛+1∈𝐻²₀(Ω) and, consequently, 𝐿(𝑤_𝑛+1,𝑤_𝑛+1)∈𝐻⁻²(Ω), ∇²_𝑘𝑤_𝑛+1∈𝐻⁻²(Ω) [19]. It means that the coupling equation 𝛽(𝑤_𝑛+1,𝐹_𝑛+1,𝜇)=0 is solvable. Consider the following difference: Δ𝐽_𝑛=𝐽𝑤_𝑛+1,𝐹_𝑛+1−𝐽𝑤_𝑛,𝐹_𝑛=12Δ𝑤_𝑛+1−𝑤_𝑛,Δ𝑤_𝑛+1+𝑤_𝑛+12Δ𝐹_𝑛+1−𝐹_𝑛,Δ𝐹_𝑛+1+𝐹_𝑛−𝑞,𝑤_𝑛+1−𝑤_𝑛. (25)
Owing to {𝑤_𝑛,𝐹_𝑛}∈𝑀, {𝑤_𝑛+1,𝐹_𝑛+1}∈𝑀, (25) gives Δ𝐽_𝑛=Φ𝑤_𝑛,𝐹_𝑛+1,𝛿𝑤2‖Δ𝛿𝑤‖²_0,Ω+12‖Δ𝛿𝐹‖²_0,Ω, (26) where 𝛿𝑤=𝑤_𝑛+1−𝑤_𝑛, 𝛿𝐹=𝐹_𝑛+1−𝐹_𝑛. Taking (23) into account, one observes that 𝛿𝑤 serves as a general solution to the boundary value problem: Δ²𝛿𝑤=−𝜌_𝑛Φ𝑤_𝑛,𝐹_𝑛,𝛿𝑤∈𝐻²₀(Ω). (27)
Further on it means that 𝛿𝑤=−𝜌_𝑛𝐺Φ𝑤_𝑛,𝐹_𝑛, (28) where 𝐺[•]∶𝐻⁻²(Ω)→𝐻²₀(Ω) stands for the linear bounded operator being inversed to operator Δ²(•). Therefore Δ𝐽_𝑛=−𝜌_𝑛Φ𝑤_𝑛,𝐹_𝑛,𝐺Φ𝑤_𝑛,𝐹_𝑛+12‖Δ𝛿𝑤‖²_0,Ω+12‖Δ𝛿𝐹‖²_0,Ω. (29)
Let us proceed to the second order terms. Taking in (11) 𝜇=𝛿𝑤 and applying (28), one gets ‖Δ𝛿𝑤‖²_0,Ω=−𝜌_𝑛Φ𝑤_𝑛,𝐹_𝑛,𝛿𝑤=𝜌²_𝑛Φ𝑤_𝑛,𝐹_𝑛,𝐺Φ𝑤_𝑛,𝐹_𝑛. (30)
Let us estimate the last term in (29). Since {𝑤_𝑛,𝐹_𝑛}∈𝑀 and {𝑤_𝑛+1,𝐹_𝑛+1}∈𝑀, then for 𝛿𝐹 the following equation should be satisfied: Δ𝛿𝐹,Δ𝜇+𝐿𝑤_𝑛,𝛿𝑤,𝜇+∇²_𝑘+1𝛿𝑤,𝜇2𝐿(𝛿𝑤,𝛿𝑤),𝜇=0𝛿𝐹∈𝐻²₀(Ω),∀𝜇∈𝐻²₀(Ω). (31)
This, in particular, yields ‖Δ𝛿𝐹‖_0,Ω≤𝑐₇𝐿𝑤_𝑛,𝛿𝑤_𝐿1(Ω)+‖𝐿(𝛿𝑤,𝛿𝑤)‖_𝐿¹_(Ω)+∇²_𝑘𝛿𝑤_𝐿1(Ω). (32)
However, 𝑤_𝑛 belongs to the bounded set in 𝐻²₀(Ω) for arbitrary 𝑛. It implies that ‖Δ𝛿𝐹‖_0,Ω≤𝑐₈‖Δ𝛿𝑤‖²_0,Ω or equivalently ‖Δ𝛿𝐹‖²_0,Ω≤𝑐₉𝜌⁴_𝑛Φ𝑤_𝑛,𝐹_𝑛,𝐺Φ𝑤_𝑛,𝐹_𝑛². (33)
Substituting (30), (33) into (29), and taking into account both positive defiantness (in the sense of (Φ(𝑤_𝑛,𝐹_𝑛),𝐺[Φ(𝑤_𝑛,𝐹_𝑛)])≥𝛼‖Φ(𝑤_𝑛,𝐹_𝑛)‖²) and the constraints of operator 𝐺[•], one gets Δ𝐽_𝑛≤−𝜌_𝑛𝑐₁₀Φ𝑤_𝑛,𝐹_𝑛²×𝜌−1+^𝑛2+𝑐₁𝜌³_𝑛2Φ𝑤_𝑛,𝐹_𝑛². (34)
The latter estimation shows that the values 𝜌_𝑛≠0 are responsible for the satisfaction of inequality (24). For this purpose 𝜌_𝑛 should be chosen in the following way: 𝜌_𝑛2+𝑐₁𝜌³_𝑛2Φ𝑤_𝑛,𝐹_𝑛²≤1−𝜀. (35)
It can always be done, since 0<𝜀<1.
Taking 𝜌_𝑛 in accordance with the algorithm applied so far, the following estimations are obtained on each step: Δ𝐽_𝑛≤−𝜌_𝑛𝜀Φ𝑤_𝑛,𝐹_𝑛², (36) which means that for the arbitrarily taken 𝑛 we have 𝐽_𝑛+1−𝐽_𝑛≤0. Since functional 𝐽 is bounded from below, the last inequality yields for 𝑛→∞Δ𝐽_𝑛→0. Besides, (36) gives Φ𝑤_𝑛,𝐹_𝑛²≤−Δ𝐽_𝑛𝜀𝜌_𝑛. (37)
Let us emphasize that the so far introduced algorithm of the choice of 𝜌_𝑛 guarantees that for arbitrary 𝑛 we have 𝜌_𝑛≥𝜌₀>0. In fact, because Δ𝐽_𝑛≤0, then 𝐽𝑤_𝑛,𝐹_𝑛≤𝐽𝑤₀,𝐹₀=𝐴. (38)
Owing to (38), norms ‖𝑤_𝑛‖_2,Ω,‖𝐹_𝑛‖_2,Ω are bounded. Therefore, also the norm ‖Φ(𝑤_𝑛,𝐹_𝑛)‖ is bounded. In addition, taking (37) into account, we have ‖Φ(𝑤_𝑛,𝐹_𝑛)‖→0 for 𝑛→∞, and consequently, also (Φ(𝑤_𝑛,𝐹_𝑛),𝜇)→0 for 𝑛→∞for all 𝜇∈𝐻²₀(Ω). The occurrence of a convergent subseries follows now from a bound of norms ‖𝑤_𝑛‖_2,Ω,‖𝐹_𝑛‖_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 Ω∈𝑅²: Δ²𝑤𝑥,𝑦=𝑔𝑥,𝑦. (39)

On the space boundary the following boundary conditions are given: 𝑤|_Γ=0,Δ𝑤−𝜒𝜕𝑤𝜕𝑛_Γ=0, (40) where 𝜒 denotes the curvature of boundary Γ.

Let us introduce the following new function 𝑀_𝑤=Δ𝑤. Substituting this function into (39) the following system of two Poisson-type equations is obtained: Δ𝑀_𝑤𝑥,𝑦=𝑔𝑥,𝑦,Δ𝑤𝑥,𝑦=𝑀_𝑤𝑥,𝑦,𝑥,𝑦∈Ω. (41)

Boundary conditions have the following form: 𝑤(𝑥,𝑦)_{𝑥,𝑦∈Γ}=0,𝑀_𝑤(𝑥,𝑦)_{𝑥,𝑦∈Γ}=0. (42)

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 𝜉={𝑤,𝑀_𝑤}: 𝐸=𝜉∈𝜌^∞(Ω)∣𝜉_Γ=0, (43) where 𝜌^∞(Ω) is the set of functions infinitely many times differentiable on Ω∈𝑅². Closure of set (43) in norm 𝐻²(Ω) is a subspace in 𝐻²(Ω) which is denoted by 𝑉(Ω). It is clear that 𝑉(Ω)=𝐻²(Ω)∩𝐻¹₀(Ω).

It is known (see [13]) that a solution to problems (39) and (40) is equivalent to minimization on 𝑉(Ω) of the following functional: 𝐽(𝑣)=12_ΩΔ𝜉²𝑑Ω−_Ω1𝑔𝜉𝑑Ω−2_Γ𝜒𝜕𝜉𝜕𝑛²𝑑𝑠. (44)

Hybrid Variation Problem Formulation. We assume that instead of functional (44) the following one is minimized: Φ=1𝜉,𝜓,𝛼2_Ω𝜓²𝑑Ω−_Ω1𝑔𝜉𝑑Ω−2_Γ𝜒|𝛼|²𝑑𝑠, (45) on such triads (𝜉,𝜓,𝛼)∈𝑉(Ω)×𝐿²(Ω)×𝐿²(Γ) that their elements are coupled through equalities −Δ𝜉=𝜓,(𝜕𝑤/𝜕𝑛)|_Γ=𝛼.

Let us define the space of the following functions: 𝑃(Ω)=𝜉,𝜓,𝛼∈𝐻¹₀(Ω)×𝐿²(Ω)×𝐿²(Γ)∣∀𝜇∈𝐻¹(Ω),𝛽𝜉,𝜓,𝛼,𝜇=0, (46) where the bilinear form 𝛽[⋅,⋅] is defined as 𝛽=1𝜉,𝜓,𝛼,𝜇2_Ω∇𝜉∇𝜇𝑑Ω−_Ω𝜓𝜇𝑑Ω−_Γ𝛼𝜇𝑑𝑠. (47)

Theorem 4. If the space Ω is convex and has a Lipschitz continuous boundary Γ, then, first, the map (𝜉,𝜓,𝛼)∈𝑃(Ω)→|𝜓|²_0,Ω, now and later |𝜉|_𝑚,Ω=(∑_|𝑘|=𝑚∫_Ω|𝜕^𝑘𝜉|²𝑑Ω)^1/2, ‖𝜉‖_𝑚,Ω=(∑_{|𝑘|≤𝑚}∫_Ω|𝜕^𝑘𝜉|²𝑑Ω)^1/2) is the norm on space 𝑃(Ω) equivalent to the real dot product form (𝜉,𝜓,𝛼)∈𝑃(Ω)→(|𝜉|²_1,Ω+|𝜓|²_0,Ω+|𝛼|²_0,Γ)^1/2 transforming 𝑃(Ω) into a Hilbert space; second, if (𝜉,𝜓,𝛼)∈𝑃(Ω), then 𝜉,𝜓,𝛼∈𝑉×𝐿²(Ω)×𝐿²(Γ),−Δ𝜉=𝜓,𝜕𝑤𝜕𝑛_Γ=𝛼, (48) and if (48) is satisfied, then (𝜉,𝜓,𝛼)∈𝑃(Ω).

Proof. In the beginning we show the second statement. Since Ω has a continuous boundary, then the following Green formula holds: _Ω∇𝜉∇𝜇𝑑Ω=−_ΩΔ𝜉𝜇𝑑Ω+_Γ𝜕𝜉𝜕𝑛𝜇𝑑𝑠,∀𝜉∈𝐻²(Ω),∀𝜇∈𝐻¹(Ω). (49) Let (𝜉,𝜓,𝛼)∈𝑃(Ω). Then 𝜉∈𝐻¹₀(Ω), 𝜓∈𝐿²(Ω), 𝛼∈𝐿²(Γ), and 𝛽[(𝜉,𝜓,𝛼),𝜇]=0 for all 𝜇∈𝐻¹(Ω).
The last condition, in particular for ∀𝜇∈𝐻¹₀(Ω) yields _Ω∇𝜉∇𝜇𝑑Ω=_Ω𝜓𝜇𝑑Ω. (50)
It follows from (50) that 𝑣 appears as a solution to the Dirichlet problem for the operator −Δ for 𝜉|_Γ=0. Since space Ω is convex, therefore 𝜉∈𝐻²(Ω) ([13], Section 7.1, page 373), and consequently 𝜉∈𝐻²(Ω)∩𝐻¹₀(Ω). Using now (50) for 𝜇∈𝐻¹₀(Ω), we get −Δ𝜉=𝜓. Using the Green formula for 𝜇∈𝐻¹(Ω), we find that (𝜕𝜉/𝜕𝑛)|_Γ=𝛼. Assume that (48) holds. We show that (𝜉,𝜓,𝛼)∈𝑃(Ω). Since 𝜉∈𝑉(Ω)⊂𝐻²(Ω) and −Δ𝜉=𝜓, (𝜕𝜉/𝜕𝑛)|_Γ=𝛼, then (49) yields ∫_Ω∇𝜉∇𝜇𝑑Ω=∫_Ω𝜓𝜇𝑑Ω+∫_Γ𝛼𝜇𝑑𝑠 for all 𝜇∈𝐻¹(Ω); that is, 𝛽[(𝜉,𝜓,𝛼),𝜇]=0 for all 𝜇∈𝐻¹(Ω). Besides 𝜉∈𝑉(Ω)⊂𝐻¹₀(Ω) and the second statement is proved.
Let us prove the first statement. Endowed with the multiplication norm 𝑃(Ω) is a Hilbert space. Let (𝜉,𝜓,𝛼)∈𝑃(Ω). Then, as it has been shown, 𝜉∈𝐻²(Ω)∩𝐻¹₀(Ω). For 𝜇=𝜉 condition 𝛽[(𝜉,𝜓,𝛼),𝜇]=0 yields 𝜉²_1,Ω≤𝐶₁𝜓_0,Ω𝜉_0,Ω. (51)
Let us introduce space 𝑀⊂𝐻¹(Ω) such that 𝐻¹(Ω)=𝐻¹₀(Ω)⊕𝑀. Besides, let us introduce the following operator 𝐵∶𝐻¹(Ω)→𝐿²(Ω) defined in the following way: for all 𝜓∈𝐿²(Ω)𝛼=𝐵𝜓∈𝐿²(Γ) we have a unique solution to the following equation: _Γ𝛼𝜇𝑑𝑠=_Ω∇𝜉∇𝜇𝑑Ω−_Ω𝜓𝜇𝑑Ω∀𝜇∈𝑀. (52)
Under the condition that 𝜉∈𝐻¹₀(Ω) satisfies the following: _Ω∇𝜉∇𝜇𝑑Ω=_Ω𝜓𝜇𝑑Ω𝜇∈𝐻¹₀(Ω). (53)
It is not difficult to verify that under the theorem conditions, 𝐵𝜓=−𝐵Δ𝜉=(𝜕𝜉/𝜕𝑛)|_Γ; that is, 𝐵 stands for the operator of external normal derivative for 𝜉∈𝐻²(Ω)∩𝐻¹₀(Ω). This operator is bounded; that is, 𝜉∈𝐻²(Ω)∩𝐻¹₀(Ω). We denote its norm by ‖𝐵‖. Then ‖𝐵‖=sup_𝑣∈𝐻²_(Ω)∩𝐻¹_0(Ω)(‖𝜕𝜉/𝜕𝑛‖_0,Γ/|Δ𝜉|_0,Ω), where ‖•‖_0,Γ denotes the norm associated with the scalar product (𝛼,𝛽)_𝐿2(Γ)=∫_Γ𝛼𝛽𝑑𝑠. Therefore, for all 𝜓∈𝐿²(Ω)‖𝛼‖_0,Γ≤‖𝐵‖|𝜓|_0,Ω. Hence, taking (52) into account we get (‖𝜉‖_1,Ω+‖𝜓‖_0,Ω+‖𝛼‖_0,Γ)≤C₂‖𝜓‖_0,Ω and the theorem is proved.

Results of this theorem allow us to transit from minimization of functional (45) on space 𝑉(Ω) to minimization of functional (45) on space 𝑃(Ω).

Theorem 5. Let 𝑤∈𝑉(Ω) be a solution to problem (44), then Φ𝑤,−Δ𝑤,𝜕𝑤𝜕𝑛⟶inf_{𝑣,𝜓,𝛼∈𝑃(Ω)}Φ𝜉,𝜓,𝛼. (54) In this case the triad (𝑤,−Δ𝑤,𝜕𝑤/𝜕𝑛)∈𝑃(Ω) is the unique solution to the problem of minimization of (54).

Proof. We prove that a symmetric bilinear form 𝑎([(𝜉,𝜓,𝛼),(𝜂,𝜑,𝛽)])=∫_Ω𝜓𝜑𝑑Ω−∫_Γ𝜒𝛼𝛽𝑑𝑠, (𝜉,𝜓,𝛼), (𝜂,𝜑,𝛽)∈𝑃(Ω) is continuous and elliptic on 𝑃(Ω).
Owing to Theorem 4, if (𝜉,𝜓,𝛼),(𝜂,𝜑,𝛽)∈𝑃(Ω), then −Δ𝜉=𝜓; (𝜕𝜉/𝜕𝑛)|_Γ=𝛼; −Δ𝜂=𝜑; (𝜕𝜂/𝜕𝑛)|_Γ=𝛽. Then we have 𝑎𝜉,𝜓,𝛼,𝜂,𝜙,𝛽=_ΩΔ𝜉Δ𝜂𝑑Ω−_Γ𝜕𝜉𝜕𝑛𝜕𝜂𝜕𝑛𝑑𝑠. (55)
For 𝜂=𝜉, 𝜑=𝜓, 𝛽=𝛼 from (55) we obtain ([13], Section 1.2, page 38) 𝑎,=𝜉,𝜓,𝛼𝜂,𝜓,𝛼_Ω𝜕²𝜉𝜕𝑥²²𝜕+2²𝜉𝜕𝑥²𝜕²𝜉𝜕𝑦²+𝜕²𝜉𝜕𝑥²²𝜕+2²𝜉𝜕𝑥𝜕𝑦²≥𝜓𝑑Ω_20,Ω. (56)
Then, elliptic property has been proved by 𝑃. Continuity of the bilinear form is evident. It means that the problem of minimization Φ𝜉^∗,𝜓^∗,𝛼^∗⟶inf_{(Ω)𝑢,𝜑,𝛽∈𝑃}Φ𝜂,𝜑,𝛽 (57) 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 (𝜉^∗,𝜓^∗,𝛼^∗)∈𝑃(Ω) is a solution to problem (57), then the following relations should hold: _Ω𝜓^∗𝜑𝑑Ω−_Ω𝜂𝑔𝑑Ω−_Γ𝜒𝛼^∗∀(Ω).𝛽𝑑𝑠=0𝜂,𝜑,𝛽∈𝑃 (58)
Since (𝜉^∗,𝜓^∗,𝛼^∗)∈𝑃(Ω), then −Δ𝜉^∗=𝜓^∗, (𝜕𝜉^∗/𝜕𝑛)|_Γ=𝛼^∗, and 𝜉^∗∈𝑃(Ω). Therefore, taking (58) into account we get ∫_ΩΔ𝜉^∗Δ𝜂𝑑Ω−∫_Γ𝜒(𝜕𝜉^∗/𝜕𝑛)(𝜕𝜂/𝜕𝑛)𝑑𝑠=∫_Γ𝜂𝑔𝑑𝑠. Therefore, 𝜉^∗ coincides with the solution 𝑤 to problems (39) and (40), and 𝜓^∗=−Δ𝑤, 𝛼^∗=(𝜕𝑤/𝜕𝑛)|_Γ.

Remark 6. Since the space is convex and its boundary is regular, then for 𝑔∈𝐻⁻¹(Ω) a solution to problems (39) and (40), 𝑤∈𝐻³(Ω)∩𝐻¹₀(Ω),Δ𝑤∈𝐻¹(Ω). (59)

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 transformation 𝐴∶𝐿²(Ω)→𝐻¹₀(Ω) in the following way: if 𝜓∈𝐿²(Ω) is a given function, then the function 𝜉=𝐴𝜓∈𝐻¹₀(Ω) is a unique solution to the equation ∫_Ω∇𝜉∇𝜇𝑑Ω=∫_Ω𝜓𝜇𝑑Ω∀𝜇∈𝐻¹₀(Ω), for 𝜉∈𝐻¹₀(Ω). This means that space 𝑃(Ω), defined by (43), can be presented in the following form: 𝑃(Ω)=𝜉,𝜓,𝛼∈𝐻¹₀(Ω)×𝐿²(Ω)×𝐿²(Γ)∣𝜉=𝐴𝜓,𝛼=𝐵𝜓. (60)

Problem (44) is equivalent to the following problem of optimal control: min_{𝜓∈𝐿2(Ω)}12_Ω𝜓²𝑑Ω−_Ω1𝑔𝜉𝑑Ω−2_Γ𝜒|𝛼|²𝑑𝑠, (61) where the state function 𝜉 and 𝛼 are coupled via control 𝜓∈𝐿²(Ω) through the following state equations: 𝜉∈𝐻¹₀(Ω),𝜉=𝐴𝜓,𝛼∈𝐿²(Γ),𝛼=𝐵𝜓. (62)

As it follows from Remark 6, although the optimal control 𝜓 is sought on 𝐿²(Ω), its regularity is higher for 𝑔∈𝐻⁻¹(Ω)𝜓∈𝐻¹(Ω). In this case the following trace is defined: 𝜓|_Γ=𝜆,𝜆∈𝑀. Furthermore, besides (62), we require that ∫_Ω∇𝜓_𝜆∇𝜇𝑑Ω=∫_Ω𝑔𝜇𝑑Ω for all 𝜇∈𝐻¹₀(Ω) and 𝜓_𝜆−𝜆∈𝐻¹₀(Ω). Then, if 𝜉_𝜆=𝐴𝜓_𝜆 and 𝛼_𝜆=𝐵𝜓_𝜆, (61) implies that min_𝜓∈𝐿²_(Ω)Φ(𝜉,𝜓,𝛼)=min_𝜆∈𝑀𝐷(𝜆), where 𝐷(𝜆)1=−2_Ω𝜓_𝜆²𝑑Ω−_Γ𝜆𝛼_𝜆−1𝑑𝑠2_Γ𝜒𝛼_𝜆²𝑑𝑠∀𝜆∈𝑀. (63)

The fundamental idea consists now in the application of a gradient method to the problem of minimization (63).

Let us take 𝑀^′ as a dual space for space 𝑀, and let ⟨⋅,⋅⟩ denote the relation of duality between spaces 𝑀 and 𝑀^′. We denote by 𝐷^′∈𝑀^′ a derivative of the functional 𝐷(𝜆). Let us introduce a map 𝑆∶𝑀→𝐻¹(Ω) in the following way: for 𝜆∈𝑀𝜑^∘_𝜆=𝑆(𝜆) is a unique function from 𝐻¹(Ω) satisfying the condition ∇𝜙_Ω^∘_𝜆∇𝜇𝑑Ω=0,∀𝜇∈𝐻¹₀(Ω),_∘𝜙_𝜆−𝜆∈𝐻¹₀(Ω). (64)

Theorem 7. For an arbitrary 𝜆∈𝑀 defined in (63), the functional is differentiable and its derivative is defined by the relation ∀𝜇∈𝑀𝐷^′(𝜆),𝜇=𝜑_Ω^∘𝜑_𝜃^∘_𝜇𝑑Ω,𝑤ℎ𝑒𝑟𝑒𝜑^∘_𝜇=𝑆𝜇,𝜑^∘_𝜃=𝑆𝜃_𝜆,𝜃_𝜆=𝜆+𝜒𝛼_𝜆,𝜃_𝜆∈𝑀. (65)

Proof. Differentiating (63) yields 𝐷^′(𝜆),𝜇=−_Ω𝜓𝜙_𝜆^∘_𝜇𝑑Ω−_Γ𝛼_𝜆−𝜇𝑑Ω_Γ𝜆+𝜒𝛼𝛽_𝜆^∘_𝜇𝑑𝑠, (66) where 𝜑^∘_𝜇=𝑆(𝜇),𝛽^∘_𝜇=𝐵𝜑^∘_𝜇.
From (66) and taking (52) into account we get 𝐷^′(𝜆),𝜇=−_Ω∇𝑣_𝜆∇𝜙^∘_𝜇𝑑Ω−_Γ𝜆+𝜒𝛼𝛽_𝜆^∘_𝜇𝑑𝑠. (67)
First term in (67) is equal to zero due to (64). Let us introduce the function 𝜑^∘_𝜃=𝑆(𝜃_𝜆), where 𝜃_𝜆 is defined by (66) and let 𝑢^∘_𝜇=𝐴𝜑^∘_𝜇. Then, (66) implies ⟨𝐷^′(𝜆),𝜇⟩=∫𝜙_Ω^∘𝜙_𝜃^∘_𝜇𝑑Ω−∫_Γ∇𝜙^∘_𝜃∇𝑢^∘_𝜇𝑑Ω.
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 {𝜆_𝑛}^∞_𝑛=0 of functions 𝜆^𝑛∈𝑀 via the following iteration scheme: 𝜇∈𝑀𝜆^𝑛+1−𝜆^𝑛,𝜇_𝑀=−𝜌_𝑛𝐷^′𝜆^𝑛,𝜇, (68) where (⋅,⋅)_𝑀 is the arbitrary scalar product in space 𝑀, 𝜌 is the positive parameter, and 𝜆⁰ is the arbitrary function of 𝑀.
Therefore, one iteration (68) corresponds to successive solutions of the following problems: (a)find for a given function 𝜆^𝑛∈𝑀 a unique function 𝜓^𝑛∈𝐻¹(Ω), satisfying the following relations: 𝜓^𝑛−𝜆^𝑛∈𝐻¹₀(Ω),∀𝜇∈𝐻¹₀(Ω)_Ω∇𝜓^𝑛∇𝜇𝑑Ω=_Ω𝑔𝜇𝑑Ω; (69)(70)(b)find a function 𝜉^𝑛∈𝐻¹₀(Ω), satisfying the following relation: ∀𝜇∈𝐻¹₀(Ω)_Ω∇𝜉^𝑛∇𝜇𝑑Ω=_Ω𝜓^𝑛𝜇𝑑Ω; (71)(c)find a function 𝛼^𝑛∈𝐿²(Γ), satisfying the following relation: ∀𝜇∈𝑀_Γ𝛼^𝑛𝜇𝑑Ω=_Ω∇𝜉^𝑛∇𝜇𝑑Ω−_Ω𝜓^𝑛𝜇𝑑Ω; (72)(d)find a function 𝜆^𝑛+1∈𝑀, satisfying the following relation: ∀𝜇∈𝑀𝜆^𝑛+1−𝜆^𝑛,𝜇_𝑀=−𝜌𝜙_Ω^∘𝜙^𝑛_𝜃^∘_𝜇𝑑Ω; (73)where 𝜑^∘𝑛_𝜃=𝑆𝜃^𝑛_𝜆, 𝜃^𝑛_𝜆=𝜆^𝑛−𝜒𝐵𝜓^𝑛, and 𝜑^∘_𝜇=𝑆𝜇.

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 map 𝐶∶𝐻¹(Ω)→𝑀 in the following way: for any function 𝜓∈𝐻¹(Ω) function 𝐶𝜓∈𝑀 is unique satisfying the following condition: ∀𝜇∈𝑀𝐶𝜓,𝜇_𝑀=_Ω𝜓𝜙^∘_𝜇𝑑Ω. (74)

Let us take ‖𝐶‖=sup_𝜓∈𝐻¹_(Ω)(|𝐶𝜓|_𝜇/|𝜓|_0,Ω), where |⋅|_𝑀 denotes the norm associated with the scalar product (⋅,⋅)_𝑀. It is clear that this norm exists, since the map 𝜇∈𝑀→𝜑^∘_𝜇∈𝐻¹(Ω) is bounded.

Theorem 8. If parameter 𝜌 satisfies the following condition: 20<𝜌<‖𝐶‖²1−𝜒_𝐿∞(Ω)‖𝑆‖‖𝐵‖², (75) then the iteration process (69)–(73) is convergent in the sense that lim_𝑛→∞𝜓^𝑛=𝜓в𝐿²(Ω),lim_𝑛→∞𝜉^𝑛=𝜉в𝐻¹₀(Ω),lim_𝑛→∞𝛼^𝑛=𝛼в𝐿²(Γ). (76)

Proof. It is sufficient to show that lim_𝑛→∞𝜓^𝑛=0 in 𝐿²(Ω) in the particular case when 𝑔=0. If we use the definition (74) of map 𝐶, then the recurrent formula (73) gives 𝜆^𝑛+1=𝜆^𝑛−𝜌𝐶𝜑^∘𝑛_𝜃,where𝜑^∘𝑛_𝜃=𝑆𝜃^𝑛_𝜆,𝜃^𝑛_𝜆=𝜆^𝑛+𝜒𝐵𝜓^𝑛, (77) and therefore 𝜆^𝑛+12_𝑀=𝜆^𝑛2_𝑀−2𝜌𝐶𝜙^∘𝑛_𝜃,𝜆^𝑛_𝑀+𝜌²𝐶𝜙^∘𝑛_𝜃²_𝑀. (78)
Consider the term (𝐶𝜑^∘𝑛_𝜃,𝜆^𝑛)_𝑀: 𝐶𝜑^∘𝑛_𝜃,𝜆^𝑛_𝑀=𝜑_Ω^∘𝑛_𝜃𝜓^𝑛𝑑Ω=−_Γ𝜆^𝑛+𝜒𝐵𝜓^𝑛𝐵𝜓^𝑛=𝑑𝑠_Ω𝜓^𝑛2𝑑Ω−_Γ𝜒𝐵𝜓^𝑛2𝑑𝑠≥𝜓^𝑛2_0,Ω. (79)
Let us estimate the norm |𝜑^∘𝑛_𝜃|: 𝜑^∘𝑛_𝜃0,Ω=𝑆(𝜆^𝑛+𝜒𝐵𝜓^𝑛)_0,Ω≤𝜓^𝑛_0,Ω+‖𝑆‖𝜒_𝐿∞(Ω)‖𝐵‖𝜓^𝑛_0,Ω≤1+𝜒_𝐿∞(Ω)‖𝑆‖‖𝐵‖𝜓^𝑛_0,Ω. (80)
The latter inequalities and (78) imply the following estimation: 𝜆^𝑛+12_𝑀−𝜆^𝑛2_𝑀≤−𝜌2𝜂−‖𝐶‖²1+𝜒_𝐿∞(Ω)‖𝑆‖‖𝐵‖²𝜌𝜓^𝑛2_0,Ω. (81)
Hence, in particular, we get lim_𝑛→∞𝜓^𝑛_0,Ω=0, (82) if 𝜌 satisfies inequalities (75). Besides, we have lim_𝑛→∞𝜉^𝑛=lim_𝑛→∞𝐴𝜓^𝑛=0в𝐻¹₀(Ω),lim_𝑛→∞𝛼^𝑛=lim_𝑛→∞𝐵𝜓^𝑛=0в𝐿²(Γ), (83) which finishes the proof.

Since convergence of the considered method is guaranteed, any choice of subspace 𝑀, satisfying the condition 𝐻¹(Ω)=𝐻¹₀(Ω)⊕𝑀 and any choice of scalar product (⋅,⋅)_𝑀 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 𝜆^𝑛+1∈𝑀.

If the scalar product in 𝑀 is defined via the following formula: 𝜆,𝜇_𝑀=_Ω𝑆𝜆⋅𝑆𝜇𝑑Ω, (84) then as 𝜆^𝑛+1 any function from 𝑀 can be taken, assuming that the following condition is satisfied: 𝜆^𝑛+1_Γ=𝜆^𝑛−𝜌𝜃^𝑛_𝜆Γ. (85)

Equality (85) can be understood in the sense of trace equality on a boundary. In fact, if (85) is satisfied, then 𝜆^𝑛+1−𝜆^𝑛,𝜇_𝑀=−𝜌_Ω𝑆𝜃^𝑛⋅𝑆𝜇𝑑Ω=−𝜌𝜑_Ω^∘𝜑^𝑛_𝜃^∘_𝜇𝑑Ω, (86) which means that conditions of Theorem 7 are satisfied.

We may choose also the following scalar product: either𝜆,𝜇_𝑀=_Ω∇𝜆∇𝜇𝑑Ω.or𝜆,𝜇_𝑀=_Ω𝜆𝜇𝑑Ω. (87)

However, in the latter case one needs to compute gradient ⟨𝐷^′(𝜆),𝜇⟩ 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 Section 3.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 Section 3.2 can be applied to solve (2) using the iterative procedure, whose convergence has been proved in Section 3.1.

For this purpose new variables 𝑀_𝑤(𝑥,𝑦) and 𝑀_𝐹(𝑥,𝑦) are introduced 𝑀_𝑤𝑥,𝑦=Δ𝑤𝑥,𝑦,𝑀_𝐹𝑥,𝑦=Δ𝐹𝑥,𝑦. (88)

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: Δ𝑀_𝑤^(𝑘)=𝑞+𝐿𝑤^(𝑘−1),𝐹^(𝑘−1),Δ𝑤^(𝑘)=𝑀_𝑤^(𝑘),Δ𝑀_𝐹^(𝑘)=𝐿𝑤^(𝑘),𝑤^(𝑘),Δ𝐹^(𝑘)=𝑀_𝐹^(𝑘),𝑥,𝑦∈Ω. (89)

Boundary conditions (4) are transformed to the following form: 𝑤|_Γ=𝑀_𝑤Γ=𝐹|_Γ=𝑀_𝐹Γ=0. (90)

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) where 𝑇 stands for a certain operator defined on set 𝐷(𝑇) of the Hilbert space 𝐿₂(Ω); 𝑔(𝑥,𝑦) is the function given for two variables 𝑥 and 𝑦; 𝜔(𝑥,𝑦) is the function of these two variables being sought; Ω(𝑥,𝑦) is the space of changes of variables 𝑥 and 𝑦.

If Ω(𝑥,𝑦)=𝑋×𝑌 (𝑋 is the certain bounded set of variables 𝑥, 𝑌 is the bounded set of 𝑦), then a solution to (91) can be given in the following form: 𝜔_𝑁𝑥,𝑦=^𝑁_𝑖=1𝑢_𝑖(𝑥)𝑣_𝑖𝑦, (92) where functions 𝑢_𝑖(𝑥) and 𝑣_𝑖(𝑦) are defined by the following system of equations: _𝑋𝑇𝜔_𝑁−𝑔𝑢₁(𝑥)⋮𝑑𝑥=0,_𝑋𝑇𝜔_𝑁−𝑔𝑢_𝑁(𝑥)𝑑𝑥=0,_𝑌𝑇𝜔_𝑁−𝑔𝑣₁𝑦⋮𝑑𝑦=0,_𝑌𝑇𝜔_𝑁−𝑔𝑣_𝑁𝑦𝑑𝑦=0. (93)

It is found in the following way: we have a certain system composed of 𝑁 functions regarding one variable, for instance, 𝑢⁰₁(𝑥),𝑢⁰₂(𝑥),…,𝑢⁰_𝑁(𝑥), and then from the first 𝑁 equations of system (93) the system of 𝑁 functions 𝑣¹₁(𝑥),𝑣¹₂(𝑥),…,𝑣¹_𝑁(𝑥) is defined. Then, the so far obtained functions represent a new choice of the functions regarding the variable 𝑥−𝑢²₁(𝑥),𝑢²₂(𝑥),…,𝑢²_𝑁(𝑥), and the latter serves to get a new set of functions regarding variable 𝑦−𝑣³₁(𝑥),𝑣³₂(𝑥),…,𝑣³_𝑁(𝑥), 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 functions 𝑢_𝑖(𝑥) and 𝑣_𝑖(𝑦) on the 𝑘th step, which, for example, corresponds to the choice of functions 𝑣^𝑘₁(𝑦),𝑣^𝑘₂(𝑦),…,𝑣^𝑘_𝑁(𝑦), we define the function 𝜔^𝑘_𝑁=^𝑁_𝑖=1𝑢_𝑖^𝑘−1(𝑥)𝑣^𝑘_𝑖𝑦, (94) taken as the approximating solution of (91) obtained by MVI.

Remark 10. Here and further on, we shall take as operator 𝑇 a certain differential operator defined on set 𝐷(𝑇) of the Hilbert space 𝐿₂(Ω). Then, on each step system (93) shall be transformed to a system of ODEs which can be solved further.

Remark 11. We call function 𝜔_𝑁(𝑥,𝑦) the 𝑁th approximation to (91) if the number of series terms in (92) is equal to 𝑁.

Let us study the case of first approximation; that is, the following solution of (91) is sought: 𝜔₁𝑥,𝑦=𝑢(𝑥)𝑣𝑦, (95) where functions 𝑢(𝑥) and 𝑣(𝑦) are defined through the illustrated way from the following system of equations: 𝑇𝑢(_𝑋𝑥)⋅𝑣𝑦−𝑔𝑢(𝑥)𝑑𝑥=0,_𝑌𝑇𝑢(𝑥)⋅𝑣𝑦−𝑔𝑢𝑦𝑑𝑦=0. (96)

Let the operator 𝑇 in (91) be positive definite. Let us introduce the following notation: 𝐻_𝑇(𝑋×𝑌) is the energy space of the operator 𝑇; [⋅,⋅] is the scalar product of elements in 𝐻_𝑇; 𝜔₀ is the exact solution to (91).

Theorem 12. If 𝑇 is a positive definite operator with the space of action 𝐷(𝑇)⊂𝐻_𝑇, then the sequence of elements 𝛼_𝑘=‖𝜔^𝑘₁(𝑥,𝑦)−𝜔₀‖_𝐻𝑇 is monotonously decreasing; that is, for arbitrary 𝑖 and 𝑗 if 𝑖≥𝑗, then 𝜔^𝑖₁−𝜔_0𝐻𝑇≤𝜔^𝑗₁−𝜔_0𝐻𝑇. (97)

Proof. We consider a subset 𝑀¹₁ of the space 𝐻_𝑇 which has the following form: 𝑀¹₁=𝜔𝑥,𝑦∣𝜔𝑥,𝑦=𝑢⁰(𝑥)𝑣𝑦,𝑢⁰(𝑥)∈𝐻_𝑇(𝑥),𝑣𝑦∈𝐻_𝑇(𝑌). (98)
It is clear, that set 𝑀¹₁ represents a subspace of space 𝐻_𝑇(𝑋×𝑌) (generally, of infinite dimension). Therefore, one may define 𝜔₀ projection onto space 𝑀¹₁. As it is known that element 𝑢⁰(𝑥)𝑣^∗(𝑦)∈𝑀¹₁ stands for the projection of 𝜔₀ onto 𝑀¹₁ if the following condition is satisfied: 𝑢⁰(𝑥)𝑣^∗𝑦−𝜔₀,𝑢⁰(𝑥)𝑣𝑦_𝐻𝑇=0 (99) for arbitrary elements 𝑢₀(𝑥)𝑣(𝑦)∈𝑀¹₁. It is clear that if 𝑢⁰(𝑥)𝑣^∗(𝑦)∈𝑀¹₁, then (99) coincides with the first equation of system (97).
Since the element 𝑢⁰(𝑥)𝑣¹(𝑦) obtained through the first step of MVI is a projection of element 𝜔₀ onto the subspace 𝑀¹₁, hence the following inequality holds: 𝜔¹₁𝑥,𝑦−𝜔_0𝐻𝑇≤𝑢⁰(𝑥)𝑣𝑦−𝜔_0𝐻𝑇 (100) for arbitrary elements 𝑢⁰(𝑥)𝑣(𝑦)∈𝑀¹₁. An analogous construction allows us to get a similar inequality for the subspaces; that is, we have 𝑀¹₂=𝜔𝑥,𝑦∣𝜔𝑥,𝑦=𝑢⁰(𝑥)𝑣¹𝑦,𝑢(𝑥)∈𝐻_𝑇(𝑥),𝑣¹𝑦∈𝐻_𝑇(𝑌). (101)
In the case corresponding to the second MVI step, 𝑢²(𝑥)𝑣¹𝑦−𝜔_0𝐻𝑇≤𝑢(𝑥)𝑣¹𝑦−𝜔_0𝐻𝑇 (102) for arbitrary elements 𝑢(𝑥)𝑣¹(𝑦)∈𝑀¹₂. It follows from (100) and (102) that ‖𝑢²(𝑥)𝑣¹(𝑦)−𝜔₀‖_𝐻𝑇≤‖𝑢⁰(𝑥)𝑣¹(𝑦)−𝜔₀‖_𝐻𝑇. Considerations similar to those so far provided and obtained for the 𝑘th MVI step prove the theorem as well as inequality (97) with the help of induction.

Remark 13. Results of Theorem 12 are extended into the case of 𝑁th approximation, and therefore inequality (97) is given in the following form: 𝜔^𝑛_𝑁𝑥,𝑦−𝜔_0𝐻𝑇≤𝜔^𝑛_𝑁𝑥,𝑦−𝜔_0𝐻𝑇,𝑚≥𝑛. (103)

In order to prove the theorem we introduce the following lemma.

Lemma 14. Let each of elements of the basis system of space 𝐻_𝑇 have the following form: 𝜃_𝑖𝑥,𝑦=𝜑_𝑖(𝑥)𝜓_𝑖𝑦,∀^𝑖𝜑_𝑖(𝑥)∈𝐻_𝑇(𝑋),𝜑_𝑖𝑦∈𝐻_𝑇(𝑌). (104)

If for the initial MVI approximation one takes any component of a certain basis function 𝜃_𝑖, that is, 𝑢⁰(𝑥)≡𝜑_𝑖(𝑥), then for an arbitrary number 𝑘 of the MVI steps the following inequality holds: 𝜔^𝑘₁𝑥,𝑦−𝜔_0𝐻𝑇≤𝑐𝜑_𝑖(𝑥)𝜓_𝑖𝑦−𝜔_0𝐻𝑇, (105) where 𝑐 is the arbitrarily taken real number.

Proof. Since 𝑢⁰(𝑥)𝑣¹(𝑦)≡𝜑_𝑖(𝑥)𝑣¹(𝑦), then Theorem 12 yields 𝜔^𝑘₁𝑥,𝑦−𝜔_0𝐻𝑇≤𝜑_𝑖(𝑥)𝑣¹𝑦−𝜔_0𝐻𝑇≤𝑐𝜑_𝑖(𝑥)𝜓_𝑖𝑦−𝜔_0𝐻𝑇. (106)
On the basis of the given lemma we formulate one of the MVI convergence criterions. Initially, we identify space 𝐻_𝑇 with space 𝑊^0𝑚₂(Ω) which is generated through a closure regarding the norm ‖𝜔‖_𝑊^𝑚₂=_Ω^𝑚_𝑘=0(𝑘)𝜕^𝑘𝜔𝜕^𝑘1𝑥𝜕^𝑘2𝑦𝑑𝑥𝑑𝑦² (107) of a set of infinitely differentiable functions 𝐶^0∞(Ω) with a compact carrier inΩ.

Theorem 15. Let each of elements of the basis system of space 𝑊^0𝑚₂(𝑋×𝑌) have the following form: 𝜃_𝑖𝑥,𝑦=𝜑_𝑖(𝑥)𝜓_𝑖𝑦, (108) where {𝜑_𝑖(𝑥)} stands for the basis system in space 𝑊^0𝑚₂(𝑋) and {𝜓_𝑖(𝑦)} in space 𝑊^0𝑚₂(𝑌), and in order to get an arbitrary 𝑁th order approximation of the MVI we take components of the elements of the basis system {𝜃_𝑖(𝑥,𝑦)} as the initial functions. Then, for sufficiently large 𝑁 the MVI gives a unique approximate solution 𝜔_𝑁, and the sequence {𝜔_𝑁} is convergent with respect to the norm of space 𝑊^0𝑚₂(𝑋×𝑌) to the exact solution 𝜔₀ irrespectively of the number of steps 𝑘. This construction can be carried out for each 𝑁th approximation; that is, 𝜔^𝑘_𝑁−𝜔_00𝑊𝑚2⟶0,𝑁⟶∞. (109)

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 arbitrary 𝑘th step. Therefore, let us consider the 𝑁th approximation of problem (91), obtained on the first step. In a way similar to considerations regarding Theorem 12, one may show that each 𝜔¹_𝑁 stands for a projection of element 𝜔₀ onto the following subspace: 𝑀¹_𝑁=𝜔𝑥,𝑦∣𝜔𝑥,𝑦=^𝑁_𝑖=1𝑢⁰(𝑥)𝑣¹𝑦, (110) where 𝑢⁰_𝑖(𝑥) denotes 𝑁 the fixed elements from system {𝜑_𝑖(𝑥)}, and for arbitrary 𝑖 and 𝑣_𝑖(𝑦) they cover the whole space 𝑊^0𝑚₂(𝑌). Therefore, 𝜔¹_𝑁=𝑃_𝑁𝜔₀, where 𝑃_𝑁 stands for the operator of an orthogonal projection onto subspace 𝑀¹_𝑁 which is bounded. Since elements of the basis system 𝜃_𝑖(𝑥,𝑦) have form (108), it is limiting dense ([20], page 191) in 𝑊^0𝑚₂(𝑋×𝑌). 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 Theorem 12 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 𝑥,𝑦,𝑧, and hence the approximate solution of MVI can be sought in the following form: 𝜔𝑥,𝑦,𝑧=𝑢(𝑥)𝑣𝑦,𝑧. (111)

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 operator 𝑇 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) where 𝜔(𝑥,𝑦) is the normal plate deflection in point 𝑥,𝑦; 𝑔(𝑥,𝑦) is the intensity of the normal load; 𝐷=𝐸ℎ³[12(1−𝜈²)]⁻¹; 𝐸,𝜈 are the Young modulus and Poisson constant, respectively; 2ℎ is the plate thickness, 𝑦=0; 1,𝑥=0;1, where 𝑦=𝑦/𝑏, 𝑥=𝑥/𝑎; 𝑎 and 𝑏 are the plate dimensions; and 𝑔(𝑥,𝑦)=𝑔𝑎²𝑏²[𝐸ℎ⁴]⁻¹. Plate space in the plane 𝑥,𝑦 is denoted by Ω, and its contour is Γ. Below, we study two types of boundary conditions.

In the case of a simple support 𝑤=𝜕²𝑤/𝜕𝑛²|_Γ=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: 𝑊(𝑥,𝑦)=𝐴sin(𝜋𝑥)sin(𝜋𝑦), which satisfies the boundary conditions. Substituting this into (112) and applying the Bubnov procedure we obtain 𝐴=(4/𝜋⁶)𝑞, and for 𝑞=1 we have 𝐴=0.00416. The deflection function 𝑤(𝑥,𝑦)=∑_𝑚,𝑛𝐴_𝑚,𝑛sin(𝑚𝜋𝑥/𝑎)sin(𝑛𝜋𝑦/𝑏) is substituted into (20) then multiplied by sin(𝑚𝜋𝑥/𝑎), sin(𝑛𝜋𝑦/𝑏), and integrated regarding the plate surface. The following deflection value is obtained: 𝑤=𝑞𝑥,𝑦𝑥24𝐷⁴−2𝑎𝑥³+𝑎³𝑥−4𝑞𝑎⁴𝜋⁵𝐷sin𝑚𝜋𝑥𝑎×_𝑚,𝑛1𝑚⁵𝑎_𝑚𝑡ℎ𝛼_𝑚+22𝑐ℎ𝛼_𝑚𝑐ℎ2𝛼_𝑚𝑦𝑏−𝛼_𝑚2𝑐ℎ𝛼_𝑚2𝑦𝑏𝑠ℎ2𝛼_𝑚𝑦𝑏. (113)

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 (Table 1).

Table 1 

Solution methods.

In the second case the plate contours are clamped; that is, 𝑤=(𝜕𝑤/𝜕𝑛)|_Γ=0. Here computations are carried out in the first approximation, where owing to Remark 16, we take sin(𝜋𝑥) 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 line 𝑦=0.5 are given in Table 2. These results coincide with the conclusion of Theorem 12 regarding monotonous series {𝑎_𝑘} behavior. The exact value of deflection equal to 0.0138 is taken from monograph [20].

Table 2 

Deflection of the plate quarter.

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 Section 3 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: 𝜕𝑤=²𝑤𝜕𝑛²𝜕=𝐹=²𝐹𝜕𝑛²=0,𝑥=𝑦=0,𝑥=𝑦=1,𝜕𝑤=²𝑤𝜕𝑛²=𝐹=𝜕𝐹𝜕𝑛=0;𝑥=𝑦=0;𝑥=𝑦=1,𝑤=𝜕𝑤𝜕𝑛=𝐹=𝜕𝐹𝜕𝑛=0;𝑥=𝑦=0;𝑥=𝑦=1. (114)(115)(116)

For simplicity, we apply the MVI using its first approximation 𝑁=1. ODEs are reduced to AE (algebraic equations) through FDM with approximation 0(ℎ²), which is solved using the Gauss method. The interval of integration [0,1] was divided into 100 parts. Relation 𝑞[𝑤(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 Δ𝑞=10, where in order to accelerate convergence of the iterative procedure, 𝑤 and 𝐹 were taken from the previous step.

Figure 2 

𝑞 versus 𝑤.

Dependence of the deflection change in the plate center on the number of iterations is shown in Figure 3 (curve (1)). Boundary conditions correspond to the plate support on flexible noncompressed ribs (114). The remaining parameters are 𝑁=20, 𝑞=60, 𝜈=0.28, and 𝜀=10⁻³. We require 16 iterations to achieve a priori given accuracy.

Figure 3 

𝑤 versus 𝑛.

One may see from Figure 3 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 𝑤 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 𝐹, 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: 𝑤𝑤=^odd+𝑤_even2. (117)

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.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

This work has been partially supported by the National Science Foundation of Poland under Grant MAESTRO2 no. 2012/04/A/ST8/00738 for years 2010–2016.