A fast accurate artificial boundary condition for the Euler-Bernoulli beam

It is difficult to convert the PDEs defined on the whole space with higher order space derivatives into initial boundary value problems by proposing boundary conditions. In this paper, we use an absorbing boundary condition method to solve the Cauchy problem for one-dimensional Euler-Bernoulli beam with fast convolution boundary condition which is derived through the Padé approximation for the square root function ⋅\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\sqrt {\cdot }$\end{document}. We also introduce a constant damping term to control the error between the resulting approximation of the Euler-Bernoulli system and the original one. Numerical examples verify the theoretical results and demonstrate the performance for the fast numerical method.


Introduction
As one of the elemental objects in solid mechanics, the Euler-Bernoulli beam has wide applications, such as in railway and bridge engineering [1,2]. The infinite beam models still exist in some special problems, for instance, the development of multiscale algorithms for beams [3] and the study of beams response subjected to moving load [4]. The following shows the infinite Euler-Bernoulli beam vibration equation without loads: where , , , and are the material density, the area of the beam cross-section, the Young's modulus, and the bending moment of inertia, respectively. In (1), represents the transverse displacement of an undriven beam.
The original system (1) is derived and set on the whole space. Therefore, the numerical study for this kind of system should be restricted to a bounded domain and one has to propose transparent (i.e., non-reflective) boundary conditions at the boundaries of the bounded domain. The problem can be set as follows: given compactly supported initial data, one searches for transparent (i.e., non-reflective) boundary conditions such that the solution computed with these boundary conditions coincides with the real solution. There are rather extensive studies [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22] on ABCs (artificial boundary condition) for heat equation, wave equation, and classical Schrödinger equation. In general, for the Cauchy problems of partial differential equation, the exact ABCs are nonlocal in time, having some temporal convolutions in the formulations. The nonlocal convolutions in the exact ABCs cause huge amount of calculation for the truncated problem in bounded domains (the computational cost for the convolution integral at the th time step is of order ). Thus, the fast algorithms are introduced to reduce the computing load. In [29,30], the matching boundary conditions (MBC) for the semi-discrete version of Euler-Bernoulli beam (1) with compact supported initial data are proposed. However, when using MBC which enable faster computations, the issue of accuracy comes up especially for the PDEs with higher order space derivatives.
In order to treat the accuracy and the efficiency issues, we propose a fast algorithm to solve the Cauchy problem for one-dimensional Euler-Bernoulli beam (1) with compact supported initial data. To this end, firstly we convert the Euler-Bernoulli beam (1) into an equivalent form by introducing a constant damping term. After that, a semi-discrete exact boundary condition for the semi-discrete system is derived through the Laplace-transform, and then a second-order time discretization is used to give a fully discrete scheme for the semi-discrete system. Using the Padé rational expansion of the square root function [27], the semi-discrete convolution kernel involved in the exact semi-discrete boundary conditions is approximated by the summation of exponentials which enable a fast convolution computation [23][24][25][26][27][28] at the boundaries. Then the computing load can be reduced from to 1 . Some kinds of approximations for kernel symbols are studied in [23][24][25][26][27][28]31]. Finally, we give some numerical tests. The proposed method is not only more accurate, it is more important that we can measure the accuracy of numerical results through numerical analysis.
The rest of this paper is organized as follows. In Section 2, we introduce the setting of the problem and convert the system into a semi-discrete system. The exact ABC for the semi-discrete system is obtained. In Section 3, a Padé approximation is introduced to approximate the convolution kernel in the exact semi-discrete ABCs, which brings a fast convolution ABC. In Section 4, time discretization is given to obtain a fully discrete scheme for the semi-discrete system. In Section 5, we estimate the error between the resulting approximation system and original system. In Section 6, numerical examples are given to demonstrate the effectiveness of the proposed numerical method.

Semi-discretization for the Euler-Bernoulli beam and its exact ABC
In this section, a four-order semi-discrete scheme in space for one-dimensional Euler-Bernoulli beam will be proposed and an exact semi-discrete ABC for the semi-discrete system will be derived.
It is convention to define Laplace transform for given function , namely, 0 where is in the right half complex plane . We consider the initial value problem for the Euler-Bernoulli beam (1) set on the whole space with constant , , , and . We choose to denote the spatial step size. Then, the mesh points can be defined as for . We then use to denote the semi-discrete numerical approximation of in the original system (1). One common spatially discrete model of the beam equation in application is the central finite-difference scheme, Thus, the approximation system can be rewritten as Let us confine the computational domain of (1) as rather than the whole real axis, we assume to be a positive integer such that with the mesh size . Accordingly, the initial data is compactly supported, i.e., 0 0 and 0 0 hold for for system (2). Then, the truncation system of (2) is defined with . Thus, we still need the values of 1 and 2 to close the truncation system. Firstly, we derive the right boundary condition. Due to the compactly supported assumption, after a Laplace transform for the control equation of (2) with , one has ( 3 ) One assumes We can pick two roots 1 and 2 such that 1 0 and 2 0 as . Then, by (4), we get Thus, the right boundary condition can be written as 1 1 . The left boundary condition can be obtained in the same manner.
Assuming the initial data is compactly supported such that 0 0 and 0 0 for , by (5), we can reduce the semi-discretized problem (2) to where is the Bessel function of the order . We refer to [32] where more details for the derivations of and can be found. We introduce and 1 such that, 1 1 and . By introducing 1 and 2 such that 1 1 and 2 1 2 , it is obvious that In the same manner, one has Finally, by (7) and (8), the semi-discrete reduced problem (6) can be reformulated to the following semi-discretized problem on a bounded domain:

Reformulation of (9)
We introduce a transform such that . In this case, (9) can be written as where 1 1 and . We will see in the next several sections this transform allows us to approximate the transformed kernel functions with high accuracy.
In this section, an algorithm for approximating the convolution , 1 , and 2 in (10) to (15) is introduced. Then, we give a fast convolution absorbing boundary condition for system (10) to (15).

Rational approximation for the symbol of the convolution kernel
In [27], the Pade approximation for the function 1 can be expressed as Using this Pade approximation, we can write a rational approximation for the square root function as: Thus, Then the rational approximation of the above symbol can be written as: We can approximate the convolution kernel 1 and 2 by 1 and 2 in the same manner as (17).

Fast convolutions of (21) to (24)
In Section 2.1, we introduced a transform such that the original system can be converted into an approximation system (19) to (24) by Pade approximation. Recalling (17), we can use fast algorithm introduced in the previous section to compute (21) to (24), namely, with 0 0. For two sequences and , we use Thus, one has the following recurrence relation, Therefore, the numerical approximation of can be implemented by a fast convolution (29) in (25), namely, 1 1 (30) with the numerical solution of . The numerical approximations of 1 and 2 can be obtained in the same manner. Therefor, the boundary integrals (21) to (24) can be evaluated by fast convolution performed in (30).
From (28) and (30), we can see the computational cost in the th time step is of order rather than . Thus, the computational cost for all the time steps is order rather than 2 .

Numerical discretization for control equation
We assume is the total computational time and is the total number of steps, and the time step. Then we can define the time points as following: 0 1 . We can rewrite the above scheme as it is the interior scheme of 1 . This time discretization scheme is equivalent to the implicit Newmark method with integral parameters being 1 2 and 1 4 .

Numerical discretization combining with boundary condition
On the boundary, (21) can be discretized as In the same way, the equivalent discretized version of (24) reads

Error between approximation system and original system
We estimate the error between approximation system (19) to (24) and original system (10) to (15) in the following theorem. where is a constant depending on , , , 0 , and 1 .
In order to prove Theorem 1, we will give some useful lemmas and propositions. Firstly we introduce the notations such that In order to estimate the errors between the kernel functions , 1 , 2 and the approximation kernel functions , 1 , 2 , we give the following proposition. Recalling defined in (16), the following error result can be found in [27].

Lemma 4 Let 0.
Then we have the following identity: Then we prove (41).
In the same way, if 0 and satisfies (47), . Thus, we prove (44). Then, by Proposition 1, we can prove the following lemma.

Lemma 5
Under assumption of Proposition 1 with 0 0, for , we have the following estimates, Hence, we prove (49). In the same manner, (50) to (53) can be proved through Proposition 1.
Let us prove the Theorem 1 by using Lemma 5.
Proof We write the solution of (19) to (24) as , such that We can define error function as 1 and 2 where , is the solution of (10) to (15). The system for 1 and 2 reads, which finishes the proof of Theorem 1.
Remark 1 Theorem 1 implies that by choosing appropriate , the semi-discrete system (10) to (15) can be approximated by system (19) to (24) with arbitrary accuracy.

Numerical results
Some numerical tests are provided to validate the theoretical results in this section. In the tests, we will take 0.01. Following Theorem 1, we will also take the number of Padé expansion terms by using the following criterion: where is given in Lemma 3. If is the total number of time step, the computation cost of the fast convolution is of order rather than 2 .

Example 1
Firstly, the same initial deflection distribution as in [29] is considered, and we take velocity as zero distribution, namely, The computation domain is also identical to those in [29], i.e., 20 and the computational time 1.5. By taking 0.05, we can convert the original system (1) into the semidiscretized system (10) to (15). Then, the numerical solution can be obtained based on the algorithm proposed in Section 4 with 10 4 . The evolution of numerical solutions for the deflection is depicted in Fig. 1. Only the positive half of the spatial domain is shown for the symmetry of the problem. It can be seen the right going wave propagates with an almost constant group velocity, and there is no significant reflection when it meets the boundary. A reference solution is computed on a much larger spatial domain to avoid any boundary effect by taking 500, and the difference between the two is shown in Fig. 2. Although the error caused by the artificial reflection can be noticed, its magnitude is reasonably small.
In the same manner, the evolution of numerical solutions for the transversal velocity and the errors are depicted in Figs. 3 and 4. We do not detect significant spurious reflection near the boundaries.
To study the performance of the proposed method, we fix the total computing time , but employ different time steps , ranging from 10 5 to 10 3 . We report the -norm errors of at different time for various values of in Table 1 and the -norm errors of in Table 2. The -norm of the numerical errors of the numerical solutions at 1.5 are then plotted against the time step in the left panel of Fig. 5. From the log-log plot, it can be seen the data points are almost on a straight line with a slope of 2. This indicates a second-order convergence order of the proposed method. For a closer look at the solution, the deflection at 1.5 is shown in Fig. 6. The solution obtained by using matching boundary condition (MBC, see [29]) is also given for comparison. At 1.5, when main waves should have left the computation domain, as the MBC solution suffers from significant artificial reflection.

Example 2
In order to observe more clearly the dispersion behavior of the system, we further consider an initial velocity problem, namely, 0 0 1 exp 2 2 cos 4 5 0 5.     The computation domain, the settings of parameters are identical to those in Example 1. The deflection evolution is shown in Fig. 7. The dispersion phenomenon becomes more obvious as the wave package propagates, but no spurious reflection is found near the boundary. The error between the numerical solution and the reference one is shown in Fig. 8. The error becomes greater with increasing time, but the magnitude is still small at the end of the computation.
In the same manner, the evolution of transversal velocity and the errors are depicted in Figs. 9 and 10. The reflection near the boundaries is also very weak.
We report the -norm errors of at different time for various values of in Table 3 and the -norm errors of in Table 4. Different time steps are employed, and the numerical errors of the displacement and velocity vectors at time 1.5 are plotted in the left panel of Fig. 11. We can see a second-order convergence order for -norm, again. The right panel of Fig. 11 clearly shows that the method has a linear time complexity.

Conclusion
We converted the Euler-Bernoulli beam on a unbounded domain into a semidiscretized form, then we use a two-order numerical scheme to solve the resulted initial-boundary value problem with fast convolutions at the boundaries where the convolution kernel is approximated by the summation of exponentials through the Padé rational expansion. A damping term was introduced such that it can be proved theoretically that the error between the resulting approximation system and original system can be arbitrarily small by choosing appropriate terms of exponentials. The theoretical analysis and the effectiveness of the numerical method are confirmed by the numerical tests.