Green functions for three-point boundary value problems governed by differential equation systems with applications to Timoshenko beams

The present paper is devoted to the issue of the Green function matrices that belongs to some three-point boundary- and eigenvalue problems. A detailed definition is given for the Green function matrices provided that the considered boundary value problems are governed by a class of ordinary differential equation systems associated with homogeneous boundary and continuity conditions. The definition is a constructive one, i.e., it provides the means needed for calculating the Green function matrices. The fundamental properties of the Green function matrices—existence, symmetry properties, etc.—are also clarified. Making use of these Green functions, a class of three-point eigenvalue problems can be reduced to eigenvalue problems governed by homogeneous Fredholm integral equation systems. The applicability of the novel findings is demonstrated through a Timoshenko beam with three supports.

The influence of an axial load on the natural frequencies of a uniform single-span beam is investigated in [10]. It is found that Galef's formula is only valid for a few types of end-conditions. The effect of end-supports on the frequency is only significant for the first few modes. The vibratory behavior of clamped-free beams with an intermediate axial force is the subject of [11] using the Hamilton principle. The frequencies show an increase as the force is moved closed to the fixed end. It is studied [12] how a beam system with tendon loading vibrates. The effect of the number and location of attachment points is evaluated.
The dynamic response of beams to a moving load is sought in [13] by a linear model. The dynamic Green function is constructed in closed form for elastic end-restraints and with this; the results are found in an exact and direct method. The effect of the load speed parameter is studied on the deflections for various, limiting boundary conditions. Apart from vibrations, it is worthy to mention some further notable results within the field of the Green function [14]. Book [15] that presents the Green theorem, introduces the concept of the Green function and applies it to electrostatic problems. Since that pioneering work, the Green function has widely been applied to various problems [16,17]. The concept of the Green function for two-point boundary value problems governed by ordinary differential equations was first introduced in paper [18]. Furthermore, books [19,20] extend the knowledge by defining the Green function for ordinary linear differential equations with their most important properties. Later, the Green function matrix was introduced as a generalization for a class of ordinary differential equation systems in [21]. For degenerated ordinary differential equation systems, new findings are given in [22]. The existence for some three-point boundary value problems associated with third-order nonlinear differential equations is detailed in [23] using Green functions. For a class of second-order ordinary differential equations, a method is proposed to find the corresponding Green functions for three-point boundary value problems in [24]. Similarly, article [25] is about a special class of third-order three-point boundary value problems. The author demonstrates how to find the Green function for such issues. An application is also presented. A non-local three-point boundary value problem is selected and examined in [26]. Existence and uniqueness of solutions are given, and the Green function is also constructed. A type of nonlinear third-order non-local boundary value problem is in the spotlight in [27]. The Schauder fixed point theorem is used for the solution. A third-order linear differential equation is investigated in [28]. The existence of the Green function is proven, and solution is given. Article [29] determines the Green function for compressed straight Euler-Bernoulli beams with three supports and solves the linear stability issue by a boundary element technique. The procedure is applicable when the three-point boundary value problem is governed by a single ordinary differential equation.
Based on the above, to the best knowledge of the authors, the Green function matrix has not been defined for three-point boundary value problems governed by ordinary differential equation systems. To this end, the article addresses this issue with some numerical examples as possible and effective applications. First, the related class of three-point boundary value problems is introduced. Then, a suitable definition and computation steps for the corresponding Green function matrix are given. Numerical examples related to nonhomogeneous Timoshenko beams with three supports are given as one possible demonstration of the applicability of the findings.

Differential operator with boundary and continuity conditions
Consider the class of eigenvalue problems governed by the homogeneous ordinary differential equation system K y = λ M y , (2.1a) where y(x) = [y 1 (x)|y 2 (x)| . . . |y n (x)], n ≥ 2 is the unknown function vector, λ is a parameter (the eigenvalue sought), 2κ and 2μ (κ > μ) are the order of the differential operators K y and M y . Let K ν (x) (n×n) , (ν = 0, 1, . . . , κ) and M ν (x) (n×n) , (ν = 0, 1, . . . , μ) be square matrices in the interval x ∈ [a, c] (c > a). The differential operators K y and M y are defined by K y = κ ν=0 (−1) ν K ν (x) y (ν) (x) (ν) , d n (. . .) dx n = (. . .) (n) ; (ν) . (2.1b) It is assumed that (K ν (x)) [M ν (x)] is differentiable continuously(κ)[μ] times and K 2κ (x) and M 2μ (x) have inverse if x ∈ [a, c]. It is also assumed that they are differentiable as many times as required. Note that 2κ, which is the order of the differential operator on the left side of (2.1a), is greater than 2μ, which is the order of the differential operator on the right side. Let x ∈ [a, c], c > a, c − a = be the interval in which the solution of differential equation Some quantities in the intervals [a, b] and [b, c] are denoted by the Latin I and I I subscripts. Accordingly, y I and y I I are the solutions to the differential equation (2.1) in the intervals I and I I .
Differential equation system (2.1) is associated with the following boundary and continuity conditions: where α νr I , β νr I , β νr I I and γ νr I I are nonzero square matrices of size n × n. Differential equation (2.1) with the boundary and continuity conditions (2.2) constitute an eigenvalue problem with λ being the eigenvalue to be found. If μ = 0, the right side of (2.1) changes as and the eigenvalue problem is called simple. From now on, it is assumed that M 0 (x) has an inverse if x ∈ [a, c]. The eigenvalue problems that provide the eigenfrequencies for the longitudinal and torsional vibrations of rods as well as for the transverse vibrations of strings and beams are all simple ones. If μ > 0, the eigenvalue problem is called generalized [19].
The scalar equations that constitute the boundary-and continuity conditions should be linearly independent of each other. It is obvious that a linear combination of the boundary conditions is also a boundary condition. By selecting suitable linear combinations, derivatives with an order higher than κ − 1 could be removed. If it is done in all possible ways, the total number of boundary conditions which do not involve derivatives higher than κ − 1 is, say, e. These boundary conditions are called essential boundary conditions. The remaining 2k − e boundary conditions are the natural boundary conditions. The column matrices u(x) and v(x) (|u(x)|, |v(x)| are not identically equal to zero if x ∈ [a, c]) are called comparison functions if they satisfy the boundary and continuity conditions and are called eigenfunctions if they, in addition to this, satisfy differential equation (2.1).

Self-adjointness
The integrals taken on the set of the comparison functions u(x), v(x) are products defined on the operators K and M. Let us now detail the product (u, v) K . Making use of (2.1b), one may write in which the integral can be manipulated into a more suitable form through integrations: Hence, It follows from Eq. (2.8a) that These results are naturally valid for the products (v, u) K and (v, u) M . The expressions K 0 (u, v) and M 0 (u, v) defined by the right sides of Eq. (2.8) are called boundary-and continuity expressions. Eigenvalue problem (2.1), (2.2) is said to be self-adjoint if the products (2.4) are commutative, i.e., it holds that Conditions (2.9) are called conditions of self-adjointness.
Assume that the three-point boundary value problem (2.1), (2.2) is self-adjoint. Then, the eigenfunctions are orthogonal to each other in general sense: k, = 1, 2, 3 . . . (2.11) In addition to this, the eigenvalues are all real numbers. These statements can be proven easily by recalling the similar proofs given for two-point boundary value problems in [30].

Sign of eigenvalues
The three-point eigenvalue problem ( Hence, . (2.13) This equation shows that the sign of λ is a function of the products (y , y ) K and (y , y ) M . Assume that (u, u) K > 0, and (u, u) M > 0 (2.14) for any comparison function u(x). Then, the three-point eigenvalue problem (2.1), (2.2) is positive definite. The Rayleigh quotient is defined by the equation (2.15) in which u(x) is a comparison function. Substituting (2.8a) and (2.8b) for the numerator and denominator in (2.13) yields . (2.16) If the three-point eigenvalue problem (2.1), (2.2) is self-adjoint, then Assume that the three-point eigenvalue problem (2.1), (2.2) is self adjoint and Then, the considered self-adjoint three point eigenvalue problem is called simple.

Determination of the eigenvalues
Let us denote the linearly independent particular solutions of the differential equation K y = λ M y by z (x, λ) ( = 1, 2, . . . , 2κ × n). With z (x, λ), the general solution is of the form This equation needs to be solved to find the eigenvalues. The determinant Δ(λ) is referred to as characteristic determinant.
If Δ(λ) is identically equal to zero, then each λ is an eigenvalue. Otherwise, function Δ(λ) has an infinite sequence of isolated zero points which can be ordered according to their magnitude:

The Green function matrix
Consider the inhomogeneous ordinary differential equation system where the differential operator of order 2κ is defined by the following equation: Here, κ ≥ 1 is a natural number, , (ν = 0, 1, . . . , 2κ), y(x) (n×1) and r(x) It is assumed that the inhomogeneous differential equation where G(x, ξ) is the Green function matrix [30] defined by the following properties: 1. The Green function has the following structure: In addition, it is 2κ times differentiable with respect to x and the derivatives are also continuous functions of x and ξ in the triangles a ≤ should be continuous for x = ξ : In contrast to this, G 2I (x, ξ) and its derivatives are all continuous functions for any are all continuous functions for any x in [a, c].
Though the function G 2I I (x, ξ) and its derivatives are also continuous for x = ξ : . . , n) be an arbitrary constant matrix where the elements are finite. For a fixed ξ ∈ [a, c], the productŷ(x) = G(x, ξ)α as a function of x (x = ξ ) should satisfy the homogeneous differential equation 6. The productŷ = G(x, ξ)α as a function of x should satisfy both the boundary conditions and the continuity conditions: These criteria should be applied to the function pairs G 1I (x, ξ), G 2I (x, ξ) and G 1I I (x, ξ), G 2I I (x, ξ) as well. Remark 1 This definition is a generalization of the definition given for three-point boundary value problems governed by ordinary differential equation in [29,31]. Remark 2 It can be proved by repeating the line of thought presented in Remark 10.4 in book [30] that vector (3.2) given in terms of the previous Green function matrix satisfies the inhomogeneous differential equation system (3.1a) and the boundary and continuity conditions (2.20).

Structure of the general solution
The definition of the Green function matrix given in the previous Section is a constructive one which means that it provides the means that are needed to calculate it. The general solution of the homogeneous differential equation system has the following form , where each column of the matrices Y satisfies the homogeneous differential equation system (4.1), C is a constant square matrix while e is a constant column matrix (a constant vector).
Recalling the fifth property of the definition and the structure of the general solution (4.1), it follows that the elements of the Green function matrix can be given in the same form the expression in the square brackets has in (4.2) -see [30].

Calculation of the Green matrix function if ξ ∈ [a, b]
The matrix of the integration constants C should be different in the two triangular domains a ≤ and where the coefficient matrices A I (ξ ), B I (ξ ) and C I (ξ ) are unknown function matrices. Representation of G 1I (x, ξ) and G 2I (x, ξ) in the forms (4.3) and (4.4) automatically ensures the fulfillment of the first and fifth properties of the definition.
Utilizing the second and third properties of the definition, the following equations can be established for calculating the elements of the matrices B I (ξ ): where 0 is an n × n zero matrix. Let us denote the νth column (ν = 1, 2, . . . , n) of the matrix B I by B I ν : ]. (4.6) The matrix is that of the unknowns for a given ν. Let the ith column (i = 1, 2, . . . , n) in the matrix Y be denoted by η i : ]. (4.8) Then, the matrix B ν is the solution of the linear equation system: while P ν is the transpose of the νth row in the matrix: Note that the coefficient matrix W is the same for each B ν .

Remark 3
The determinant of W is the corresponding Wronskian which is different from zero [32]. This ensures that equation system (4.9a) is solvable.
After having determined the matrices B I , the next step is the calculation of the matrices A I and C I by utilizing the fourth property of the definition. Let α be the νth unit vector in the n × n space: (4.10) With α T , the sixth property of the definition yields The matrices contain the unknowns for a fixed ν in α-note that here the same notational convention is used as in Eqs. (4.6) and (4.7). Making use of Eq. (4.12) we can rewrite equation system (4.11) into the following form: where U ar , U br and U cr are matrices with size n × (2κ × n). The linear equation system (4.13) is solvable for A ν and C ν if its determinant is different from zero. If the boundary and continuity conditions (2.2) are linearly independent equation system (4.13) is, in general, solvable.

Remark 4
Note that the coefficient matrix in (4.13) is the same for each A ν and C ν .

Calculation of the Green function if ξ ∈ [b, c]
The procedure is similar to the one described in the previous Subsection. Thus, it is assumed that and where the coefficients matrices A I I (ξ ), B I I (ξ ) and C I I (ξ ) are again unknown function matrices. The fourth property yields the following equations for calculating the elements of the matrices B (ξ ): (4.17a) (4.17c) By introducing the notation conventions in the same way as for Eq. (4.12), we can manipulate equation system (4.17) into the following form: If the boundary and continuity conditions (2.2) are linearly independent, the above linear equation system is also solvable for A ν and C ν . The proof of this statement is the same as for two-point boundary value problems. The latter is presented in Subsection 9.2.8 of [30].

Governing equations of Timoshenko beams with cross-sectional heterogeneity
A pinned-pinned heterogeneous Timoshenko beam of length L is considered with an intermediate roller support atb as shown in Fig. 1. It has a uniform cross section throughout its length. The centerline of the beam coincides with the axisx of the coordinate system (xŷẑ). Its origin is at the left end of the centerline. The coordinate planexẑ is a symmetry plane. The modulus of elasticity E, the shear modulus of elasticity G and the Poisson number satisfy the relations E(ŷ,ẑ) = E(−ŷ,ẑ), G(ŷ,ẑ) = G(−ŷ,ẑ) and ν(ŷ,ẑ) = ν(−ŷ,ẑ) over the cross section A, which means that they are even functions ofŷ and are independent of the coordinatê x. In this case, we speak about cross-sectional heterogeneity. Note that the E-weighted first moment is zero in this coordinate system: Let A g , I ey be the G-weighted area and E-weighted moment of inertia: The shear correction factor for Timoshenko beams with cross-sectional heterogeneity is denoted by κ y , and its detailed definition is presented in "Appendix A". Equilibrium problems of Timoshenko beams with cross-sectional heterogeneity, subjected to an axial force are governed by the ordinary differential equation system [30,33]: These equations are associated with the boundary and continuity conditions presented in Table 1.
If the beam vibrates freely, the amplitudes will be denoted in the same way as in Eq. (5.5). It can be checked with ease that the amplitude functions should satisfy the differential equations 2 ( ) = 0 Continuity conditions where ρ a is the average density over the cross section and is the ρ weighted moment of inertia. Equation (5.7) in matrix form are Note that the matrices M 0 in (5.6) and (5.9) are different from each other. Consider the differential equation system associated with the boundary and continuity conditions given in Table 1.

Remark 8 The operator L[y]
coincides with the left side of Eqs. (5.6) and (5.9) if the dimensionless axial force N is set to zero. Consequently, the three-point boundary value problem defined by Eq. (5.10) associated with the boundary and continuity conditions in Table 1 describes the mechanical behavior of a pinned-pinned Timoshenko beam with intermediate roller support if the beam is subjected to the dimensionless distributed load r.
Calculation of the Green function for the three point boundary value problem defined by Eq. (5.10) associated with the boundary and continuity conditions presented in Table 1 is detailed in "Appendix B".

Remark 9
With the Green function matrix, the integral is the dimensionless solution of the boundary value problem defined by differential equation system (5.10) and the boundary and continuity conditions of Table 1. Here, −G 11 (x, ξ) and −G 21 (x, ξ) are the dimensionless deflection and rotation at the point x due to a dimensionless unit force r 1 (ξ ) = 1 exerted on the beam at the point ξ , while −G 12 (x, ξ) and −G 22 (x, ξ) are the dimensionless deflection and rotation at the point x due to a dimensionless unit couple r 2 (ξ ) = 1 exerted on the beam at the point ξ . Figures 2, 3, 4 and 5 represent the dimensionless displacement and angle of rotation for the cross section shown in Fig. 13 provided that the corresponding data are those we used to calculate χ by Eq. (A.6). Note that the derivative −dG 22 (x, ξ = 0.75)/dx has a finite jump since the bending moment is discontinuous at this point.

Free vibration of the pinned-pinned beam with intermediate roller support
If the axial force is zero, the amplitudes of the free vibrations should satisfy the differential equitation system associated with the boundary and continuity conditions given in Table 1. The expression on the right side of (5.9) corresponds to r in solution (5.11). Hence, Introduction of the new variables Y and K as shown by the equation dξ results in an eigenvalue problem governed by the Fredholm integral equation system Remark 10 Integral equation system (5.13) coincides formally with integral equation system (9.151) in book [30]-the kernels are, however, different. The later was used to find the eigenvalues for four two-point eigenvalue problems established for pinned-pinned, fixed-pinned fixed-fixed and fixed-free Timoshenko beams. With the eigenvalues, the first two natural circular frequencies were also computed for these beams. The results are shown in Figures 9.3 to 9.6 in [30].
Eigenvalue problem (5.13) can be reduced to an algebraic eigenvalue problem using the boundary element technique. The details are presented in Subsection 9.2.13 of [30]. A FORTRAN 90 code has been developed for solving the algebraic eigenvalue problem. The eigenvalues λ are computed using IMSL Subroutine DGVLRG. With λ k , the following equation can be utilized for calculating the corresponding circular frequency: (5.14) For the beam with cross section shown in Fig. 13  Smaller r means more slender beam. In this case, there cannot be significant difference between the Euler-Bernoulli and the Timoshenko beam theories. Table 2 contains the computational results for √ λ k /π 2 as a function of b -for symmetry reasons b ∈ [0, 0.5]. The results presented in Table 2 coincide with four to five digit accuracy with the results published in [31] for pinned-pinned beams with intermediate roller support within the framework of the Euler-Bernoulli beam theory.
The effect of r on the first four eigenvalues is presented in Table 3 for r = 0.03 and r = 0.05. The results obtained for the first two eigenvalues regarding r = 0.001 443 3 (Euler-Bernoulli beam theory), r = 0.03, r = 0.05 r = 0.075, and r = 0.1 as parameters are depicted in Figs. 6 and 7.
The data in Table 4

Conclusions
The paper is devoted to the issue of how the concept of Green function matrices can be generalized for threepoint boundary value problems governed by ordinary differential equation systems. The investigations are based on the definition of the Green function matrix given for two-point boundary value problems in [30]. The definition is a constructive one since it provides the means necessary for calculating the Green function matrix. Utilizing the definition the Green function matrix is determined for pinned-pinned Timoshenko beams with an intermediate roller support (PrsP beams). After that the self-adjoint three-point eigenvalue problem that governs the free vibration of PrsP Timoshenko beams is reduced to an eigenvalue problem governed by a  homogeneous Fredholm integral equation system with the Green function matrix as its kernel. This eigenvalue problem is reduced to an algebraic eigenvalue problem which is solved numerically. The solutions for the eigenvalues are presented in tabular and graphical formats.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Funding Open access funding provided by University of Miskolc.
A Shear correction factor for heterogeneous beams Figure 12 depicts the cross section of a Timoshenko beam with cross-sectional heterogeneity. The E-weighted first moment Q ey of the shaded area denoted by A in Fig. 12 is defined by the integral The shear correction factor is given by the equation [33]: For the cross section shown in Fig. 13 Here, (A.3f) is the shear correction factor. It also holds that .
On the basis of (4.3) and (4.4), we may write and where the coefficients matrices A I (ξ ), B I (ξ ) and C I (ξ ) are the unknowns. For the sake of further calculations, it is worth partitioning the matrices Y , A I , B I and C I : It follows from Eq. (4.5) that the matrices B I should satisfy the following equations: By introducing the new variables According to Eq.
According to Eq. (4.17b), the product G(x, ξ)α should satisfy the continuity conditions at the intermediate As per Eq. (4.17a), the product G(x, ξ)α should satisfy boundary conditions at the right end of the beam. If 2 , A i , B i and C i into (B.8), (B.9) and (B.10), one arrives at the following equation system: are the solutions. The elements of matrices A 1I , A 2I , C 1I and C 2I can be given in terms of ξ as well. Matrix A 1I : On the basis of (4.14) and (4.15), it can be written that and where the coefficients matrices A I I (ξ ), B I I (ξ ) and C I I (ξ ) are the unknowns. Since B I I (ξ ) = B I (ξ ), we turn our attention to the matrices A I I (ξ ) and C I I (ξ ). It is assumed that A I I (ξ ), B I I (ξ ) and C I I (ξ ) are partitioned in the same manner as the matrices A I (ξ ), B I (ξ ) and C I (ξ ). The notation for these matrix blocks and their elements will be the same as before since this will not cause any misunderstanding. According to Eq.
According to Eq. (4.17a), the product G(x, ξ)α should also satisfy boundary conditions at the right end of the beam. If α T = [1|0] (i = 1) or α T = [0|1] (i = 2), we get: The nonzero solutions are, therefore The elements of matrices A 1I I , A 2I I , C 1I I and C 2I I can be given in terms of ξ . Matrix A 1I I : Matrix C 1I I :