^{1}

^{*}

^{2}

^{2}

Euler-Bernoulli beam equation is very important that can be applied in the field of mechanics, science and technology. Some authors have put forward many different numerical methods, but the precision is not enough high. In this paper, we will illustrate the high-precision numerical method to solve Euler-Bernoulli beam equation. Three numerical examples are studied to demonstrate the accuracy of the present method. Results obtained by our method indicate new algorithm has the following advantages: small computational work, fast convergence speed and high precision.

Under the free vibration, the equation of transverse motion of a uniform Euler-Bernoulli beam is determined by a partial differential equation [

E I ∂ 4 u ( x , t ) ∂ x 4 + m ∂ 2 u ( x , t ) ∂ t 2 + r e ∂ u ( x , t ) ∂ t + r i ∂ ∂ t ( ∂ 4 u ( x , t ) ∂ x 4 ) = 0 (1)

where E I , m , r e , and r i are the beams flexural rigidity and mass per unit length, external (air) damping coefficient, and Kelvin-Voigt internal damping coefficient [

In [

E I ( x ) ∂ 4 u ( x , t ) ∂ x 4 + m ( x ) ∂ 2 u ( x , t ) ∂ t 2 + r e ( x ) ∂ u ( x , t ) ∂ t + r i ( x ) u ( x , t ) = f ( x , t ) (2)

In this paper, we consider the following variable coefficient Euler-Bernoulli beam equation

∑ i = 0 4 d i ( x ) ∂ i u ( x , t ) ∂ x i + d 5 ( x ) ∂ 2 u ( x , t ) ∂ t 2 + d 6 ( x ) ∂ u ( x , t ) ∂ t + d 7 ( x ) ∂ ∂ t ( ∂ 4 u ( x , t ) ∂ x 4 ) = f ( x , t , u ) , (3)

the initial conditions are expressed as:

u ( x , 0 ) = ϕ ( x ) , u t ( x , 0 ) = ψ ( x ) , x ∈ [ 0 , l ] , (4)

and the boundary conditions are expressed as:

u ( 0 , t ) = p 0 ( t ) , u ( l , t ) = p 1 ( t ) , u x ( 0 , t ) = p 2 ( t ) , u x ( l , t ) = p 3 ( t ) , t ≥ 0 , (5)

where, u ( x , t ) represents an unknown function at position x and time t. d i ( x ) , i = 0 , 1 , ⋯ , 7 is a known function. If f ( x , t , u ) is nonlinear of u, the equation is a nonlinear Partial differential equation. If f ( u ) is linear of u, then the Equation (3) is linear Partial differential equation.

The model (3) can be written as:

∑ i = 0 4 d i ( x ) ∂ i u ( x , t ) ∂ x i + d 5 ( x ) ∂ 2 u ( x , t ) ∂ t 2 + d 6 ( x ) ∂ u ( x , t ) ∂ t + d 7 ( x ) ∂ ∂ t ( ∂ 4 u ( x , t ) ∂ x 4 ) = f ( x , t ) , (6)

We introduce the barycentric Lagrange interpolation collocation method [

We consider a regular region W = [ 0, l ] × [ 0, T ] , the interval [ 0, l ] is divided into M different nodes and the interval [ 0, T ] is also divided into N different nodes. The nodes on the interval [ 0, l ] and [ 0, T ] constitute two groups of column vectors respectively. They are defined as:

x 0 = [ x 1 , x 2 , ⋯ , x M ] T , t 0 = [ t 1 , t 2 , ⋯ , t N ] T , (7)

then the matrixs X, T above the region Ω can be generated by the above column vectors x 0 , t 0 , and the form of matrixs X and T are as follows:

X T = [ ( x 0 ) T , ( x 0 ) T , ⋯ , ( x 0 ) T ] , T = [ t 0 , t 0 , ⋯ , t 0 ] . (8)

The matrix X, T are respectively written as column vectors x, t of N × M dimension:

x = [ X 1 , X 2 , ⋯ , X M × N ] T , t = [ T 1 , T 2 , ⋯ , T M × N ] T . (9)

The relations between the components x i , y j of the Formula (11) and the components X K , Y K of Formula (13) are as follows:

X k = X ( i − 1 ) N + j = x i , T k = T ( i − 1 ) N + j = t j , i = 1 , 2 , ⋯ , M ; j = 1 , 2 , ⋯ , N ; k = 1 , 2 , ⋯ , M × N . (10)

The barycentric interpolation of u ( x , t ) at nodes ( x i , t j ) can be written as [

u ( x , t ) = ∑ i = 1 M ∑ j = 1 N ξ i ( x ) η j ( t ) v i j , i = 1 , 2 , ⋯ , M ; j = 1 , 2 , ⋯ , N , (11)

where, ξ i ( x ) = ∏ k = 1 , k ≠ i M ( x − x k ) ∏ k = 1 , k ≠ i M ( x i − x k ) , i = 1 , 2 , ⋯ , M , η j ( t ) = ∏ k = 1 , k ≠ j N ( t − t k ) ∏ k = 1 , k ≠ j N ( t j − t k ) , j = 1 , 2 , ⋯ , N .

Use Formula (11), the l + k order partial derivative of function v ( x , t ) at nodes can be expressed as:

∂ l + k u ∂ x l ∂ t k = ∑ i = 1 M ∑ j = 1 N ξ i ( l ) ( x ) η j ( k ) ( t ) u i j , l , k = 1 , 2 , ⋯ . (12)

At nodes ( x p , t q ) , the function values of partial derivative are defined as:

u ( l , k ) ( x p , t q ) : = ∂ l + k u ( x p , t q ) ∂ x l ∂ t k = ∑ i = 1 M ∑ j = 1 N ξ i ( l ) ( x p ) η j ( k ) ( t q ) u p q , p = 1 , 2 , ⋯ , M ; q = 1 , 2 , ⋯ , N . (13)

The function values of the Formula (11) and the Formula (12) at the node form column vectors u , u ( l , k ) , and they are as follows:

u = [ u 1 , u 2 , ⋯ , u M × N ] T , u ( l , k ) = [ u 1 ( l , k ) , u 2 ( l , k ) , ⋯ , u M × N ( l , k ) ] T ,

u p = u ( X p , T p ) , u p ( l , k ) = u ( l , k ) ( X p , T p ) , p = 1 , 2 , ⋯ , M × N .

(13) can be written in following matrix form [

u ( l , k ) = D ( l , k ) u . (14)

In the above formula, D ( l , k ) = C ( l ) ⊗ D ( k ) is the Kronecker product of matrix C ( l ) and D ( k ) , which is also called ( l , k ) order partial differential matrix at nodes { ( x i , t j ) , i = 1 , 2 , ⋯ , M ; j = 1 , 2 , ⋯ , N } . C ( l ) is l order differential matrix on x direction nodes, and D ( k ) is k order differential matrix on t direction nodes.

Denote:

C ( 0 ) = I M , D ( 0 ) = I N , (15)

where, I M and I N are M order unit matrix and N order unit matrix respectively.

By using the notation (11), (12), (13), the calculation formula of barycentric Lagrange interpolation collocation method for model (3) can be written in following matrix form:

[ ∑ i = 0 4 d i a g ( d i ( x ) ) D ( i ,0 ) + ∑ i = 1 2 d i a g ( d 4 + i ( x ) ) D ( 0, i ) + d i a g ( d 7 ( x ) ) D ( 4,1 ) ] u = f . (16)

Here, diag is a symbol of diagonal matrix composed of vectors.

Given initial value v 0 , we can construct following linear iterative format:

∑ i = 0 4 d i ( x ) ∂ i u n ( x , t ) ∂ x i + d 5 ( x ) ∂ 2 u n ( x , t ) ∂ t 2 + d 6 ( x ) ∂ u n ( x , t ) ∂ t + d 7 ( x ) ∂ ∂ t ( ∂ 4 u n ( x , t ) ∂ x 4 ) = f ( x , t , u n − 1 ) , (17)

So the model (17) can be written as the following matrix:

[ ∑ i = 0 4 d i a g ( d i ( x ) ) D ( i ,0 ) + ∑ i = 1 2 d i a g ( d 4 + i ( x ) ) D ( 0, i ) + d i a g ( d 7 ( x ) ) D ( 4,1 ) ] u = f 0 . (18)

The discretization of initial boundary conditions requires the use of the barycenter interpolation Formula (11). By acting on the initial boundary conditions, the discrete formula of initial boundary conditions is given:

∑ i = 1 M ∑ j = 1 N ξ i ( x p ) η j ( 0 ) u i j = ϕ ( x p ) , ∑ i = 1 M ∑ j = 1 N ξ i ( x p ) η ′ j ( 0 ) u i j = φ ( x p ) , ∑ i = 1 M ∑ j = 1 N ξ i ( 0 ) η j ( t q ) u i j = p 0 ( t q ) , ∑ i = 1 M ∑ j = 1 N ξ i ( l ) η j ( t q ) u i j = p 1 ( t q ) , ∑ i = 1 M ∑ j = 1 N ξ ′ i ( 0 ) η j ( t q ) u i j = p 2 ( t q ) , ∑ i = 1 M ∑ j = 1 N ξ ′ i ( l ) η j ( t q ) u i j = p 3 ( t q ) . (19)

The boundary conditions are replaced by displacement method.

In this section, following the guidance of the discussions in Section 2, we will select appropriate free parameters and present some numerical simulations for preceding cases, which implies that our current method is a satisfactory and efficient algorithm.

Example 1 We consider the vibration Euler-Bernoulli beam equation of fixed-supported at both ends.

{ ∂ 2 ∂ x 2 ( E I ( x ) ∂ 2 u ( x , t ) ∂ x 2 ) + ∂ 2 u ( x , t ) ∂ t 2 = f ( x , t ) , u ( x , 0 ) = u t ( x , 0 ) = 0 , x ∈ [ 0 , 1 ] , u ( 0 , t ) = u ( 1 , t ) = u x ( 0 , t ) = u x ( 1 , t ) = 0 , t ≥ 0 , (20)

1) where E I ( x ) = 1 + sin π x , and the source term f ( x , t ) is determined by (21) consistent with the chosen solution. The exact solution of beam deflection is as follows:

u ( x , t ) = x ( 1 − x ) sin ( 4 π x ) t 2 e − t . (21)

By the proposed algorithm, we obtain the numerical solution and the absolute error which are given in Figures 1-3 and

2) where E I ( x ) = 1 + 1 4 exp ( − 40 ( x − 1 3 2 ) ) , and the source term f ( x , t ) is

determined by (22) consistent with the chosen solution. The exact solution of beam deflection is as follows:

x | t | Numerical solution | Exact solution | Absolute errors Rational interpolation | Absolute errors Lagrange interpolation |
---|---|---|---|---|---|

0.0990 | 0.0251 | −0.0019 | −0.0018 | 1.1803E−12 | 5.3171E−13 |

0.7939 | 1.0000 | −0.0315 | −0.0315 | 2.2247E−13 | 2.0698E−13 |

0.1654 | 1.9010 | 0.0651 | 0.0651 | 8.7291E−15 | 4.4977E−12 |

u ( x , t ) = sin ( 4 π x ) 3 t 2 exp ( − t ) . (22)

By the proposed algorithm, we obtain the numerical solutions which are given in Figures 4-6 and

Example 2 We consider a cantilever Euler-Bernoulli beam equation.

{ ∂ 2 ∂ x 2 ( E I ( x ) ∂ 2 u ( x , t ) ∂ x 2 ) + ∂ 2 u ( x , t ) ∂ t 2 = f ( x , t ) , u ( x , 0 ) = u t ( x , 0 ) = 0 , x ∈ [ 0 , 1 ] , u ( 0 , t ) = u x ( 0 , t ) = E I ( x ) ∂ 2 u ( 1 , t ) ∂ x 2 = 0 , ∂ ∂ x ( E I ( x ) ∂ 2 u ∂ x 2 ) ( 1 , t ) = − 6 π 3 t 2 e − t , t ≥ 0 , (23)

By the proposed algorithm, we obtain the numerical solution and the absolute error which are given in Figures 7-14 and

Example 3 Consider the model (3) with d 0 ( x ) = x , d 1 ( x ) = d 2 ( x ) = d 3 ( x ) = 0 , d 4 ( x ) = E I ( x ) = 1 + 1 4 exp ( − 40 ( x − 1 3 2 ) ) , d 5 ( x ) = x , d 6 ( x ) = sin x , α 1 = 0 , α 2 = 0 , α 4 = 0 , f ( x , t , u ) = g ( x , t ) + u 2 . The exact solution for this problem is

u ( x , t ) = 4 3 e 1 2 x − 2 3 t .

The numerical solution and the absolute error diagram of Example 3 are given in

x | t | Numerical solution | Exact solution | Absolute errors Rational interpolation | Absolute errors Lagrange interpolation |
---|---|---|---|---|---|

0.1284 | 0.3765 | 0.0058 | 0.0058 | 2.3209E−11 | 2.2052E−10 |

0.5522 | 1.6234 | 0.4990 | 0.4990 | 3.5232E−10 | 4.2077E−09 |

0.8715 | 1.9749 | 0.0327 | 0.0327 | 1.6892E−10 | 6.1629E−10 |

x | t | Numerical solution | Exact solution | Absolute errors Rational interpolation | Absolute errors Lagrange interpolation |
---|---|---|---|---|---|

0.5000 | 0.5661 | 0.1819 | 0.1819 | 1.5702E−11 | 4.72800E−11 |

0.6545 | 1.0000 | 0.2545 | 0.2545 | 1.4067E−11 | 3.0668E−10 |

0.9567 | 2.0000 | 0.0013 | 0.0013 | 2.2005E−09 | 1.6833E−08 |

x | t | Numerical solution | Exact solution | Absolute errors Present method ( M = 8 , N = 10 ) | Absolute errors Variational iteration method |
---|---|---|---|---|---|

0.5 | 0.2 | 1.4983 | 1.4983 | 0.0040E−05 | 0.004732E−02 |

0.5 | 0.4 | 1.3113 | 1.3113 | 0.0249E−05 | 0.006379E−02 |

0.5 | 0.8 | 1.0044 | 1.0044 | 0.1163E−05 | 0.013248E−02 |

0.5 | 1.0 | 0.8790 | 0.8790 | 0.1743E−05 | 0.000969E−02 |

1.0 | 0.2 | 1.9239 | 1.9239 | 0.0025E−05 | 0.006076E−02 |

1.0 | 0.6 | 1.4736 | 1.4736 | 0.0392E−05 | 0.019163E−02 |

1.0 | 0.8 | 1.2896 | 1.2896 | 0.0781E−05 | 0.017012E−02 |

1.5 | 0.2 | 2.4703 | 2.4703 | 0.0008E−05 | 0.007802E−02 |

1.5 | 0.4 | 2.1620 | 2.1620 | 0.0065E−05 | 0.010516E−02 |

1.5 | 0.8 | 1.6559 | 1.6559 | 0.0359E−05 | 0.021844E−02 |

1.5 | 1.0 | 1.4492 | 1.4492 | 0.0601E−05 | 0.001597E−02 |

In this paper, we use the method of barycentric Lagrange interpolation collocation method to solve Euler-Bernoulli beam equation. Some calculation results of this method and comparison with other methods show that this method has high accuracy and convergence. From this article, we can find that our method can be applied to solving such demographic models. So, we can extend this method to a winder area in the future.

The authors declare no conflicts of interest regarding the publication of this paper.

Zhang, H.L., Chen, L.W. and Fu, L. (2021) Numerical Solution of Euler-Bernoulli Beam Equation by Using Barycentric Lagrange Interpolation Collocation Method. Journal of Applied Mathematics and Physics, 9, 594-605. https://doi.org/10.4236/jamp.2021.94043