Numerical Solution of Euler-Bernoulli Beam Equation by Using Barycentric Lagrange Interpolation Collocation Method

Haolu Zhang¹, Lianwang Chen^2*, Lei Fu²

¹School of Civil Engineering, Inner Mongolia University of Technology, Hohhot, China.
²College of Science, Inner Mongolia University of Technology, Hohhot, China.

DOI: 10.4236/jamp.2021.94043                

Abstract

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.

Keywords

Barycentric Interpolation Collocation Method, Euler-Bernoulli Beam Equation, Linearized Iterative

1. Mathematical Modeling

The Mathematical Model of Euler-Bernoulli Beam Equation

Under the free vibration, the equation of transverse motion of a uniform Euler-Bernoulli beam is determined by a partial differential equation [1] [2], as shown in

$EI\frac{{\partial }^{4}u\left(x\mathrm{,}t\right)}{\partial x^4}+m\frac{{\partial }^{2}u\left(x\mathrm{,}t\right)}{\partial t^2}+r_e\frac{\partial u\left(x\mathrm{,}t\right)}{\partial t}+r_i\frac{\partial }{\partial t}\left(\frac{{\partial }^{4}u\left(x\mathrm{,}t\right)}{\partial x^4}\right)=0$ (1)

where $EI\mathrm{,}m\mathrm{,}{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 [3], respectively. The transverse displacement $u\left(x\mathrm{,}t\right)$ of the beam varies with the position x and time t.

In [4], define $u\left(x\mathrm{,}t\right)$ as the transverse deflection of the beam, $f\left(x\mathrm{,}t\right)$ as the generic arbitrary dynamic loads that distribute along the beam axis and t is time. The axial deformation is not considered here, namely the axial deformation is assumed to be zero. Using the Hamilton principle and employing the Euler-Bernoulli beam theory, variable coefficient Euler-Bernoulli beam equation can be obtained as

$EI\left(x\right)\frac{{\partial }^{4}u\left(x\mathrm{,}t\right)}{\partial x^4}+m\left(x\right)\frac{{\partial }^{2}u\left(x\mathrm{,}t\right)}{\partial t^2}+r_e\left(x\right)\frac{\partial u\left(x\mathrm{,}t\right)}{\partial t}+r_i\left(x\right)u\left(x\mathrm{,}t\right)=f\left(x\mathrm{,}t\right)$ (2)

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

$\begin{array}{l}{\displaystyle \underset{i=0}{\overset{4}{\sum }}}{d}_i\left(x\right)\frac{{\partial }^{i}u\left(x\mathrm{,}t\right)}{\partial x^i}+d_5\left(x\right)\frac{{\partial }^{2}u\left(x\mathrm{,}t\right)}{\partial t^2}+d_6\left(x\right)\frac{\partial u\left(x\mathrm{,}t\right)}{\partial t}\\ +d_7\left(x\right)\frac{\partial }{\partial t}\left(\frac{{\partial }^{4}u\left(x\mathrm{,}t\right)}{\partial x^4}\right)=f\left(x,t,u\right),\end{array}$ (3)

the initial conditions are expressed as:

$u\left(x,0\right)=\varphi \left(x\right),u_t\left(x,0\right)=\psi \left(x\right),x\in \left[0,l\right],$ (4)

and the boundary conditions are expressed as:

$u\left(0,t\right)=p_0\left(t\right),u\left(l,t\right)=p_1\left(t\right),u_x\left(0,t\right)=p_2\left(t\right),u_x\left(l,t\right)=p_3\left(t\right),t\ge 0,$ (5)

where, $u\left(x\mathrm{,}t\right)$ represents an unknown function at position x and time t. $d_i\left(x\right),i=0,1,\cdots ,7$ is a known function. If $f\left(x\mathrm{,}t\mathrm{,}u\right)$ is nonlinear of u, the equation is a nonlinear Partial differential equation. If $f\left(u\right)$ is linear of u, then the Equation (3) is linear Partial differential equation.

2. The Barycentric Lagrange Interpolation Collocation Method

2.1. In Model (1), If $f\left(x,t,u\right)=f\left(x,t\right)$, the Model Is a Linear Problem

The model (3) can be written as:

$\begin{array}{l}{\displaystyle \underset{i=0}{\overset{4}{\sum }}}{d}_i\left(x\right)\frac{{\partial }^{i}u\left(x,t\right)}{\partial x^i}+d_5\left(x\right)\frac{{\partial }^{2}u\left(x,t\right)}{\partial t^2}+d_6\left(x\right)\frac{\partial u\left(x,t\right)}{\partial t}\\ +d_7\left(x\right)\frac{\partial }{\partial t}\left(\frac{{\partial }^{4}u\left(x,t\right)}{\partial x^4}\right)=f\left(x,t\right),\end{array}$ (6)

We introduce the barycentric Lagrange interpolation collocation method [5] - [10] to solve this problem:

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

$x^0={\left[x_1,x_2,\cdots ,x_M\right]}^{\text{T}},t^0={\left[t_1,t_2,\cdots ,t_N\right]}^{\text{T}},$ (7)

then the matrixs X, T above the region $\Omega$ can be generated by the above column vectors $x^0\mathrm{,}{t}^0$, and the form of matrixs X and T are as follows:

$X^{\text{T}}=\left[\begin{array}{c}{\left(x^0\right)}^{\text{T}},{\left(x^0\right)}^{\text{T}},\cdots ,{\left(x^0\right)}^{\text{T}}\end{array}\right],T=\left[t^0,t^0,\cdots ,t^0\right].$ (8)

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

$x={\left[X_1,X_2,\cdots ,X_{M\times N}\right]}^{\text{T}},t={\left[T_1,T_2,\cdots ,T_{M\times N}\right]}^{\text{T}}.$ (9)

The relations between the components $x_i\mathrm{,}{y}_j$ of the Formula (11) and the components $X_K\mathrm{,}{Y}_K$ of Formula (13) are as follows:

$\begin{array}{l}{X}_k=X_{\left(i-1\right)N+j}=x_i,T_k=T_{\left(i-1\right)N+j}=t_j,\\ i=1,2,\cdots ,M;j=1,2,\cdots ,N;k=1,2,\cdots ,M\times N.\end{array}$ (10)

The barycentric interpolation of $u\left(x\mathrm{,}t\right)$ at nodes $\left(x_i\mathrm{,}{t}_j\right)$ can be written as [5] - [10]:

$u\left(x,t\right)={\displaystyle \underset{i=1}{\overset{M}{\sum }}}{\displaystyle \underset{j=1}{\overset{N}{\sum }}}{\xi }_i\left(x\right){\eta }_j\left(t\right)v_{ij},i=1,2,\cdots ,M;j=1,2,\cdots ,N,$ (11)

where, ${\xi }_i\left(x\right)=\frac{{\displaystyle \underset{k=1,k\ne i}{\overset{M}{\prod }}}\left(x-x_k\right)}{{\displaystyle \underset{k=1,k\ne i}{\overset{M}{\prod }}}\left(x_i-x_k\right)},i=1,2,\cdots ,M$, ${\eta }_j\left(t\right)=\frac{{\displaystyle \underset{k=1,k\ne j}{\overset{N}{\prod }}}\left(t-t_k\right)}{{\displaystyle \underset{k=1,k\ne j}{\overset{N}{\prod }}}\left(t_j-t_k\right)},j=1,2,\cdots ,N$.

Use Formula (11), the $l+k$ order partial derivative of function $v\left(x\mathrm{,}t\right)$ at nodes can be expressed as:

$\frac{{\partial }^{l+k}u}{\partial x^l\partial t^k}={\displaystyle \underset{i=1}{\overset{M}{\sum }}}{\displaystyle \underset{j=1}{\overset{N}{\sum }}}{\xi }_i^{\left(l\right)}\left(x\right){\eta }_j^{\left(k\right)}\left(t\right)u_{ij},l,k=1,2,\cdots .$ (12)

At nodes $\left(x_p\mathrm{,}{t}_q\right)$, the function values of partial derivative are defined as:

$\begin{array}{l}{u}^{\left(l,k\right)}\left(x_p,t_q\right):=\frac{{\partial }^{l+k}u\left(x_p,t_q\right)}{\partial x^l\partial t^k}\\ ={\displaystyle \underset{i=1}{\overset{M}{\sum }}}{\displaystyle \underset{j=1}{\overset{N}{\sum }}}{\xi }_i^{\left(l\right)}\left(x_p\right){\eta }_j^{\left(k\right)}\left(t_q\right)u_{pq},p=1,2,\cdots ,M;q=1,2,\cdots ,N.\end{array}$ (13)

The function values of the Formula (11) and the Formula (12) at the node form column vectors $u\mathrm{,}{u}^{\left(l\mathrm{,}k\right)}$, and they are as follows:

$u={\left[u_1,u_2,\cdots ,u_{M\times N}\right]}^{\text{T}}$, $u^{\left(l,k\right)}={\left[u_1^{\left(l,k\right)},u_2^{\left(l,k\right)},\cdots ,u_{M\times N}^{\left(l,k\right)}\right]}^{\text{T}}$,

$u_p=u\left(X_p,T_p\right)$, $u_p^{\left(l,k\right)}=u^{\left(l,k\right)}\left(X_p,T_p\right)$, $p=1,2,\cdots ,M\times N$.

(13) can be written in following matrix form [5] - [10]:

$u^{\left(l\mathrm{,}k\right)}=D^{\left(l\mathrm{,}k\right)}u\mathrm{.}$ (14)

In the above formula, $D^{\left(l\mathrm{,}k\right)}=C^{\left(l\right)}\otimes D^{\left(k\right)}$ is the Kronecker product of matrix $C^{\left(l\right)}$ and $D^{\left(k\right)}$, which is also called $\left(l\mathrm{,}k\right)$ order partial differential matrix at nodes $\left\{\left(x_i,t_j\right),i=1,2,\cdots ,M;j=1,2,\cdots ,N\right\}$. $C^{\left(l\right)}$ is l order differential matrix on x direction nodes, and $D^{\left(k\right)}$ is k order differential matrix on t direction nodes.

Denote:

$C^{\left(0\right)}=I_M,D^{\left(0\right)}=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:

$\left[{\displaystyle \underset{i=0}{\overset{4}{\sum }}}diag\left(d_i\left(x\right)\right)D^{\left(i\mathrm{,0}\right)}+{\displaystyle \underset{i=1}{\overset{2}{\sum }}}diag\left(d_{4+i}\left(x\right)\right)D^{\left(\mathrm{0,}i\right)}+diag\left(d_7\left(x\right)\right)D^{\left(\mathrm{4,1}\right)}\right]u=f\mathrm{.}$ (16)

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

2.2. If the Model Is a Nonlinear Problem, We Can Use the Solution Method of the Linear Problem Above

Given initial value $v_0$, we can construct following linear iterative format:

$\begin{array}{l}{\displaystyle \underset{i=0}{\overset{4}{\sum }}}{d}_i\left(x\right)\frac{{\partial }^i{u}_n\left(x,t\right)}{\partial x^i}+d_5\left(x\right)\frac{{\partial }^2{u}_n\left(x,t\right)}{\partial t^2}+d_6\left(x\right)\frac{\partial u_n\left(x,t\right)}{\partial t}\\ +d_7\left(x\right)\frac{\partial }{\partial t}\left(\frac{{\partial }^4{u}_n\left(x,t\right)}{\partial x^4}\right)=f\left(x,t,u_{n-1}\right),\end{array}$ (17)

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

$\left[{\displaystyle \underset{i=0}{\overset{4}{\sum }}}diag\left(d_i\left(x\right)\right)D^{\left(i\mathrm{,0}\right)}+{\displaystyle \underset{i=1}{\overset{2}{\sum }}}diag\left(d_{4+i}\left(x\right)\right)D^{\left(\mathrm{0,}i\right)}+diag\left(d_7\left(x\right)\right)D^{\left(\mathrm{4,1}\right)}\right]u=f_0\mathrm{.}$ (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:

$\begin{array}{l}{\displaystyle \underset{i=1}{\overset{M}{\sum }}}{\displaystyle \underset{j=1}{\overset{N}{\sum }}}{\xi }_i\left(x_p\right){\eta }_j\left(0\right)u_{ij}=\varphi \left(x_p\right),{\displaystyle \underset{i=1}{\overset{M}{\sum }}}{\displaystyle \underset{j=1}{\overset{N}{\sum }}}{\xi }_i\left(x_p\right){{\eta }^{\prime }}_j\left(0\right)u_{ij}=\phi \left(x_p\right),\\ {\displaystyle \underset{i=1}{\overset{M}{\sum }}}{\displaystyle \underset{j=1}{\overset{N}{\sum }}}{\xi }_i\left(0\right){\eta }_j\left(t_q\right)u_{ij}=p_0\left(t_q\right),{\displaystyle \underset{i=1}{\overset{M}{\sum }}}{\displaystyle \underset{j=1}{\overset{N}{\sum }}}{\xi }_i\left(l\right){\eta }_j\left(t_q\right)u_{ij}=p_1\left(t_q\right),\\ {\displaystyle \underset{i=1}{\overset{M}{\sum }}}{\displaystyle \underset{j=1}{\overset{N}{\sum }}}{{\xi }^{\prime }}_i\left(0\right){\eta }_j\left(t_q\right)u_{ij}=p_2\left(t_q\right),{\displaystyle \underset{i=1}{\overset{M}{\sum }}}{\displaystyle \underset{j=1}{\overset{N}{\sum }}}{{\xi }^{\prime }}_i\left(l\right){\eta }_j\left(t_q\right)u_{ij}=p_3\left(t_q\right).\end{array}$ (19)

The boundary conditions are replaced by displacement method.

3. Numerical Simulations

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.

$\{\begin{array}{l}\frac{{\partial }^2}{\partial x^2}\left(EI\left(x\right)\frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}\right)+\frac{{\partial }^{2}u\left(x,t\right)}{\partial t^2}=f\left(x,t\right),\\ u\left(x,0\right)=u_t\left(x,0\right)=0,x\in \left[0,1\right],\\ u\left(0,t\right)=u\left(1,t\right)=u_x\left(0,t\right)=u_x\left(1,t\right)=0,t\ge 0,\end{array}$ (20)

1) where $EI\left(x\right)=1+\sin \pi x$, and the source term $f\left(x\mathrm{,}t\right)$ is determined by (21) consistent with the chosen solution. The exact solution of beam deflection is as follows:

$u\left(x,t\right)=x\left(1-x\right)\sin \left(4\pi x\right)t^2{\text{e}}^{-t}.$ (21)

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

2) where $EI\left(x\right)=1+\frac{1}{4}\exp \left(-40\left(x-{\frac{1}{3}}^2\right)\right)$, and the source term $f\left(x\mathrm{,}t\right)$ is

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

Figure 1. Exact solution for Example 1 (i).

Figure 2. Numerical solution for Example 1 (i).

Figure 3. Absolute errors for Example 1 (i).

Table 1. Numerical comparison of absolute errors of Example 1 (i) with M = 31, N = 15.

$u\left(x,t\right)=\sin {\left(4\pi x\right)}^3{t}^2\exp \left(-t\right).$ (22)

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

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

$\{\begin{array}{l}\frac{{\partial }^2}{\partial x^2}\left(EI\left(x\right)\frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}\right)+\frac{{\partial }^{2}u\left(x,t\right)}{\partial t^2}=f\left(x,t\right),\hfill \\ u\left(x,0\right)=u_t\left(x,0\right)=0,x\in \left[0,1\right],\hfill \\ \begin{array}{l}u\left(0,t\right)=u_x\left(0,t\right)=EI\left(x\right)\frac{{\partial }^{2}u\left(1,t\right)}{\partial x^2}=0,\\ \frac{\partial }{\partial x}\left(EI\left(x\right)\frac{{\partial }^{2}u}{\partial x^2}\right)\left(1,t\right)=-6{\pi }^3{t}^2{\text{e}}^{-t},t\ge 0,\end{array}\hfill \end{array}$ (23)

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

Example 3 Consider the model (3) with $d_0\left(x\right)=x$, $d_1\left(x\right)=d_2\left(x\right)=d_3\left(x\right)=0$, $d_4\left(x\right)=EI\left(x\right)=1+\frac{1}{4}\exp \left(-40\left(x-{\frac{1}{3}}^2\right)\right)$, $d_5\left(x\right)=x$, $d_6\left(x\right)=\sin x$, ${\alpha }_1=0$, ${\alpha }_2=0$, ${\alpha }_4=0$, $f\left(x,t,u\right)=g\left(x,t\right)+u^2$. The exact solution for this problem is

$u\left(x,t\right)=\frac{4}{3}{\text{e}}^{\frac{1}{2}x-\frac{2}{3}t}.$

The numerical solution and the absolute error diagram of Example 3 are given in Figure 15 and Figure 16 and Table 4, respectively. As can be seen from the Figure, the method can be used to obtain a smaller absolute error.

Figure 4. Exact solution for Example 1 (ii).

Figure 5. Numerical solution for Example 1 (ii).

Figure 6. Absolute errors for Example 1 (ii).

Table 2. Numerical comparison of absolute errors of Example 1 (ii) with M = 31, N = 15.

Figure 7. Exact solution for Example 2.

Figure 8. Numerical solution for Example 2.

Figure 9. Absolute errors by using Lagrange interpolation at $M=21,N=15$ for Example 2.

Table 3. Numerical comparison of absolute errors of Example 2 with M = 31, N = 15.

Figure 10. Absolute errors by using Lagrange interpolation at $M=31,N=15$ for Example 2.

Figure 11. Absolute errors by using Lagrange interpolation at $M=41,N=15$ for Example 2.

Figure 12. Absolute errors by using rational interpolation at $M=21,N=15$ for Example 2.

Figure 13. Absolute errors by using rational interpolation at $M=31,N=15$ for Example 2.

Figure 14. Absolute errors by using rational interpolation at $M=41,N=15$ for Example 2.

Figure 15. Numerical solution obtained by the present method for Experiment 3 with $M=50,N=50$.

Figure 16. Absolute error for Experiment 3 with $M=8,N=10$.

Table 4. Comparison of the absolute errors for Example 3.

4. Conclusion

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.

Conflicts of Interest

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

References

[1]	Rao, S.S. (2019) Vibration of Continuous Systems. John Wiley Sons, Inc., Hoboken.
[2]	Clough, R. and Penzien, J. (1993) Dynamics of Structures. McGraw-Hill, New York.
[3]	Banks, H. and Inman, D. (1991) On Damping Mechanisms in Beams. Journal of Applied Mechanics, 58, 716-723. https://doi.org/10.1115/1.2897253
[4]	Yu, H.T., Yang, Y.S. and Yuan, Y. (2018) Analytical Solution for a Finite Euler-Bernoulli Beam with Single Discontinuity in Section under Arbitrary Dynamic Loads. Applied Mathematics Model, 60, 571-580. https://doi.org/10.1016/j.apm.2018.03.046
[5]	Li, S.P. and Wang, Z.Q. (2013) Barycentric Interpolation Collocation Method for Solving Elastic Problems. Journal of Central South University, 44, 2031-2040.
[6]	Li, S.P. and Wang, Z.Q. (2018) Barycentric Interpolation Collocation Method for Solving Nonlinear Vibration Problems. Noise Vibration Control, 28, 49-52.
[7]	Li, S.P. and Wang, Z.Q. (2015) Barycentric Interpolation Collocation Method for Nonlinear Problems. National Defense Industry Press, Beijing.
[8]	Li, S.P. and Wang, Z.Q. (2012) High-Precision Non-Grid Center of Gravity Interpolation Collocation Method: Algorithm, Program and Engineering Application. Science Press, Beijing.
[9]	Liu, H.Y., Huang, J., Pan, Y.B. and Zhang, J.P. (2018) Barycentric Interpolation Collocation Methods for Solving Linear and Nonlinear High-Dimensional Fredholm Integral Equations. Computers Mathematics with Applications, 327, 141-154. https://doi.org/10.1016/j.cam.2017.06.004
[10]	Liu, F.F., Wang, Y.L. and Li, S.G. (2018) Barycentric Interpolation Collocation Method for Solving the Coupled Viscous Burgers’ Equations. International Journal of Computer Mathematics, 95, 2162-2173. https://doi.org/10.1080/00207160.2017.1384546