Iterative Method for Solving the Second Boundary Value Problem for Biharmonic-Type Equation

Abstract

Solving boundary value problems (BVPs) for the fourth-order differential equations by the reduction of them to BVPs for the second-order equations with the aim to use the achievements for the latter ones attracts attention from many researchers. In this paper, using the technique developed by ourselves in recent works, we construct iterative method for the second BVP for biharmonic-type equation, which describes the deflection of a plate resting on a biparametric elastic foundation. The convergence rate of the method is established. The optimal value of the iterative parameter is found. Several numerical examples confirm the efficiency of the proposed method.

1. Introduction

Solving BVPs for the fourth-order differential equations by the reduction of them to BVPs for the second-order equations with the aim to use a lot of efficient algorithms for the latter ones attracts attention from many researchers. Namely, for the biharmonic equation with the Dirichlet boundary condition, there is intensively developed the iterative method, which leads the problem to two problems for the Poisson equation at each iteration (see e.g., [1–3]). Recently, Abramov and Ul’yanova [4] proposed an iterative method for the Dirichlet problem for the biharmonic-type equation, but the convergence of the method is not proved. In our previous works [5, 6] with the help of boundary or mixed boundary-domain operators appropriately introduced, we constructed iterative methods for biharmonic and biharmonic-type equations associated with the Dirichlet boundary condition. For the biharmonic-type equation with Neumann boundary conditions in [7] an iterative method also was proposed. It is proved that the iterative methods are convergent with the rate of geometric progression. In this paper we develop our technique in the above-mentioned papers for the second BVP for the biharmonic-type equation. Namely, we consider the following problem Δ²𝑢−𝑢𝑎Δ=𝑢𝑔+₁𝑏𝑢=𝑓inΩ,(1.1)Δ𝑢=𝑔₂onΓ,(1.2),onΓ,(1.3)

where Δ is the Laplace operator, Ω is a bounded domain in 𝑅^𝑛(𝑛≥2),Γ is the sufficiently smooth boundary of Ω, and 𝑎,𝑏 are nonnegative constants. This problem has not yet considered in [4].

It should be noticed that when 𝑎=0,𝑏=0 (1.1) is the equation for a thin plate, and problem (1.1)–(1.3) is decomposed into two consecutive problems for the Poisson equations.

In this paper we suppose that 𝑎≥0,𝑏≥0,𝑎+𝑏>0. Then (1.1) describes the deflection of a plate resting on biparametric elastic foundation. For solving this equation several methods such as the boundary element, the finite element methods [8–10], domain/boundary element technique [11], and boundary differential integral equation (BDIE) method [12] were used. It should be noticed that at present the boundary element method is intensively developed and is used for solving more complex problems of plates and shells (see e.g., [13–15]).

In this paper we use completely different approach to (1.1). Two cases will be treated in dependence on the sign of 𝑎²−4𝑏. In the case if 𝑎²−4𝑏≥0 we can immediately decompose the problem into two problems for second-order equations. In the opposite case we propose an iterative method for reducing problem (1.1)–(1.3) to a sequence of second-order problems. The convergence of the method is established and verified on examples.

2. Case 𝑎²−4𝑏≥0

In this case we always can lead the original problem (1.1)–(1.3) to two problems for second-order equations. To do this, we put 𝜇=12𝑎+√𝑎²−4𝑏.(2.1) Then problem (1.1)–(1.3) is reduced to the following problems: 𝜇𝑣Δ=𝑣1−𝑏𝑣=𝑓,inΩ,𝜇𝑔²1−𝑔₁,onΓ,𝜇Δ𝑢𝑢=−𝑔𝑢₁=𝑣,inΩ,,onΓ.(2.2) These Dirichlet problems can be solved by known methods such as finite element method, boundary element method, or finite difference method. Some fast Poisson solvers in [16, 17] can be applied for the above problems.

3. Case 𝑎²−4𝑏<0

This case is very important in mechanics because (1.1) describes the bending plate on elastic foundation (see [18]).

3.1. Construction of Iterative Method on Continuous Level

As in [6], we set 𝜑Δ=𝑢=−𝑏𝑣𝑢,.((33..12))

Then problem (1.1)–(1.3) leads to the following second-order problems Δ𝑣−𝑣𝑎=𝑣=𝑔₂𝑓+𝜑,inΩ,Δ𝑢𝑢==𝑔₁𝑣,,oinnΩΓ,,,onΓ,(3.3) where all the functions 𝑢,𝑣, and 𝜑 are unknown but they are related with each other by (3.2).

Now consider the following iterative process for finding 𝜑 and simultaneously for finding 𝑣,𝑢. (1)Given 𝜑⁽⁰⁾∈𝐿₂(Ω), for example, 𝜑⁽⁰⁾=0inΩ.(3.4)(2)Knowing 𝜑^(𝑘)(𝑥) on Ω(𝑘=0,1,…) solve consecutively two problems Δ𝑣^(𝑘)−𝑣𝑎𝑣^(𝑘)=𝑓+𝜑^(𝑘)inΩ,_(𝑘)Δ𝑢_(𝑘)=𝑔₂onΓ,(3.5)𝑢=𝑣^(𝑘)inΩ,_(𝑘)=𝑔₁onΓ.(3.6)(3)Compute the new approximation 𝜑^(𝑘+1)=1−𝜏_𝑘+1𝜑^(𝑘)−𝑏𝜏_𝑘+1𝑢^(𝑘)inΩ,(3.7) where 𝜏_𝑘+1 is an iterative parameter to be chosen later.

3.2. Investigation of Convergence

In order to investigate the convergence of the iterative process (3.4)–(3.7), firstly we rewrite (3.7) in the canonical form of two-layer iterative scheme [19] 𝜑^(𝑘+1)−𝜑^(𝑘)𝜏_𝑘+1+𝜑^(𝑘)+𝑏𝑢^(𝑘)=0.(3.8)

Now, we introduce the operator 𝐴 defined by the formula 𝐴𝜑=𝑢,(3.9) where 𝑢 is the function determined from the problems Δ𝑣Δ𝑣𝑢−𝑢==𝑎=𝑣00𝑣=oo𝜑innnΓΩΓin.,,Ω,((33..1101))

The properties of the operator 𝐴 will be investigated in the sequel.

Now, let us return to the problem (3.3). We represent their solutions in the form 𝑢=𝑢₁+𝑢₂,𝑣=𝑣₁+𝑣₂,(3.12) where 𝑢₁,𝑢₂,𝑣₁,𝑣₂ are the solutions of the problems Δ𝑣₁𝑣−𝑎𝑣₁=𝜑inΩ,₁Δ𝑢₁=0onΓ,𝑢=𝑣₁inΩ,₁Δ𝑣₂=0onΓ,(3.13)𝑣−𝑎𝑣₂=𝑓inΩ,₂Δ𝑢₂=𝑔₂onΓ,𝑢=𝑣₂inΩ,₂=𝑔₁onΓ.(3.14) According to the definition of 𝐴 we have 𝐴𝜑=𝑢₁.(3.15) Since 𝜑 should satisfy the relation (3.2), taking into account the representation (3.12) we obtain the equation (𝐼+𝑏𝐴)𝜑=𝐹,(3.16) where 𝐼 is the identity operator, and 𝐹=−𝑏𝑢₂.(3.17)

Thus, we have reduced the original problem (1.1)–(1.3) to the operator equation (3.16), whose right-hand side is completely defined by the data functions 𝑓,𝑔, and ℎ, and coefficients 𝑎,𝑏.

Proposition 3.1 3.1. The iterative process (3.4)–(3.7) is a realization of the two-layer iterative scheme 𝜑^(𝑘+1)−𝜑^(𝑘)𝜏_𝑘+1+(𝐼+𝑏𝐴)𝜑^(𝑘)=𝐹,𝑘=0,1,2,…(3.18) for the operator equation (3.16).

Proof. Indeed, if in (3.5), (3.6) we put 𝑢^(𝑘)=𝑢^1(𝑘)+𝑢₂,𝑣^(𝑘)=𝑣^1(𝑘)+𝑣₂,(3.19) where 𝑣₂,𝑢₂ are the solutions of problem (3.14), then we get Δ𝑣^1(𝑘)𝑣−𝑎𝑣^1(𝑘)=𝜑^(𝑘)inΩ,_1(𝑘)Δ𝑢_1(𝑘)=0onΓ,𝑢=𝑣^1(𝑘)inΩ,_1(𝑘)=0onΓ.(3.20) From here it is easy to see that 𝐴𝜑^(𝑘)=𝑢^1(𝑘).(3.21) Therefore, taking into account the first relation of (3.19), the above equality, and (3.17), from (3.8) we obtain (3.18). Thus, the proposition is proved.

Proposition 3.1 enables us to lead the investigation of convergence of processs (3.4)–(3.7) to the study of scheme (3.18). For this reason we need some properties of the operator 𝐴.

Proposition 3.2. The operator 𝐴 defined by (3.9)–(3.11) is linear, symmetric, positive, and compact operator in the space 𝐿₂(Ω).

Proof. The linearity of 𝐴 is obvious. To establish the other properties of 𝐴 let us consider the inner product (𝐴𝜑,𝜑) for two arbitrary functions 𝜑 and 𝜑 in 𝐿₂(Ω). Recall that the operator 𝐴 is defined by (3.9)–(3.11). We denote by 𝑢 and 𝑣 the solutions of (3.10) and (3.11), where instead of 𝜑 there stands 𝜑. We have 𝐴𝜑,𝜑==_Ω𝑢𝜑𝑑𝑥=_Ω𝑢Δ𝑣−𝑎𝑣𝑑𝑥_Ω𝑢Δ𝑣𝑑𝑥−𝑎_Ω𝑢𝑣𝑑𝑥.(3.22) Applying the Green formula for the functions 𝑢 and 𝑣, vanishing on the boundary Γ, we obtain _Ω𝑢Δ𝑣𝑑𝑥=_Ω𝑣Δ𝑢𝑑𝑥=_Ω𝑣𝑣𝑑𝑥,_Ω𝑢𝑣=_Ω𝑢Δ𝑢=−_Ω∇𝑢⋅∇𝑢𝑑𝑥.(3.23) Hence, 𝐴𝜑,𝜑=_Ω𝑣𝑣𝑑𝑥+𝑎_Ω∇𝑢⋅∇𝑢𝑑𝑥.(3.24) Clearly, (𝐴𝜑,𝜑)=𝐴𝜑_Ω,𝜑=𝐴𝜑,𝜑,𝑣²𝑑𝑥+𝑎_Ω||∇𝑢||²𝑑𝑥≥0(3.25) are due to 𝑎≥0. Moreover, it is easy seen that (𝐴𝜑,𝜑)=0 if and only if 𝜑=0. So, we have shown that the operator 𝐴 is symmetric and positive in 𝐿₂(Ω).
It remains to show the compactness of 𝐴. As is well known that if 𝜑∈𝐿₂(Ω) then problem (3.10) has a unique solution 𝑣∈𝐻²(Ω), therefore, problem (3.11) has a unique solution 𝑣∈𝐻⁴(Ω). Consequently, the operator 𝐴 maps 𝐿₂(Ω) into 𝐻⁴(Ω). In view of the compact imbedding 𝐻⁴(Ω) into 𝐿₂(Ω) it follows that 𝐴 is compact operator in 𝐿₂(Ω).
Thus, the proof of Proposition 3.2 is complete.

Due to the above proposition the operator 𝐵=𝐼+𝑏𝐴(3.26) is linear, symmetric, positive definite, and bounded operator in the space 𝐿₂(Ω). More precisely, we have 𝛾₁𝐼<𝐵≤𝛾₂𝐼,(3.27) where 𝛾₁=1,𝛾₂=1+𝑏‖𝐴‖.(3.28) Notice that since the operator 𝐴 is defined by (3.9)–(3.11) its norm ‖𝐴‖ depends on the value of 𝑎 but not on 𝑏 in (1.1).

From the theory of elliptic problems [20] we have the following estimates for the functions 𝑣,𝑢 given by (3.10), (3.11): ‖𝑣‖_𝐻²_(Ω)≤𝐶₁‖𝜑‖_𝐿2(Ω),‖𝑢‖_𝐻⁴_(Ω)≤𝐶₂‖𝑣‖_𝐻²_(Ω),(3.29) where 𝐶₁,𝐶₂ are constants. Therefore, ‖𝑢‖_𝐻⁴_(Ω)≤𝐶₁𝐶₂‖𝜑‖_𝐿2(Ω).(3.30) Before stating the result of convergence of the iterative process (3.5)–(3.7) we assume the following regularity of the data of the original problem (1.1)–(1.3): 𝑓∈𝐿₂(Ω),𝑔₁∈𝐻^7/2(Γ),𝑔₂∈𝐻^5/2(Γ).(3.31) Then the problem (1.1)–(1.3) has a unique solution 𝑢∈𝐻⁴(Ω). For the function 𝑢₂ determined by (3.14) we have also 𝑢₂∈𝐻⁴(Ω).

Theorem 3.3. Let 𝑢 be the solution of problem (1.1)–(1.3) and 𝜑 be the solution of (3.16). Then, if {𝜏_𝑘+1} is the Chebyshev collection of parameters, constructed by the bounds 𝛾₁ and 𝛾₂ of the operator 𝐵, namely 𝜏₀=2𝛾₁+𝛾₂,𝜏_𝑘=𝜏⁰𝜌₀𝑡_𝑘+1,𝑡^𝑘=cos2𝑘−1𝜌2𝑀𝜋,𝑘=1,…,𝑀,₀=1−𝜉1+𝜉,𝜉=𝛾𝛾¹₂,(3.32) we have ‖‖𝑢^(𝑀)−𝑢‖‖_𝐻4(Ω)≤𝐶𝑞_𝑀,(3.33) where 𝐶=𝐶₁𝐶₂‖‖𝜑⁽⁰⁾−𝜑‖‖_𝐿2(Ω),(3.34) with 𝐶₁,𝐶₂ being the constant in (3.30), 𝑞_𝑀=2𝜌^𝑀11+𝜌^12𝑀,𝜌₁=1−√𝜉1+√𝜉.(3.35) In the case of the stationary iterative process, 𝜏_𝑘=𝜏₀(𝑘=1,2,….) we have ‖‖𝑢^(𝑘)−𝑢‖‖_𝐻4(Ω)≤𝐶𝜌^𝑘0,𝑘=1,2,….(3.36)

Proof. According to the general theory of the two-layer iterative schemes (see [21]) for the operator equation (3.16) we have ‖‖𝜑^(𝑀)−𝜑‖‖_𝐿2(Ω)≤𝑞_𝑀‖‖𝜑⁽⁰⁾−𝜑‖‖_𝐿2(Ω),(3.37) if the parameter {𝜏_𝑘+1} is chosen by the formulae (3.32) and ‖‖𝜑^(𝑘)−𝜑^∗‖‖_𝐻4(Ω)≤𝜌^𝑘0‖‖𝜑⁽⁰⁾−𝜑^∗‖‖_𝐿2(Ω),𝑘=1,2,…(3.38) if 𝜏_𝑘=𝜏₀(𝑘=1,2,…). In view of these estimates the corresponding estimates (3.33) and (3.36) follow from (3.30) applied to the problems Δ𝑣^(𝑘)−𝑣−𝑣𝑎𝑣^(𝑘)−𝑣=𝜑^(𝑘)−𝜑,inΩ,_(𝑘)Δ𝑢−𝑣=0onΓ,_(𝑘)𝑢−𝑢=𝑣^(𝑘)−𝑣,inΩ,_(𝑘)−𝑢=0onΓ,(3.39) which are obtained in the result of the subtraction of (3.3) from (3.5) and (3.6), respectively. The theorem is proved.

Remark 3.4. From the above theorem it is easy to see that for each fixed value of 𝑎 the numbers 𝜌₀ and 𝑞_𝑀 characterizing the rate of convergence of the iterative method decrease with the decrease of 𝑏. So, the smaller 𝑏 is, the faster the iterative process converges. In the special case when 𝑏=0 the mentioned above numbers also are zero, hence the iterative process converges at once and the original problem (1.1)–(1.3) is decomposed into two second-order problems.

3.3. Computation of the Norm ‖𝐴‖

As we see in Theorem 3.3 for determining the iterative parameter 𝜏 we need the bounds 𝛾₁ and 𝛾₂ of the operator 𝐵, and in its turn the latter bound requires to compute ‖𝐴‖. Therefore, below we consider the problem of computation ‖𝐴‖.

Suppose the domain Ω=[0,𝑙₁]×[0,𝑙₂] in the plane 𝑥𝑂𝑦. In this case by Fourier method we found that the system of functions 𝑒_𝑚𝑛(𝑥,𝑦)=2√𝑙₁𝑙₂sin𝑚𝜋𝑥𝑙₁sin𝑛𝜋𝑦𝑙₂(𝑚,𝑛=1,2,…)(3.40) is the eigenfunctions of the spectral problem Δ𝑢𝑢==0𝜆𝑢oninΓΩ,(3.41) corresponding to the eigenvalues 𝜆_𝑚𝑛=−𝜋²𝑚𝑙²²¹+𝑛²𝑙²².(3.42) Moreover, this system is orthogonal and complete in 𝐿₂(Ω).

Now let a function 𝜑∈𝐿₂(Ω) have the expansion 𝜑(𝑥,𝑦)=^∞_𝑚_,𝑛=1𝜑_𝑚𝑛𝑒_𝑚𝑛(𝑥,𝑦).(3.43) Then we have ‖𝜑‖²=(𝜑,𝜑)=^∞_𝑚_,𝑛=1||𝜑_𝑚𝑛||².(3.44) Representing the solution 𝑣 of (3.10) in the form of the double series 𝑣(𝑥,𝑦)=^∞_𝑚_,𝑛=1𝑣_𝑚𝑛𝑒_𝑚𝑛(𝑥,𝑦)(3.45) we found 𝑣_𝑚𝑛=𝜑^𝑚𝑛𝜆_𝑚𝑛−𝑎.(3.46) Next, we seek the solution of (3.11) in the form 𝑢(𝑥,𝑦)=^∞_𝑚_,𝑛=1𝑢_𝑚𝑛𝑒_𝑚𝑛(𝑥,𝑦).(3.47) Then from (3.45) we find 𝑢_𝑚𝑛=𝑣^𝑚𝑛𝜆_𝑚𝑛.(3.48) From the definition of the operator ‖𝐴‖ by (3.9)–(3.11) we have (𝐴𝜑,𝜑)=(𝑢,𝜑)=^∞_𝑚_,𝑛=1𝑢_𝑚𝑛𝑒_𝑚𝑛,^∞_𝑚_,𝑛=1𝜑_𝑚𝑛𝑒_𝑚𝑛=_∞_𝑚_,𝑛=1𝑢_𝑚𝑛𝜑_𝑚𝑛=^∞_𝑚_,𝑛=1||𝜑_𝑚𝑛||²=𝜆_𝑚𝑛𝜆_𝑚𝑛−𝑎_∞_𝑚_,𝑛=1||𝜑_𝑚𝑛||²≤𝜋²𝑚²𝑙¹⁻²+𝑛²𝑙²⁻²𝜋²𝑚²𝑙¹⁻²+𝑛²𝑙²⁻²+𝑎_∞_𝑚_,𝑛=1||𝜑_𝑚𝑛||²𝜋²𝑙¹⁻²+𝑙²⁻²𝜋²𝑙¹⁻²+𝑙²⁻²+𝑎(3.49) due to the orthogonality of the system {𝑒_𝑚𝑛} and (3.46), (3.48), and (3.42). Thus, there holds the estimate (𝐴𝜑,𝜑)≤1𝜋²𝑙¹⁻²+𝑙²⁻²𝜋²𝑙¹⁻²+𝑙²⁻²+𝑎(𝜑,𝜑).(3.50) The sign of equality occurs for 𝜑=𝑒₁₁(𝑥,𝑦). Since 𝐴 is shown to be symmetric and positive in 𝐿₂(Ω) it follows: ‖𝐴‖=sup_𝜑≠0(𝐴𝜑,𝜑)(𝜑,𝜑)=𝜋²𝑙¹⁻²+𝑙²⁻²𝜋1²𝑙¹⁻²+𝑙²⁻²+𝑎.(3.51)

4. Numerical Realization of the Iterative Method

In the previous section we proposed and investigated an iterative method for problem (1.1)–(1.3) in the case if 𝑎²−4𝑏<0. Now we study numerical realization of the method.

For simplicity we consider the case, where the domain Ω is a rectangle, Ω=[0,𝑙₁]×[0,𝑙₂] in the plane 𝑥𝑂𝑦. In this domain we construct an uniform grid Ω_ℎ=(𝑥,𝑦),𝑥=𝑖ℎ₁,𝑦=𝑗ℎ₂,0≤𝑖≤𝑁₁,0≤𝑗≤𝑁₂,(4.1) where ℎ₁=𝑙₁/𝑁₁,ℎ₂=𝑙₂/𝑁₂.

Denote by Ω_ℎ the set of inner nodes, by Γ_ℎ the set of boundary nodes of the grid, and by 𝑣_ℎ,𝑢_ℎ,… the grid functions defined on Ω_ℎ.

Now consider a discrete version of the iterative method (3.4)–(3.7) when 𝜏_𝑘≡𝜏. (1)Given a starting 𝜑^ℎ(0), for example, 𝜑^ℎ(0)=0inΩ_ℎ.(4.2)(2)Knowing 𝜑^ℎ(𝑘) on Ω_ℎ(𝑘=0,1,…) solve consecutively two difference problems Λ𝑣^ℎ(𝑘)−𝑣𝑎𝑣^ℎ(𝑘)_ℎ(𝑘)=𝑓_ℎ+𝜑^ℎ(𝑘)inΩ_ℎ,Λ𝑢_ℎ(𝑘)=𝑔_2ℎonΓ_ℎ,(4.3)𝑢_ℎ(𝑘)=𝑣^ℎ(𝑘)inΩ_ℎ,=𝑔_1ℎonΓ_ℎ,(4.4) where Λ is the discrete Laplace operator, (Λ𝑣)_𝑖𝑗=𝑣^{𝑖−1,𝑗}−2𝑣_𝑖𝑗+𝑣_𝑖+1,𝑗ℎ²¹+𝑣^{𝑖,𝑗−1}−2𝑣_𝑖𝑗+𝑣_𝑖,𝑗+1ℎ²².(4.5)(3)Compute the new approximation 𝜑^ℎ(𝑘+1)=1−𝜏_ℎ𝜑^ℎ(𝑘)−𝑏𝜏_ℎ𝑢^ℎ(𝑘)inΩ_ℎ.(4.6)It is easy to see that the convergence of the above iterative method is related to the discrete version 𝐴_ℎ of the operator 𝐴, defined by the formula 𝐴_ℎ𝜑_ℎ=𝑢_ℎ,(4.7) where 𝑢_ℎ is determined from the difference problems Λ𝑣_ℎ𝑣−_ℎ𝑎𝑣_ℎ=𝜑_ℎinΩ_ℎ,Λ𝑢_ℎ=0onΓ_ℎ,𝑢_ℎ=𝑣_ℎinΩ_ℎ,=0onΓ_ℎ.(4.8) Using the results of the spectral problem for the discrete Laplace operator Λ (see [19]) we find the bounds of 𝐴_ℎ: 1𝛽₂𝛽₂+𝑎𝐼≤𝐴^ℎ≤1𝛽₁𝛽₁+𝑎𝐼,(4.9) where 𝛽₁=4ℎ²¹sin²𝜋ℎ₁2𝑙₁+4ℎ²²sin²𝜋ℎ₂2𝑙₂𝛽₂,=4ℎ²¹cos²𝜋ℎ₁2𝑙₁+4ℎ²²cos²𝜋ℎ₂2𝑙₂.(4.10) Therefore, for the operator 𝐵_ℎ, the discrete version of 𝐵, we obtain the estimate 𝛾^ℎ1𝐼≤𝐵_ℎ≤𝛾^ℎ2𝐼,(4.11) where 𝛾^ℎ1=1+𝑏𝛽₂𝛽₂+𝑎,𝛾^ℎ2=1+𝑏𝛽₁𝛽₁+𝑎.(4.12) Hence, we choose 𝜏_ℎ=2𝛾^ℎ1+𝛾^ℎ2,(4.13) which is the optimal parameter of the iterative process (4.2)–(4.6).

Now we study the deviation of 𝑢^ℎ(𝑘) from 𝑢^(𝑘) obtained by the iterative process (3.4)–(3.7). In the future for short we will write ‖⋅‖ instead of ‖⋅‖_∞.

Proposition 4.1. For any 𝑘=0,1,… there holds the estimate ‖‖𝜑^ℎ(𝑘)−𝜑^(𝑘)‖‖=𝑂ℎ²,‖‖𝑢^ℎ(𝑘)−𝑢^(𝑘)‖‖=𝑂ℎ²,(4.14) where ℎ²=ℎ²¹+ℎ²²,  𝑢^(𝑘),𝜑^(𝑘) are computed by the process (3.4)–(3.7) and 𝑢^ℎ(𝑘),𝜑^ℎ(𝑘) are computed by (4.2)–(4.6).

Proof. We shall prove this proposition by induction in 𝑘.
For 𝑘=0 we have ||𝜑^ℎ(0)−𝜑⁽⁰⁾||=0 and the second estimate in (4.14) is valid due to the the second-order accuracy of the iterative schemes (4.3) and (4.4) for the problems (3.5) and (3.6) (see [21]).
Now suppose (4.14) is valid for 𝑘−1≥0. We shall show that it also is valid for 𝑘. For this purpose we recall that 𝜑^(𝑘) is computed by the formula 𝜑^(𝑘)=(1−𝜏)𝜑^(𝑘−1)−𝑏𝜏𝑢^(𝑘−1)onΩ,(4.15) where 𝜏=2(2+𝑏‖𝐴‖)(4.16) and 𝜑^ℎ(𝑘) is computed by the formula 𝜑^ℎ(𝑘)=1−𝜏_ℎ𝜑^{ℎ(𝑘−1)}−𝑏𝜏_ℎ𝑢^{ℎ(𝑘−1)}onΩ_ℎ,(4.17)𝜏_ℎ being given by (4.13) and (4.12).
From (4.10)–(4.13), (4.16), and (3.51) it is possible to obtain the estimate 𝜏_ℎ=𝜏+𝑂ℎ².(4.18) Next, subtracting (4.15) from (4.17) and taking into account the above formula we get 𝜑^ℎ(𝑘)−𝜑^(𝑘)=(1−𝜏)𝜑^{ℎ(𝑘−1)}−𝜑^(𝑘−1)+𝜏𝑏𝑢^{ℎ(𝑘−1)}−𝑢^(𝑘−1)+𝑂ℎ².(4.19) By the assumptions of the induction ‖‖𝑢^{ℎ(𝑘−1)}−𝑢^(𝑘−1)‖‖=𝑂ℎ²,‖‖𝜑^{ℎ(𝑘−1)}−𝜑^(𝑘−1)‖‖=𝑂ℎ²(4.20) from (4.19) it follows ‖𝜑^ℎ(𝑘)−𝜑^(𝑘)‖=𝑂(ℎ²). Now, having in mind this estimate due to the second-order approximation of the difference operators in (4.3) and (4.4) we get the second estimate in (4.14). Thus, the proof of the proposition is complete.

In realization of the discrete iterative process (4.2)–(4.6) we shall stop iterations when ‖𝑢^ℎ(𝑘)−𝑢^{ℎ(𝑘−1)}‖<TOL, where TOL is a some given accuracy. Then for the total error of the discrete solution 𝑢^ℎ(𝑘) there there holds the following theorem.

Theorem 4.2. For the total error of the discrete solution 𝑢^ℎ(𝑘) from the exact solution 𝑢 of the original problem (1.1)–(1.3) there holds the estimate ‖‖𝑢^ℎ(𝑘)−𝑢‖‖≤TOL+𝐶₂ℎ²+𝐶₁𝜌^0𝑘−1,(4.21) where 𝐶₁,𝐶₂ are some constants and 𝜌₀ is the number in (3.36).

Proof. We have the following estimate ‖‖𝑢^ℎ(𝑘)−𝑢‖‖≤‖‖𝑢^ℎ(𝑘)−𝑢^{ℎ(𝑘−1)}‖‖+‖‖𝑢^{ℎ(𝑘−1)}−𝑢^(𝑘−1)‖‖+‖‖𝑢^(𝑘−1)−𝑢‖‖.(4.22) Since the space 𝐻⁴(Ω) is continuously embedded to the space 𝐶(Ω) (see [20]) from (3.36) we have ‖‖𝑢^(𝑘−1)−𝑢‖‖≤𝐶₁𝜌^0𝑘−1(4.23) for some constant 𝐶₁. Now using Proposition 4.1 and the above estimate, from (4.22) we obtain (4.21). Thus, the theorem is proved.

Below we report the results of some numerical examples for testing the convergence of the iterative method. In all examples we test the iterative method for some values of 𝑎 and 𝑏 with TOL=10⁻⁵. The results of convergence of the method are given in tables, where 𝑘 is the number of iterations, 𝐸 is the error of approximate solution 𝑢^ℎ(𝑘), 𝐸=||𝑢^ℎ(𝑘)−𝑢||.

Example 4.3. We take an exact solution 𝑢=sin𝑥sin𝑦 in the rectangle [0,𝜋]×[0,𝜋] and corresponding boundary conditions. The right-hand side of (1.1) in this case is 𝑓=(4+2𝑎+𝑏)sin𝑥sin𝑦.
The results of convergence in the case of the uniform grids including 65×65,129×129, and 257×257 nodes for 𝑎=0,𝑎=0.5 and 𝑎=1 are presented in Tables 1, 2, and 3, respectively.

Table 1 

Convergence of the method in Example 4.3 for 𝑎=0.

Table 2 

Convergence of the method in Example 4.3 for 𝑎=0.5.

Table 3 

Convergence of the method in Example 4.3 for 𝑎=1.

Example 4.4. We take an exact solution 𝑢=𝑥sin𝑦+𝑦sin𝑥 in the rectangle [−𝜋,𝜋]×[−𝜋,𝜋] with corresponding boundary conditions. The the right-hand side of (1.1) is 𝑓=(1+𝑎+𝑏)(𝑥sin𝑦+𝑦sin𝑥).
The results of convergence in the case of the grids including 65×65,129×129, and 257×257 nodes for 𝑎=0,𝑎=0.5, and 𝑎=1 are presented in Tables 4, 5, and 6, respectively.
In the two above examples the grid step sizes are ℎ=𝜋/64,𝜋/128,𝜋/256. Therefore, ℎ²=0.0024,6.0239𝑒−4,  1.5060𝑒−4. According to the estimate (4.21) the total error of the discrete approximate solution depends on ℎ². The columns 𝐸1,𝐸2,𝐸3 in Tables 1–6 show this fact. It is interesting to notice that in these tables 𝐸1/𝐸2,𝐸2/𝐸3≈4 and (64×64)/(128×128)=(128×128)/(256×256)=4. It means that if the number of grid nodes increases in 4 times then it is expected that the accuracy of the approximate solution increases in the same times. From the tables we also observe that the number of iterations increases with the increase of the parameter 𝑏 for fixed values of parameter 𝑎. This confirms Remark 3.4. We also remark that for each pair of the parameters 𝑎and𝑏 the number of iterations for achieving an accuracy corresponding to the grid step size (or density of grid) does not depend on the grid step size.
In general the error of discrete approximate solution strongly depends on the step size of grid. So, it is not expected to get an approximate solution of higher accuracy on a grid of low density. However, in some exceptional cases we can obtain very accurate approximate solution on sparse grid. Below is an example showing this fact.

Table 4 

Convergence of the method in Example 4.4 for 𝑎=0.

Table 5 

Convergence of the method in Example 4.4 for 𝑎=0.5.

Table 6 

Convergence of the method in Example 4.4 for 𝑎=1.

Example 4.5. We take an exact solution 𝑢=(𝑥²−4)(𝑦²−1) in the rectangle [−2,2]×[−1,1]. The the right-hand side of (1.1) is 𝑓=8−2𝑎(𝑥²+𝑦²−5)+𝑏(𝑥²−4)(𝑦²−1).
The results of convergence in the case of the grids including 65×65 and 129×129 nodes are presented in Tables 7 and 8, respectively.
From Tables 7 and 8 we see a high accuracy even on the grid 65×65. The reason of this fact is that for quadratic function the approximation error of the central difference scheme is zero for any grid step size.
The above three numerical examples demonstrate the fast convergence of the iterative method (3.4)–(3.7) for problem (1.1)–(1.3).
Below, we consider an example for examining the variation of the solution of Prob. (1.1)–(1.3) in dependence of the parameters 𝑎 and 𝑏.

Table 7 

Convergence of the method in Example 4.5 for grid 65×65.

Table 8 

Convergence of the method in Example 4.5 for grid 129×129.

Example 4.6. We take the the right-hand side of (1.1) 𝑓=1 and the boundary conditions 𝑔₁=𝑔₂=0 and use the proposed method for finding approximate solution for different values of 𝑎 and 𝑏. In all experiments we use the grid of 65×65 nodes in the domain Ω=[−2,2]×[−1,1] and TOL=10⁻⁵. It turns out that the number of iterations in all cases of 𝑎 and 𝑏 does not exceed 7 and the value of the solution at any fixed point is decreasing with the growth of 𝑎 and 𝑏. This fact is obvious from Figure 1. The graph of the solution in the case of 𝑎=0.5 and 𝑏=1 is given in Figure 2.

Figure 1 

Variation of the value of 𝑢 in the middle point with 𝑏 for several values of 𝑎.

Figure 2 

Solution in the case of 𝑎=0.5 and 𝑏=1.

5. Concluding Remark

In the paper an iterative method was proposed for reducing the second problem for biharmonic-type equation to a sequence of Dirichlet problems for second-order equations. The convergence of the method was proved. In the case when the computational domain is a rectangle the optimal iterative parameter was given. Several numerical examples in this case show fast convergence of the method. When the computational domain consists of rectangles the proposed iterative method can be applied successfully if combining with the domain decomposition method.

Acknowledgments

This work is supported by Vietnam National Foundation for Science and Technology Development (NAFOSTED) under the Grant 102.99–2011.24. The authors would like to thank the anonymous reviewers sincerely for their helpful comments and remarks to improve the original manuscript.