Numerical solution of the static beam problem by Bernoulli collocation method

Accepted Manuscript

Numerical solution of the static beam problem by Bernoulli collocation method Quanwei Ren, Hongjiong Tian PII: DOI: Reference:

S0307-904X(16)30268-2 10.1016/j.apm.2016.05.018 APM 11171

To appear in:

Applied Mathematical Modelling

Received date: Revised date: Accepted date:

8 April 2015 16 March 2016 17 May 2016

Please cite this article as: Quanwei Ren, Hongjiong Tian, Numerical solution of the static beam problem by Bernoulli collocation method, Applied Mathematical Modelling (2016), doi: 10.1016/j.apm.2016.05.018

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights • We use Bernoulli polynomial approximation to reduce the problem to a finite dimensional system.

CR IP T

• We provide an iterative process to solve this system together with a collocation method.

AC

CE

PT

ED

M

AN US

• Numerical experiments verify that our numerical method has high accuracy and efficiency.

1

ACCEPTED MANUSCRIPT

Quanwei Rena , Hongjiong Tianb,∗ a

CR IP T

Numerical solution of the static beam problem by Bernoulli collocation method ✩

AN US

Department of Mathematics, Shanghai Normal University, 100 Guilin Road, Shanghai 200234, P. R. China b Department of Mathematics, Shanghai Normal University; Division of Computational Science, E-Institute of Shanghai Universities; Scientific Computing Key Laboratory of Shanghai Universities; 100 Guilin Road, Shanghai 200234, P. R. China

Abstract

ED

M

We propose a numerical scheme to obtain an approximate solution of a boundary value problem for fourth order nonlinear integro-differential equation of Kirchhoff type. We first reduce the problem to a nonlinear finite dimensional system based on the Bernoulli polynomial approximation and then solve it by an iterative process together with a collocation method. Convergence of the iterative process and error estimates of the approximate solution are provided. Numerical experiments are conducted to illustrate the performance of the proposed method.

PT

Keywords: Bernoulli polynomial, collocation method, nonlinear fourth order integro-differential equation, iterative process, error estimate. 1. Introduction

AC

CE

Integro-differential equation is an equation that both integral(s) and derivative(s) of the unknown function are involved and frequently arises in fluid dynamics, biology, chemical kinetics, economics, potential theory and many others (see [2, 10] and the references therein). ✩

The work of this author is supported in part by E-Institutes of Shanghai Municipal Education Commission (No. E03004), NSF of China (No. 11471217), and NSF of Shanghai (No. 16ZR1424900). ∗ Corresponding author.

Preprint submitted to Applied Mathematical Modelling

May 27, 2016

ACCEPTED MANUSCRIPT

CR IP T

In this paper, we consider the following nonlinear fourth order integro-differential equation of Kirchhoff type [11, 38] ( RL y (4) (t) − εy 00 (t) − L2 0 [y 0 (t)]2 dty 00 (t) = f (t), 0 < t < L, (1.1) y(0) = y(L) = y 00 (0) = y 00 (L) = 0,

M

AN US

where y(t) represents the static deflection of the beam, ε and L are positive constants, f (t) is a continuous function defined on the interval [0, L]. This equation models the bending equilibrium of an extensible beam (of length L) which is simply supported at t = 0 and attached to a fixed nonlinear torsional spring at t = L. f (t) stands for the force exerted by the foundation. The mathematical RL modelling for elastic beams with nonlocal terms of the type 0 [y 0 (t)]2 dty 00 (t) can be found in [3, 18]. Due to the presence of an integral over [0, L], the equation is not a pointwise identity and therefore is called a nonlocal problem. In 2010, Dang and Luan [9] proved that Eq. (1.1) has a unique classical solution under the assumptions that ε > 0 and f (t) is continuous and nonpositive (or nonnegative) on the interval [0, L]. Eq. (1.1) is a generalization of the stationary problem associated with Z L 2 ! 2 ∂y E ∂ 2 y EI ∂ 4 y H dx ∂ y = 0, + − + (1.2) 2 4 ∂t ρA ∂x ρ 2ρL 0 ∂x ∂x2

CE

PT

ED

which was proposed by Woinowsky-Krieger [38] as a model for the deflection of an extensible beam of length L with hinged ends, where H, E, ρ, I and A represent the tension at rest, the Young elasticity modulus, density, cross-sectional moment of inertia, and cross-sectional area, respectively. The nonlinear term in Eq. (1.2) is a correction to the classical Euler-Bernoulli equation ∂ 2 y EI ∂ 4 y + = 0, ∂t2 ρA ∂x4

AC

which does not consider the changes of the tension induced by the variation of the length during the deflection. This kind of correction was proposed by Kirchhoff [13] to generalize D’Alembert’s equation with clamped ends. An interesting discussion on the various models for beam equations was given by Arosio [3]. The mathematical analysis of Eq. (1.2) and its generalizations are referred to [4, 5, 7, 12, 14, 16, 17, 19, 25]. In the study of the problem (1.1), a significant difficulty lies in a nonlinear term under the integral sign and hence many numerical methods have been developed in the last decades. For Eq. (1.1) and its generalizations, as well as 3

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

for equations similar to (1.1), the construction of numerical algorithms is usually a combination of two approximate methods, one of which reduces the problem to a finite-dimensional algebraic system and the other is some iterative process for solution of the finite-dimensional one. For example, Shin [28] developed a finite element approximation method to approximate the problem (1.1) based on the space of piecewise linear polynomials on the uniform grid and proposed an iterative scheme to solve a system of nonlinear equations during direct implementation. Later, Ohm, Lee and Shin [21] obtained the error estimates for this finite element approximation. In 2009, Peradze [24] used a Galerkin approximation to reduce the problem (1.1) to a system of cubic algebraic equations which are solved by means of the nonlinear Jacobi iteration process and estimated the algorithm error. Dang and Luan [9] reduced the problem to two second order equations and a nonlinear equation and then solved by a difference method of second order together with Newton or Newton-type iterative methods. Temimia, Ansari and Siddiqui [30] used an optimal homotopy asymptotic method to obtain a semianalytical solution. Peradze [23] proposed a new type of iteration method with Green function to solve the problem (1.1) directly and estimated the algorithm accuracy by the method of a priori inequalities. Recently, the so-called Bernoulli matrix method has been proposed to find approximate solutions of differential and integro-differential equations. The main characteristic of the Bernoulli matrix method is based on the operational matrices of differentiation and integration. The advantage of this technique is that by the fundamental matrix relations the Bernoulli series solution of differential and integral equations is reduced to a system of algebraic equations with unknown Bernoulli coefficients. For example, Bhrawy, Tohidi and Soleymani [6] constructed a Bernoulli matrix approach for high-order linear Fredholm integrodifferential equations. Tohidi, Bhrawy and Erfani [31] provided a numerical method based on Bernoulli operational matrix of derivatives for solving generalized pantograph equations. Tohidi, Erfani, Gachpazan and Shateyi [32] proposed a Tau method together with the Bernoulli polynomial approximation to solve nonlinear Lane-Emden type equations. Tohidi and Kili¸cman [33] proposed a collocation method based on the Bernoulli operational matrix for solving boundary value problems which arise from the problems in calculus of variation. Toutounian, Tohidi and Shateyi [35] developed a collocation method together with the Bernoulli polynomial approximations for high-order linear complex differential equations in rectangular domains. Toutounian and Tohidi [37, 34] proposed Bernoulli matrix methods for solving second order linear partial differential equations and onedimensional matrix hyperbolic equations of the first order. Tohidi and Zak [36] 4

ACCEPTED MANUSCRIPT

AN US

CR IP T

developed a Bernoulli matrix approach for solving the second-order linear matrix partial differential equations. In this paper, we construct a numerical method based on the Bernoulli polynomial approximation for Problem (1.1). In Section 2, a brief overview of Bernoulli polynomials and their relevant properties are given. In Section 3, we reduce the problem to a nonlinear finite dimensional system using the Bernoulli polynomial approximation and then provide an iterative process to solve such nonlinear system together with a collocation method. The convergence of the iterative scheme and some error estimates of the numerical scheme are discussed in Section 4 and Section 5, respectively. In Section 6, we illustrate the performance of our proposed method by some numerical examples and compare it with a finite element method and a finite difference method. Finally, some concluding remarks are given in Section 7. 2. Bernoulli polynomials and operational matrix of differentiation Bernoulli polynomials [20, 26] are widely used in number theory and classical analysis and are usually defined by the generating function ∞

ED

M

X tn text B (x) = , n et − 1 n=0 n!

|t| < 2π,

and thus the Bernoulli numbers Bn (0) can be obtained by the expansion ∞

PT

X t tn B (0) = . n et − 1 n=0 n!

The Bernoulli polynomials for n ≥ 1 possess the properties [8]

AC

CE

n   X n Bn (x) = Bk (0)xn−k , k k=0

and

Z

0

1

Bn (x)dx = 0,

Bn (x + 1) − Bn (x) = nxn−1 ,

Bn0 (x) = nBn−1 (x).

5

ACCEPTED MANUSCRIPT

B 0 (x)T

CR IP T

Define B(x) = [B0 (x), B1 (x), · · · , BN (x)]. It follows from the properties that the first derivative of B(x) can be expressed in the matrix form      B0 (x) B00 (x) 0 0 0 ··· 0 0 0       B10 (x)   1 0 0 · · · 0 0 0   B1 (x)        B 0 (x)   0 2 0 · · ·   0 0 0  2   B2 (x)     (2.1)   , = . . . . .. .. .. .. ..       .. .. .. . . . . . . .        B 0 (x)   0 0 0 · · · N − 1 0 0   B   N −1   N −1 (x)    0 BN (x) BN (x) 0 0 0 ··· 0 N 0 {z } | {z }| {z } | M

AN US

B(x)T

where M ∈ R(N +1)×(N +1) is called the operational matrix of differentiation. Moreover, the k-th derivative of B(x) can be expressed as B (k) (x) = B(x)(M T )k ,

(2.2)

M

3. Bernoulli collocation method

k = 0, 1, 2, · · · .

ED

The main objective of this section is to design a numerical scheme using the truncated Bernoulli series to approximate the solution of the model problem and provide an iterative process to solve the resulted nonlinear system. We first transform the interval [0, L] into the reference domain [0, 1] by the coordinate transformation

PT

x=

t , L

∀ t ∈ [0, L].

(3.1)

AC

CE

Set yˆ(x) = y( Lt ) and fˆ(x) = f ( Lt ). So Eq. (1.1) becomes ( R1 0 yˆ(4) (x) − (εL2 + 2 0 [ˆ y (x)]2 dx)ˆ y 00 (x) = L4 fˆ(x), 0 < x < 1, yˆ(0) = yˆ(1) = yˆ00 (0) = yˆ00 (1) = 0.

We assume that the solution yˆ(x) of Eq. (3.1) and the right term fˆ(x) are approximated by the truncated Bernoulli series with the first N + 1 terms yˆ(x) ' yˆN (x) =

N X n=0

6

an Bn (x) = B(x)A,

(3.2)

ACCEPTED MANUSCRIPT

and

N X

fˆ(x) ' fˆN (x) =

fˆn Bn (x) = B(x)Fˆ ,

(3.3)

n=0

(k)

yˆN (x) = B (k) (x)A = B(x)(M T )k A,

CR IP T

respectively, where A = [a0 , a1 , · · · , aN ]T and Fˆ = [fˆ0 , fˆ1 , · · · , fˆN ]T . Using the property (2.2), the k-th derivative of yˆN (x) can be expressed as k = 0, 1, 2 · · · .

(3.4)

AN US

The approximate solution yˆN (x) is obtained by solving the following system ( R1 0 (4) 00 (x) = L4 fˆN (x), 0 < x < 1, yN yˆN (x) − [εL2 + 2 0 (ˆ yN (x))2 dx]ˆ (3.5) 00 00 yˆN (0) = yˆN (1) = yˆN (0) = yˆN (1) = 0.

M

Due to the integral part of (3.5), a system of nonlinear equations will be obtained and also a full size matrix will be introduced when Newton’s method is adopted to solve the nonlinear system. To avoid this type of difficulty, Semper [27] proposed the following iterative procedure. Denote Z 1 0 (yN (x))2 dx. (3.6) α=2 0

0 < x < 1,

(3.7)

PT

ED

Then Eq. (3.5) becomes ( (4) 00 (x) = L4 fˆN (x), yˆN (x) − (εL2 + α)ˆ yN 00 00 (1) = 0. (0) = yˆN yˆN (0) = yˆN (1) = yˆN

AC

CE

Take α0 = 0. For k = 1, 2, . . . , Eq. (3.7) can be solved by the following iterative process ( (4) 00 yˆN,k (x) − (εL2 + αk−1 )ˆ yN,k (x) = L4 fˆN (x), 0 < x < 1, (3.8) 00 00 yˆN,k (0) = yˆN,k (1) = yˆN,k (0) = yˆN,k (1) = 0, where

αk−1 = 2

Z

0

1

0 [ˆ yN,k−1 (x)]2 dx.

(3.9)

This iterative process will be stopped when |αk − αk−1 | < Tol, where “Tol” is a given tolerance. 7

ACCEPTED MANUSCRIPT

CR IP T

Let xi = Ni , i = 0, 1, · · · , N − 4. The collocation method for solving (3.8) is to seek a polynomial yeN,k (x) of degree N defined on [0, 1] such that ( (4) 00 yeN,k (xi ) − (εL2 + αk−1 )e yN,k (xi ) = L4 fˆN (xi ), i = 0, 1, · · · , N − 4, (3.10) 00 00 (1) = 0. (0) = yeN,k yeN,k (0) = yeN,k (1) = yeN,k For this purpose, let yeN,k (x) =

N X

ak,n Bn (x) = B(x)Ak ,

n=0

Ak = [ak,0 , ak,1 , · · · , ak,N ]T .

AN US

Using the property (2.2), the first equation of (3.10) becomes   B(xi ) (M T )4 − (εL2 + αk−1 )(M T )2 Ak = L4 fˆN (xi ), i = 0, 1, . . . , N − 4, (3.11)

B1 (x0 ) B1 (x1 ) .. .

··· ··· .. .

BN (x0 ) BN (x1 ) .. .

ED

  B=  

B0 (x0 ) B0 (x1 ) .. .

B0 (xN −4 ) B1 (xN −4 ) · · · BN (xN −4 )

PT

where 

M

which can be written in a compact form   B (M T )4 − (εL2 + αk−1 )(M T )2 Ak = F,

For simplicity, we denote



  ,  



  F=  

L4 B(x0 )fˆ0 L4 B(x1 )fˆ1 .. . L4 B(xN −4 )fˆN −4

(3.12) 

  .  

CE

Wk−1 = B[(M T )4 − (εL2 + αk−1 )(M T )2 ].

So (3.12) becomes

(3.13)

AC

Wk−1 Ak = F.

Meanwhile, the corresponding boundary conditions are approximated by B(0)Ak = 0, B(1)Ak = 0, B(0)(M T )2 Ak = 0, B(1)(M T )2 Ak = 0, (3.14)

which can also be expressed as U Ak = 0, 8

(3.15)

ACCEPTED MANUSCRIPT

where U = [B(0), B(1), B(0)(M T )2 , B(1)(M T )2 ]T .

Uj = [uj0 , uj1 , uj2 , · · · , ujN ],

CR IP T

Denote the (i, j)-th element of Wk−1 by wi,j (i = 0, 1, . . . , N − 4, j = 0, 1, . . . , N ) and the j-th row of U by j = 0, 1, 2, 3.

Consequently, the solution of Eq. (3.10) is determined by solving the following system fk−1 Ak = Fe, W (3.16) where

0

wN −4 u01 u11 u21 u31

··· ··· .. .

w02 w12 .. . 1

wN −4 u02 u12 u22 u32

2



w0N w1N .. .

· · · wN −4,N ··· u0N ··· u1N ··· u2N ··· u3N

L4 B(x0 )fˆ0 L4 B(x1 )fˆ1 .. .

             4  ˆ  , Fe =  L B(xN −4 )fN −4   0       0     0   0

fk−1 ] =rank[W fk−1 , Fe] = N + 1, then we have If rank[W

PT



AN US

       wN −4 =   u00   u10   u 20  u30

w01 w11 .. .

M

fk−1 W

w00 w10 .. .

ED





       .       

fk−1 )−1 F˜ . Ak = (W

AC

CE

After the coefficient vector Ak = [ak,0 , ak,1 , · · · , ak,N ] is determined, we then obtain the approximate k-th iterative solution yˆN,k (x) '

N X

ak,n Bn (x) = B(x)Ak .

n=0

Remark 1. The proposed Bernoulli collocation method can be used not only to solve the model (1.1) with the Sturm-Liouville type boundary conditions, but also can be extended to deal with other types of boundary conditions by replacing (3.14) correspondingly. We give some examples as follows.

9

ACCEPTED MANUSCRIPT

• For the clamped boundary conditions y(0) = y(1) = y 0 (0) = y 0 (1) = 0, (3.14) is replaced with

CR IP T

B(0)Ak = 0, B(1)Ak = 0, B(0)M T Ak = 0, B(1)M T Ak = 0. • For the free boundary conditions y 00 (0) = y 00 (1) = y 000 (0) = y 000 (1) = 0, (3.14) is replaced with

B(0)(M T )2 Ak = 0, B(1)(M T )2 Ak = 0, B(0)(M T )3 Ak = 0, B(1)(M T )3 Ak = 0.

AN US

• For the clamped-free boundary conditions such as y(0) = y(1) = y 000 (0) = y 000 (1) = 0, (3.14) is replaced with B(0)Ak = 0, B(1)Ak = 0, B(0)(M T )3 Ak = 0, B(1)(M T )3 Ak = 0. 4. Convergence of iteration

Assume that y(x) ∈ C 2 [0, 1] and denote

M

R1 ky(x)k∞ = maxx∈[0,1] |y(x)|, ky(x)k2 = 0 y 2 (x)dx, R1 R1 |y(x)|21 = 0 (y 0 (x))2 dx, |y(x)|22 = 0 (y 00 (x))2 dx.

ED

We recall the following inequalities which will be used in the convergence analysis. Lemma 4.1. [29] If y(x) ∈ C 2 [0, 1] and y(0) = y(1) = 0, then

AC

CE

PT

1 ky(x)k∞ ≤ |y(x)|1 , 2 1 ky(x)k ≤ √ |y(x)|1 , 6 1 |y(x)|1 ≤ √ |y(x)|2 . 6

(4.1) (4.2) (4.3)

Lemma 4.2. Assume that yˆ(x) is the solution of (3.1). Then we have the following inequality " # 12 1 L4 |ˆ y (x)|1 ≤ p kfˆ(x)k 2 . (4.4) 4 3(εL2 + 6) 10

ACCEPTED MANUSCRIPT

0

It follows from (4.2), (4.3) and Schwartz inequality that

L4 (εL2 + 6 + 2|ˆ y (x)|21 )|ˆ y (x)|1 ≤ √ kfˆ(x)k. 6

CR IP T

Proof. Multiplying both sides of Eq. (3.1) by yˆ(x) and then integrating the obtained equality with respect to x from 0 to 1, we obtain Z 1 4 2 2 2 2 fˆ(x)ˆ y (x)dx. y (x)|1 = L |ˆ y (x)|2 + (εL + 2|ˆ y (x)|1 )|ˆ

AN US

Applying the inequality 2ab ≤ a2 + b2 , ∀a, b ∈ R to the left side, we get p L4 2 2(εL2 + 6) |ˆ y (x)|21 ≤ √ kfˆ(x)k, 6 which yields the desired result. This completes the proof. From Lemma 4.2, we have the following corollary.

M

Corollary 4.1. Assume that yˆN (x) is the solution of Eq. (3.5). Then the following result holds # 21 " 1 L4 kfˆN (x)k 2 . (4.5) |ˆ yN (x)|1 ≤ p 2 4 3(εL + 6)

ED

Lemma 4.3. Assume that yˆN,k (x) is the solution of (3.8). Then we have |ˆ yN,k (x)|1 ≤ √

L4 kfˆN (x)k. 6(εL2 + 6)

(4.6)

CE

PT

Proof. Multiplying both sides of Eq. (3.8) by yˆN,k (x) and integrating it with respect to x from 0 to 1, we get Z 1  2  2 2 4 2 fˆN (x)ˆ yN,k (x)dx. |ˆ yN,k (x)|2 + εL + 2|ˆ yN,k−1 (x)|1 |ˆ yN,k (x)|1 = L 0

AC

Using (4.2), (4.3) and Schwartz inequality, we obtain

Thus,

 2  L4 εL + 6 + 2|ˆ yN,k−1 (x)|21 |ˆ yN,k (x)|1 ≤ √ kfˆN (x)k. 6

L4 |ˆ yN,k |1 ≤ √ kfˆN (x)k. 6(εL2 + 6) This completes the proof. 11

ACCEPTED MANUSCRIPT

√

4

L4 3(εL2 +6)

 12

1 2

kfˆ(x)k , c2 =

Theorem 4.1. Assume that 0 < q1 = generated by (3.9) converges to α.

√

L4 6(εL2 +6)

2(ˆ c1 +c2 )c2 6−εL2 −2ˆ c21

kfˆN (x)k, cˆ1 =



√

4

L4 3(εL2 +6)

 12

< 1. Then the sequence {αk }∞ k=0

Proof. Subtracting (3.9) from (3.6) yields Z 1 0 0 0 0 (ˆ yN (x) − yˆN,k (x))(ˆ yN (x) + yˆN,k (x)) dx. |α − αk | ≤ 2 0

AN US

It follows from Schwartz inequality and Corollary 4.1 that

|α − αk | ≤ 2|ˆ yN (x) − yˆN,k (x)|1 (|ˆ yN (x)|1 + |ˆ yN,k |1 ) ≤ 2(ˆ c1 + c2 )|ˆ yN (x) − yˆN,k (x)|1 . Subtracting (3.8) from (3.7), we have

(4.7)

00 00 00 [ˆ yN (x) − yˆN,k (x)](4) = (εL2 + α)(ˆ yN (x) − yˆN,k (x)) + (α − αk−1 )ˆ yN,k (x).

M

Multiplying both sides of the above equation by yˆN (x) − yˆN,k (x) and integrating it with respect to x from 0 to 1, we can obtain

ED

|ˆ yN (x)−ˆ yN,k (x)|22 ≤ (εL2 +2ˆ c21 )|ˆ yN (x)−ˆ yN,k (x)|21 +(α−αk−1 )|ˆ yN,k (x)|1 |ˆ yN (x)−ˆ yN,k (x)|1 . Using (4.3), we have

PT

|ˆ yN (x) − yˆN,k (x)|1 ≤

1 (α − αk−1 )|ˆ yN,k (x)|1 . 6 − εL2 − 2ˆ c21

(4.8)

CE

It follows from Lemmas 4.2 and 4.3 and (4.7) that |α − αk | ≤ q1 |α − αk−1 |, which implies that the iterative sequence is convergent. This completes the proof. 5. Error estimation

AC

We need some classical results on the Bernoulli polynomial approximation. By using the result on expansion of an arbitrary function in Bernoulli polynomials [15], if g(x) has a continuous derivative of order N on [0, 1] then for any x ∈ [0, 1] we have the relation g(x) = gN (x) + RN [g](x), (5.1) where gN (x) and RN [g](x) are given in the following lemma. 12

1

kfˆN (x)k 2 .

CR IP T

Denote c1 =



ACCEPTED MANUSCRIPT

(5.2)

CR IP T

Lemma 5.1. [1, 15] Suppose that g(x) ∈ C N [0, 1]. Then we have R1 P Bj (x) (j−1) (1) − g (j−1) (0)], gN (x) = 0 g(x)dx + N j=1 j! [g R 1 ∗ RN [g](x) = − N1 ! 0 BN (x − t)g (N ) (t)dt,

∗ where BN (x) = BN (x − [x]) and [x] denotes the largest integer not greater than x.

AN US

Lemma 5.2. [1] Assume that g ∈ C N [0, 1]. Then the coefficients gn in (5.2) satisfy Z 1 1 (n) g (x)dx, n = 0, 1, 2, · · · , N, gn = n! 0 which implies 1 max |g (n) (x)|. |gn | ≤ n! x∈[0,1]

Remark 2. Some new theorems on Bernoulli polynomial approximation are given in [33, 35].

M

For brevity, we denote q2 = 6 + εL2 − 2(c1 + cˆ1 )ˆ c1 , and q3 =

1 . 6−εL2 −2ˆ c21

ED

Theorem 5.1. Assume that yˆ(x), yˆN (x) and yˆN k (x) are the solutions of Eqs.(3.1), (3.5) and (3.8), respectively, and q2 > 0. Then we have L4 −1 √ q kRN [fˆ](x)k∞ + c2 q3k−1 |α − α0 |, 6 2 4 yˆN,k (x)k∞ ≤ 2L√6 q2−1 kRN [fˆ](x)k∞ + c2 q3k−1 |α − α0 |.

|ˆ y (x) − yˆN,k (x)|1 ≤

(5.4)

PT

kˆ y (x) −

(5.3)

CE

Proof. We first find an upper bound for the error between the exact solution yˆ(x) and the approximate solution yˆN (x). Subtracting Eq. (3.5) from Eq. (3.1), we obtain

AC

[ˆ y (x) − yˆN (x)](4) − εLh2 [ˆ y (x) − yˆN (x)]00 i R1 0 R1 0 00 − 2 0 (ˆ y (x))2 dxˆ y 00 (x) − 0 (ˆ yN (x))2 dxˆ yN (x) = L4 RN [fˆ](x),

and hence

(4)

[ˆ y (x) − yˆN (x)]

00

Z

1

00 − εL [ˆ y (x) − yˆN (x)] − 2 (ˆ y 0 (x))2 dx[ˆ y 00 (x) − yˆN (x)] 0  Z 1 0 2 0 2 00 + [(ˆ y (x)) − (ˆ yN (x)) ]dxˆ yN (x) = L4 RN [fˆ](x). 2

0

13

ACCEPTED MANUSCRIPT

0

It follows from (4.2) and (4.3) and Schwartz inequality that

CR IP T

Multiplying both sides of the above equation by yˆ(x) − yˆN (x) and integrating it with respect to x from 0 to 1, we derive   Z 1 2 2 0 2 |ˆ y (x) − yˆN (x)|2 + εL + 2 [ˆ y (x)] dx |ˆ y (x) − yˆN (x)|21 0 Z 1 Z 1 0 0 0 0 0 0 − 2 [ˆ y (x) + yˆN (x)][ˆ yN (x) − yˆ (x)]dx yˆN (x)[ˆ y 0 (x) − yˆN (x)]dx 0 0 Z 1 [ˆ y (x) − yˆN (x)]dx. = L4 RN [fˆ](x)

AN US

  L4 6 + εL2 + 2|ˆ y (x)|21 − 2[|ˆ y (x)|1 + |ˆ yN (x)|1 ]|ˆ yN (x)|1 |ˆ y (x)−ˆ yN (x)|1 ≤ √ kRN [fˆ](x)k∞ . 6

Application of Lemma 4.1 and Corollary 4.1 yields 

 L4 6 + εL2 − 2(c1 + cˆ1 )ˆ c1 |ˆ y (x) − yˆN (x)|1 ≤ √ kRN [fˆ](x)k∞ . 6

M

Combining (4.8) with (5.5), we get

(5.5)

ED

|ˆ y (x) − yN,k (x)|1 ≤ |ˆ y (x) − yˆN (x)|1 + |ˆ yN (x) − yˆN,k (x)|1 4 L ≤ √ q2−1 kRN [fˆ](x)k∞ + c2 q3k−1 |α − α0 |. 6

PT

The inequality (5.4) then follows from (4.1). This completes the proof. 6. Numerical experiments

CE

We illustrate the performance of our proposed Bernoulli collocation method (denoted by BCM) by some numerical examples and compare it with the finite element method (denoted by FEM [28]) and the finite difference method (denoted by FDM [9]). The iterative process will be stopped when |αk∗ − αk∗ −1 | < Tol.

AC

Example 6.1. Consider the fourth order nonlinear integro-differential equation [9, 28] ( Rπ y (4) (t) − 2y 00 (t) − π2 0 [y 0 (t)]2 dty 00 (t) = −4 sin t, 0 < t < π, (6.1) y(0) = y(π) = y 00 (0) = y 00 (π) = 0.

The exact solution is y(t) = − sin t. 14

ACCEPTED MANUSCRIPT

Table 1: Numerical solutions of Problem (6.1).

t

Exact solution

BCM N =15

FEM

π 10 2π 10 3π 10 4π 10 5π 10 6π 10 7π 10 8π 10 9π 10

-0.30901699 -0.58778525 -0.80901699 -0.95105652 -1.00000000 -0.95105652 -0.80901699 -0.58778525 -0.30901699

-0.30901699 -0.58778525 -0.80901699 -0.95105651 -0.99999999 -0.95105651 -0.80901699 -0.58778525 -0.30901699

FDM

N =40

N =20

AN US

N =20

CR IP T

We solve this problem by BCM, FEM and FDM and define their error function as eN,k∗ (t) = |y(t) − y N,k∗ (t)|, where y(t) is the exact solution and y N,k∗ (t) is the computed approximate solution. In Table 1, we list the exact solution, our approximate solution with N = 15 and Tol= 10−12 , the approximate solutions by FEM and FDM N = 20, 40. In Table 2, we record the numbers of iterations, the maximum errors and the values of αk∗ for different N and Tol of our proposed method. We also plot the errors of BCM with N = 9, 11 and FDM with

-0.30901698 -0.58778522 -0.80901695 -0.95105646 -0.99999994 -0.95105646 -0.80901695 -0.58778522 -0.30901698

-0.30965218 -0.58899344 -0.81067993 -0.95301141 -1.00205550 -0.95301141 -0.81067993 -0.58899344 -0.30965218

-0.30917580 -0.58808733 -0.80943276 -0.95154528 -1.00051392 -0.95154528 -0.80943276 -0.58808733 -0.30917580

ED

M

-0.30901669 -0.58778467 -0.80901620 -0.95105558 -0.99999901 -0.95105558 -0.80901620 -0.58778467 -0.30901669

N =40

Tol

6 8 10 12 14

10−6 10−6 10−8 10−8 10−10

k∗

keN,k∗ (t)k∞

αk ∗

α

26 22 29 29 36

4.5827e-002 2.5784e-003 7.8762e-005 1.4536e-006 1.8089e-008

1.087203524 0.995483813 1.000134873 0.999997538 1.000000031

0.903868143 1.009355650 0.999500253 1.000018016 0.999999534

AC

CE

N

PT

Table 2: Numerical results of our proposed method for Problem (6.1).

N = 40, 80 in Figure 1 (left), and the errors of BCM with N = 14, 16 and FEM with N = 40, 80 in Figure 1 (right), respectively. We observe that BCM has higher accuracy than FEM and FDM. 15

ACCEPTED MANUSCRIPT

−8

7

−3

x 10

1.2 BCM for N=14 BCM for N=16 FEM for N=40 FEM for N=80

6

x 10

BCM for N=9 BCM for N=11 FDM for N=40 FDM for N=80

1

5 0.8

0.6

3 0.4 2 0.2

1

0

0

0.5

1

1.5

2

2.5

0

3

0

t

0.5

1

CR IP T

Error

Error

4

1.5

2

2.5

3

t

AN US

Figure 1: Error functions of BCM and FEM (left), BCM and FDM (right) for (6.1).

Example 6.2. Consider the following static beam problem [30] ( R1 y (4) (t) − y 00 (t) − 0 [y 0 (t)]2 dty 00 (t) = 1, 0 < t < 1, y(0) = y(1) = y 00 (0) = y 00 (1) = 0.

(6.2)

ED

M

The exact solution of this problem is unknown. We solve (6.2) by BCM, FEM and FDM and use the following residual function adopted in [22] to reflect the accuracy Z L (4) 2 00 0 2 00 [y N,k∗ (t)] dxy N,k∗ (t) − f (t) . RN,k∗ (t) = y N,k∗ (t) − εy N,k∗ (t) − L 0

CE

PT

The numbers of iterations of BCM, FEM and FDM with different Tol are recorded in Figure 2 (left) and the residual functions of these methods with Tol= 10−10 are plotted in Figure 2 (right). From the numerical results, we may conclude that BCM is more efficient and accurate than FEM and FDM for Problem (6.2).

AC

Example 6.3. Consider the following nonlinear integro-differential equation ( R1 y (4) (t) − y 00 (t) − 2 0 [y 0 (t)]2 dty 00 (t) = −t, 0 < t < 1, (6.3) y(0) = y(1) = y 00 (0) = y 00 (1) = 0. The exact solution of this problem is also unknown. We solve Problem (6.3) by BCM, FEM and FDM. The numbers of iterations of BCM, FEM and FDM with different Tol are recorded in Figure 3 (left) and their residual functions with 16

ACCEPTED MANUSCRIPT

35

30

BCM FEM FDM

BCM for (N,k*)=(15,3)

−2

10

FEM for (N,k*)=(40,34) FDM for (N,k*)=(40,34)

−4

10

−6

15

−8

10

−10

10

10

5

10

0

−12

−14

Tol=1.0E−6

Tol=1.0E−8

Tol=1.0E−10

10

0

0.1

0.2

0.3

0.4

0.5 t

CR IP T

10

20

RN,k*(t)

Numbers of iterations

25

0.6

0.7

0.8

0.9

1

AN US

Figure 2: Numbers of iterations (left) and residual functions (right) of BCM, FEM and FDM for Problem (6.2).

M

Tol= 10−10 are plotted in Figure 3 (right). From the numerical results, we may also conclude that BCM is more efficient and accurate than FEM and FDM for Problem (6.3).

PT

ED

Example 6.4. Consider the model of nonlinear beam on elastic bearings [16]  (4) R1 y (t) − y 00 (t) − 2 0 [y 0 (t)]2 dty 00 (t) = − 107 t2 + 214 t + 6, 0 < t < 1,  35 105    y 0 (0) = y 0 (1) = 0, (6.4) 000  y (0) = 24y(0) + 8,    000 y (1) = −24y(1) − 8,

AC

CE

The exact solution is y(t) = 41 t4 − 13 t3 + 41 . This type of model is used to describe the bending equilibrium state of a beam of length 1 subjected to a force function and the elastic forces at the ends given by bearings (see [16]). Following the procedure in Section 3, we can obtain its approximate solution by BCM for different N and Tol. The computed results are reported in Table 3 and error functions for different N are plotted in Figure 4. We observe that our proposed method is convergent and perform very accurate for Problem (6.4). 7. Concluding remarks In this paper, we propose a numerical scheme for solving a class of nonlinear fourth order integro-differential equations of Kirchhoff type which arises in the 17

ACCEPTED MANUSCRIPT

35

−2

BCM FEM FDM

FEM for (N,k*)=(40,33)

−6

10

20

RN,k*(t)

Numbers of itertions

FDM for (N,k*)=(40,33)

−4

10

25

15

−8

10

−10

10

10

10

5

10

0

BCM for (N,k*)=(15,3)

10

−12

−14

−16

Tol=1.0E−6

Tol=1.0E−8

Tol=1.0E−10

10

0

0.1

0.2

0.3

0.4

0.5 t

CR IP T

30

0.6

0.7

0.8

0.9

1

AN US

Figure 3: Numbers of iterations (left) and residual functions (right) of BCM, FEM and FDM for Problem (6.3).

Table 3: Numerical results for Problem (6.4).

Tol

k∗

keN,k∗ (t)k∞

7 9 11 13 15

10−7 10−7 10−9 10−9 10−11

5 5 6 6 8

6.9397e-010 9.2742e-010 3.0507e-011 3.0494e-011 1.2780e-013

αk∗

α

0.01904761961 0.01904761979 0.01904761902 0.01904761902 0.01904761905

0.01904761905 0.01904761905 0.01904761905 0.01904761905 0.01904761905

PT

ED

M

N

AC

CE

study of transverse vibrations of a hinged beam. The main idea is to reduce the problem to a nonlinear algebraic system based on the Bernoulli polynomials and their operational matrix of differentiation and then solve the nonlinear system by an iterative process together with a collocation method. Several numerical examples have been implemented to illustrate the performance of our proposed method. The comparison with FEM and FDM shows that our proposed method is more accurate, efficient and easy to use. Moreover, our proposed method can be extended to solve boundary value problems for integro-differential equation with nonlinear boundary conditions and to the problems with nonpositive (or nonnegative) f (y, y 0 , y 00 , y 000 ) under the same boundary conditions.

18

ACCEPTED MANUSCRIPT

−9

10

−10

10

−11

−12

10

−13

10

e7,5(t)

−14

10

e11,6(t) e15,8(t)

−15

10

0

0.1

0.2

0.3

0.4

0.5 t

0.6

0.7

0.8

0.9

1

CR IP T

eN,k*(t)

10

AN US

Figure 4: Error functions of BCM with N = 7, 11, 15 for Problem (6.4).

Acknowledgments

M

The authors are grateful to the referees for their helpful suggestions and valuable comments, which have greatly improved the presentation of the paper. References

ED

[1] M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, National Bureau of Standards, Wiley, New York, 1972.

PT

[2] J. Anderson, Computational Fluid Dynamics, Volume 206, Springer, 1995.

CE

[3] A. Arosio, A geometrical nonlinear correction to the Timoshenko beam equation, Nonlinear Anal. TMA, 47(2001) 729–740. [4] J. M. Ball, Initial-boundary value problems for an extensible beam, J. Math. Anal. Appl., 42(1973) 61–90.

AC

[5] J. M. Ball, Stability theory for an extensible beam, J. Differential Equations, 14(1973) 399–418.

[6] A.H. Bhrawy, E. Tohidi and F. Soleymani, A new Bernoulli matrix method for solving high-order linear and nonlinear Fredholm integro-differential equations with piecewise intervals, Appl. Math. Comput., 219(2012) 482– 497. 19

ACCEPTED MANUSCRIPT

[7] P. Biler, Remark on the decay for damped string and beam equations, Nonlinear Anal. TMA, 10(1986) 839–842.

CR IP T

[8] F. A. Costabile and F. Dell’accio, Expansion over a rectangle of real functions in Bernoulli polynomials and applications, BIT Numer. Math., 41(2001) 451– 464. [9] Q. A. Dang and V. T. Luan, Iterative method for solving a nonlinear fourth order boundary value problem, Comput. Math. Appl., 60(2010) 112–121.

AN US

[10] C. E. Elmer and E. S. Van Vleck, Traveling wave solutions for bistable differential-difference equations with periodic diffusion. SIAM J. Appl. Math., 61(2001) 1648–1679.

[11] E. Feireisl, Exponential attractors for nonautonomous systems: Long-time behaviour of vibrating beams, Math. Meth. Appl. Sci., 15(1992) 287-297.

M

[12] B. Guo and W. Guo, Adaptive stabilization for a Kirchhoff-type nonlinear beam under boundary output feedback control, Nonlinear Anal. TMA, 66(2007) 427– 441.

ED

[13] G. Kirchhoff, Vorlesungen u ¨ber Mathematische Physik, Volume 1, Teubner, 1883.

PT

[14] S. Kou´emou-Patcheu, On a global solution and asymptotic behaviour for the generalized damped extensible beam equation, J. Differential Equations, 135(1997) 299–314. [15] V. I. Krylov, Approximate calculation of integrals, Dover publications, Mineola, New York, 1962.

CE

[16] T. F. Ma, Existence results for a model of nonlinear beam on elastic bearings, Appl. Math. Lett., 13(2000) 11–15.

AC

[17] T. F. Ma, Existence results and numerical solutions for a beam equation with nonlinear boundary conditions, Appl. Numer. Math., 47(2003) 189–196.

[18] T. F. Ma, Positive solutions for a nonlocal fourth order equation of Kirchhoff type, Discrete Contin. Dyn. Syst., Supplement (2007) 694–703.

20

ACCEPTED MANUSCRIPT

[19] T. F. Ma and A. L. M. Martinezb, Positive solutions for a fourth order equation with nonlinear boundary conditions, Math. Comput. Simulation, 80(2010) 2177–2184.

CR IP T

[20] P. Natalini and A. Bernardini, A generalization of the Bernoulli polynomials, J. Appl. Math., 3(2003) 155–163.

[21] M. R. Ohm, H. Y. Lee and J. Y. Shin, Error estimates of finite-element approximations for a fourth-order differential equation, Comput. Math. Appl., 52(2006) 283–288.

AN US

[22] F.A. Oliveira, Collocation and residual correction, Numer. Math, 36(1980) 27–31.

[23] J. Peradze, An iteration method for the Kirchhoff static beam, Bull. TICMI, 16(2012) 27–33. [24] J. Peradze, A numerical algorithm for a Kirchhoff-type nonlinear static beam, J. Appl. Math. , 2009, Art. ID 818269, 12 pp.

M

[25] D. C. Pereira, Existence, uniqueness and asymptotic behavior for solutions of the nonlinear beam equation, Nonlinear Anal. TMA, 14(1990) 613–623.

ED

[26] G. Rz¸adkowski and S. Lepkowski, A generalization of the Euler-Maclaurin summation formula: an application to numerical computation of the fermidirac integrals, J. Sci. Comput., 35(2008) 63–74.

PT

[27] B. Semper, Finite element methods for suspension bridge models, Comput. Math. Appl., 26(1993) 77–91.

CE

[28] J. Y. Shin, Finite-element approximation of a fourth-order differential equation, Comput. Math. Appl., 35(1998) 95–100.

AC

[29] Z. Sun, Numerical Solution of Partial Differential Equations, Science Press, Beijing, 2005. [30] H. Temimi, A. R. Ansari and A. M. Siddiqui, An approximate solution for the static beam problem and nonlinear integro-differential equations, Comput. Math. Appl., 62(2011) 3132–3139.

21

ACCEPTED MANUSCRIPT

[31] E. Tohidi, A.H. Bhrawy and K. Erfani, A collocation method based on Bernoulli operational matrix for numerical solution of generalized pantograph equation, Appl. Math. Modelling, 37(2013) 4283–4294.

CR IP T

[32] E. Tohidi, Kh. Erfani, M. Gachpazan and S. Shateyi, A new Tau method for solving nonlinear Lane-Emden type equations via Bernoulli operational matrix of differentiation, J. Appl. Math., 2013, Article ID 850170, 9 pp.

[33] E. Tohidi and A. Kili¸cman, A collocation method based on the Bernoulli operational matrix for solving nonlinear BVPs which arise from the problems in calculus of variation, Math. Probl. Eng., 2013, Art. ID 757206, 9 pp.

AN US

[34] E. Tohidi and F. Toutounian, Convergence analysis of Bernoulli matrix approach for one-dimensional matrix hyperbolic equations of the first order, Comput. Math. Appl., 68(2014) 1–12.

M

[35] E. Tohidi, F. Toutounian and S. Shateyi, A collocation method based on the Bernoulli operational matrix for solving high-order linear complex differential equations in a rectangular domain, Abstr. Appl. Anal., 2013, Art. ID 823098, 12 pp.

ED

[36] E. Tohidi and M. K. Zak, A new matrix approach for solving second-order linear matrix partial differential equations, Mediterr. J. Math., 2015, DOI 10.1007/s00009-015-0542-2.

PT

[37] F. Toutounian and E. Tohidi, A new Bernoulli matrix method for solving second order linear partial differential equations with the convergence analysis, Appl. Math. Comput., 223(2013) 298–310.

AC

CE

[38] S. Woinowsky-Krieger, The effect of an axial force on the vibration of hinged bars, J. Appl. Mech., 17(1950) 35–36.

22