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

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.


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][6][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 2 The Scientific World Journal 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][11][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][14][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 q(x, y, t) y z x Figure 1: A studied plate. 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.

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 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: Equations governing the deflection ( , , ) and stress function ( , ) have the following form [15]: The following operators are introduced: Here and further on the following nondimensional quantities are introduced: The Scientific World Journal 3 = ( /ℎ)√ / ; = / , 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]: where Γ stands for the space boundary occupied by the plate.
The following initial conditions are attached to (2): 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).

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.

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: On the first step of the iterative procedure the following biharmonic equation for a given load ( , ) is solved: The value of (1) ( , ) is substituted into the right-hand side of equation system (6), and as a result a biharmonic equation for (1) ( , ) 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 2 (Ω) refer to a Sobolev space of functions = { , } such that where 2 (Ω) denotes the space of functions being summed with a square in Ω. Let 2 0 (Ω) denote the closure of functions from (Ω) (space of functions of class ∞ in Ω, having compact carrier in Ω) in norm 2 (Ω): Since space Ω is bounded and its boundary Γ is efficiently regular, then map → ‖Δ ‖ 0,Ω defines the norm in 2 0 (Ω) being equivalent to the norm generated by spaces 2 (Ω).
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: Points of the functional minimum (12) are weak solutions to problems (6) and (4).
Let us use the following notation (Φ( , ), ) = Equation (17) can be given in the following form: and it is clear that Φ( , ) ∈ −2 (Ω). 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.
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). 6 The Scientific World Journal

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 Ω ∈ 2 : On the space boundary the following boundary conditions are given: 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: Boundary conditions have the following form: 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.
Results of this theorem allow us to transit from minimization of functional (45) on space (Ω) to minimization of functional (45) on space (Ω).

Remark 6.
Since the space is convex and its boundary is regular, then for ∈ −1 (Ω) a solution to problems (39) and (40), 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 −Δ.
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 : → 1 (Ω) in the following way: for ∈ ∘ = ( ) is a unique function from 1 (Ω) satisfying the condition Theorem 7. For an arbitrary ∈ defined in (63), the functional is differentiable and its derivative is defined by the relation Proof. Differentiating (63) yields Ω. 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: where (⋅, ⋅) is the arbitrary scalar product in space , is the positive parameter, and 0 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 ∈ 1 (Ω), satisfying the following relations: (b) find a function ∈ 1 0 (Ω), satisfying the following relation: (c) find a function ∈ 2 (Γ), satisfying the following relation: (d) find a function +1 ∈ , satisfying the following relation: 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 : 1 (Ω) → in the following way: for any function ∈ 1 (Ω) function ∈ is unique satisfying the following condition: Let us take ‖ ‖ = sup ∈ 1 (Ω) (| | /| | 0,Ω ), where | ⋅ | denotes the norm associated with the scalar product (⋅, ⋅) . It is clear that this norm exists, since the map ∈ → ∘ ∈ 1 (Ω) is bounded. Theorem 8. If parameter satisfies the following condition: then the iteration process (69)-(73) is convergent in the sense that Proof. It is sufficient to show that lim → ∞ = 0 in 2 (Ω) in the particular case when = 0. If we use the definition (74) of map , then the recurrent formula (73) gives and therefore The Scientific World Journal 9 Consider the term ( ∘ , ) : Let us estimate the norm | ∘ |: The latter inequalities and (78) imply the following estimation: Hence, in particular, we get if satisfies inequalities (75). Besides, we have which finishes the proof.
Since convergence of the considered method is guaranteed, any choice of subspace , satisfying the condition 1 (Ω) = 1 0 (Ω) ⊕ 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: then as +1 any function from can be taken, assuming that the following condition is satisfied: Equality (85) can be understood in the sense of trace equality on a boundary. In fact, if (85) is satisfied, then which means that conditions of Theorem 7 are satisfied.
We may choose also the following scalar product: 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.

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. 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: Boundary conditions (4) are transformed to the following form: 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.

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: where stands for a certain operator defined on set ( ) of the Hilbert space 2 (Ω); ( , ) 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: where functions ( ) and V ( ) are defined by the following system of equations: It is found in the following way: we have a certain system composed of functions regarding one variable, for instance, 0 1 ( ), 0 2 ( ), . . . , 0 ( ), and then from the first equations of system (93) the system of functions V 1 1 ( ), V 1 2 ( ), . . . , V 1 ( ) is defined. Then, the so far obtained functions represent a new choice of the functions regarding the variable − 2 1 ( ), 2 2 ( ), . . . , 2 ( ), and the latter serves to get a new set of functions regarding variable − V 3 1 ( ), V 3 2 ( ), . . . , V 3 ( ), 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 V ( ) on the th step, which, for example, corresponds to the choice of functions 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 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 ( , ) 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: where functions ( ) and V( ) are defined through the illustrated way from the following system of equations: 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 ; 0 is the exact solution to (91).

Theorem 12. If is a positive definite operator with the space of action ( ) ⊂
, then the sequence of elements = ‖ 1 ( , ) − 0 ‖ is monotonously decreasing; that is, for arbitrary and if ≥ , then Proof. We consider a subset 1 1 of the space which has the following form: It is clear, that set 1 1 represents a subspace of space ( × ) (generally, of infinite dimension). Therefore, one may define 0 projection onto space 1 1 . As it is known that element 0 ( )V * ( ) ∈ 1 1 stands for the projection of 0 onto 1 1 if the following condition is satisfied: for arbitrary elements 0 ( )V( ) ∈ 1 1 . It is clear that if 0 ( )V * ( ) ∈ 1 1 , then (99) coincides with the first equation of system (97).
Since the element 0 ( )V 1 ( ) obtained through the first step of MVI is a projection of element 0 onto the subspace 1 1 , hence the following inequality holds: for arbitrary elements 0 ( )V( ) ∈ 1 1 . An analogous construction allows us to get a similar inequality for the subspaces; that is, we have In the case corresponding to the second MVI step, for arbitrary elements ( )V 1 ( ) ∈ 1 2 . It follows from (100) and (102) that ‖ 2 ( )V 1 ( ) − 0 ‖ ≤ ‖ 0 ( )V 1 ( ) − 0 ‖ . 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: In order to prove the theorem we introduce the following lemma.
If for the initial MVI approximation one takes any component of a certain basis function , that is, 0 ( ) ≡ ( ), then for an arbitrary number of the MVI steps the following inequality holds: where is the arbitrarily taken real number.
Proof. Since 0 ( )V 1 ( ) ≡ ( )V 1 ( ), then Theorem 12 yields On the basis of the given lemma we formulate one of the MVI convergence criterions. Initially, we identify space with space 0 2 (Ω) which is generated through a closure regarding the norm 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 1 stands for a projection of element 0 onto the following subspace: where 0 ( ) denotes the fixed elements from system { ( )}, and for arbitrary and V ( ) they cover the whole space 0 2 ( ). Therefore, 1 = 0 , where stands for the operator of an orthogonal projection onto subspace 1 which is bounded. Since elements of the basis system ( , ) have form (108), it is limiting dense ( [20], page 191) in 0 2 ( × ). 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: 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.

Numerical
Results. Let us consider the following rectangular plate where ( , ) is the normal plate deflection in point , ;  denoted by Ω, and its contour is Γ. Below, we study two types of boundary conditions. In the case of a simple support = 2 / 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: ( , ) = sin( ) sin( ), which satisfies the boundary conditions. Substituting this into (112) and applying the Bubnov procedure we obtain = (4/ 6 ) , 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: 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).
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]. (6). A simultaneous use of the iterative procedure is described in Section 3 and the method of (1) decrease of the order of the system twice (from the 8th to 4th order);

Iterative Procedure of Linearization and Variational Iteration for System
(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: For simplicity, we apply the MVI using its first approximation = 1. ODEs are reduced to AE (algebraic equations) through FDM with approximation 0(ℎ 2 ), 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.
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 −3 . We require 16 iterations to achieve a priori given accuracy.
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  (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: 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.

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

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.