Numerical Analysis for Sturm–Liouville Problems with Nonlocal Generalized Boundary Conditions

: For the generalized Sturm–Liouville problem (GSLP), a new formulation is undertaken to reduce the number of unknowns from two to one in the target equation for the determination of eigenvalue. The eigenparameter-dependent shape functions are derived for using in a variable transformation, such that the GSLP becomes an initial value problem for a new variable. For the uniqueness of eigenfunction an extra condition is imposed, which renders the right-end value of the new variable available; a derived implicit nonlinear equation is solved by an iterative method without using the differential; we can achieve highly precise eigenvalues. For the nonlocal Sturm–Liouville problem (NSLP), we consider two types of integral boundary conditions on the right end. For the first type of NSLP we can prove sufficient conditions for the positiveness of the eigenvalue. Negative eigenvalues and multiple solutions may exist for the second type of NSLP. We propose a boundary shape function method, a two-dimensional fixed-quasi-Newton method and a combination of them to solve the NSLP with fast convergence and high accuracy. From the aspect of numerical analysis the initial value problem of ordinary differential equations and scalar nonlinear equations are more easily treated than the original GSLP and NSLP, which is the main novelty of the paper to provide the mathematically equivalent and simpler mediums to determine the eigenvalues and eigenfunctions.

As a special case of Equations ( 1)-(3), Aliyev and Kerimov [1] studied the following spectral problem: w ′ (0) sin β = w(0) cos β, 0 ≤ β < π, (5) w ′ (1) = (aλ 2 + bλ + c)w( 1), where a ̸ = 0, b, c and β are real constants.In terms of a target function f (λ) = w ′ (1; λ) − (aλ 2 + bλ + c)w(1; λ), (7) and root functions system, the authors proved that for the basisness in L 2 , the part of the root space is quadratically close to systems of sines and cosines.In general, the target function is different from the characteristic function, which is usually not attainable in the closed form.The above work was extended in [10] to a boundary condition linearly dependent on the eigenparameter.In this paper, we also set the target equation f (λ) = a 2 (λ)w(1; λ) + b 2 (λ)w ′ (1; λ) = 0 as a mathematical tool to determine the eigenvalues.However, we emphasize the numerical aspect of the spectral problem of Equations ( 1)-( 3), which is not addressed in [1].Currently, the Sturm-Liouville problems with eigenparameter-dependent boundary conditions and transmission conditions are also important issues in many studies [11,12].Recently, many engineering issues and mathematical schemes have been presented resolving the GSLPs, e.g., an analytical method to acquire the sharp estimates for the lowest positive periodic eigenvalue and all Dirichlet eigenvalues of a GSLP [13]; a class of generalized discontinuous GSLPs with boundary conditions rationally dependent on the eigenparameter was solved by using operator theoretic formulation under the new inner product [14]; the GSLP had infinite eigenvalues and corresponding eigenfunctions existed such that the sequence of eigenvalues was increasing as shown in [15]; a thorough formulation of the weighted residual collocation method was based on the Bernstein polynomials for a class of Sturm-Liouville boundary value problems [16]; the infinitely dimensional minimization problem of the lowest positive Neumann eigenvalue for the SLP [17], and an optimal design of a structure was described by a GSLP with a spectral parameter in the boundary conditions [18].
The integral type nonlocal boundary conditions (BCs) are an area of the fast-developing differential equations theory, which may rise up when the solution on the boundary is connected with the values inside the domain.There are some works on the forced Duffing equation with integral boundary conditions [19][20][21][22], which are effective methods for solving the boundary value problem (BVP) with linear integral BCs [22][23][24].
As defined in [25], the boundary shape function (BSF) automatically satisfies the BCs, which includes the solution of BVP as a special member since the solution must exactly satisfy the specified BCs.The boundary shape function method (BSFM) has been adopted to solve some BVPs [25] with conventional local BCs.Recently, Liu [26] provided a unified BSFM to resolve the eigenvalue problems of SLP, GSLP and periodic SLP.However, the method based on the BSF is not yet developed for solving the NSLP endowed with nonlocal BCs.The SLP with nonlocal BCs has been investigated in [27][28][29][30][31], and has been surveyed in [32].Those papers were mainly focused on theoretical analyses, and the numerical algorithm to solve the NSLP is not addressed.
Sequentially, we construct two λ-dependent shape functions and a variable transformation for the GSLP in Section 2. In Section 3, an initial value problem for a new variable with the aid of the normalization technique is derived.To match the right-end BC an implicit nonlinear equation is derived to obtain the eigenvalues by an iterative scheme.
Five testing examples are revealed in Section 4. In Section 5, we discuss the first type NSLP to prove λ > 0, and develop the boundary shape function method (BSFM) for a new variable to iteratively solve the eigenvalue, wherein two examples are given.In Section 6, the second type of NSLP is studied and solved iteratively by the BSFM for a new variable.Two examples are given to display the multiple solutions.To obtain the eigenvalue precisely, we develop a two-dimensional fixed-quasi-Newton method (FQNM) in Section 7 for solving the first type NSLP in the original variable, and the second type NSLP in the new variable.The conclusions are drawn in Section 8.

Boundary Shape Function
We explore a novel method in [33] for solving a nonlinear equation, which is combined to the boundary shape function method [25] to solve Equations ( 1)- (3).For this purpose we need to develop a new mathematical method for transforming Equations ( 1)-( 3) to an initial value problem with a new idea of the eigenparameter-dependent boundary shape function.
Let s 1 (x, λ) and s 2 (x, λ) be shape functions of x and λ: If ) and s 2 (x, λ) satisfy Equations ( 8) and ( 9), then the function w(x), given by w automatically satisfies the boundary conditions ( 2) and ( 3), where T(x, λ) is given by Proof.Applying L λ to Equation (11), and using the first ones in Equations ( 8) and ( 9), we have was taken into account.Hence, we proved Equation (2).Similarly, applying R λ to Equation (11) and using the second ones in Equations ( 8) and ( 9), leads to was taken into account.Thus, Equation (3) was proven.
An eigenparameter-dependent boundary shape function is such a function that it satisfies the eigenparameter-dependent boundary conditions in Equations ( 2) and (3), automatically.Theorem 1 is crucial that the eigenparameter-dependent boundary shape function can be constructed from Equations (11) and (12) easily, with the help of two shape functions s 1 (x, λ) and s 2 (x, λ) and a new function v(x).

A Definite Initial Value Problem
One drawback in Equation (11) is that the right-end values of v(1) and v ′ (1) in a 2 (λ)v(1) + b 2 (λ)v ′ (1) appeared in the translation function T(x, λ) are themselves unknown values.In this situation the variable transformation in Equation ( 11) is hard to work in the present form, since there are three unknown values λ, v(1) and v ′ (1) in T(x, λ) to be determined; v(0) and v ′ (0) in T(x, λ) are some given initial values.We can improve this drawback by considering one extra condition in the following two normalization conditions: where A 0 and B 0 are the given nonzero constants.Hence, w(x) is unique.If the condition w(0) = A 0 is given, we need to specify the value of B 0 , say B 0 = 1; in contrast, if the condition w ′ (0) = B 0 is given, we need to specify the value of A 0 , say 2), we can adopt w(0) = A 0 or w ′ (0) = B 0 as a normalization condition.In general, we take Theorem 2. If a 1 (λ) = 0 and Equation ( 13) is adopted, then the translation function T(x, λ) in Equation ( 12) reads as If a 1 (λ) ̸ = 0 and Equation ( 14) is imposed, then T(x, λ) is given by Proof.Inserting into Equation ( 12), renders If a 1 (λ) = 0, the analysis can be succeeded as follows.Inserting Equation (18) into Equation (11), taking x = 0 and using w(0) = A 0 , we have by Equation (10), 15) is obtained by Equations ( 19) and ( 18) with a 1 (λ) = 0.
Theorem 2 is interesting that by using the normalization condition the translation function T(x, λ) can be reduced to a function that merely involves λ as an unknown value.The reduction from three unknown values λ, v(1) and v ′ (1) to one unknown value λ is a key point for the success of the iterative algorithm to be presented in Section 3.2.Theorem 3. If a 1 (λ) = 0 and Equation ( 13) is adopted, then Equations ( 1)-( 3) can be transformed to where the new variable v(x) is subjected to If a 1 (λ) ̸ = 0 and Equation ( 14) is imposed, then Equations ( 1)-( 3) can be transformed to which is subjected to For Equation (21) or Equation (23), v(0) = c 1 and v ′ (0) = c 2 are the specified initial conditions with c 1 and c 2 given constants.
Although the proof of Theorem 3 is straightforward, it is very important that we can fully transform Equations ( 1)-( 3) to an initial value problem of the second-order ordinary differential equation (ODE) for the new variable v(x), and a nonlinear equation as the target equation to determine the eigenvalue λ.

An Iterative Algorithm
Liu et al. [33] proposed an iterative scheme (LHL) to find the root of f (x) = 0: where with f (r) = 0 and f ′ (r) ̸ = 0. Equation ( 25) has third-order convergence.Choose x 0 and x 2 with r ∈ (x 0 , x 2 ) to render f (x 0 ) f (x 2 ) < 0.Then, we take x 1 = (x 0 + x 2 )/2 and approximate a and b in Equation ( 26) by using the finite differences: The procedures for solving f (x) = 0 are given as follows: (i) Take x 0 and x 2 such that f (x 0 ) f (x 2 ) < 0. (ii) Obtain a and b by Equation ( 27).(iii) For k = 0, 1, . .., iterate The idea to solve Equations ( 1)-( 3) is now simplifying to the shooting method to solve Equations ( 21) and ( 22) or Equations ( 23) and (24).For each given value of λ and the given initial values v(0) = c 1 and v ′ (0) = c 2 , we can integrate Equation ( 21) by employing the fourth-order Runge-Kutta method (RK4).At the endpoint x = 1 of the integration, we can obtain v(1) and v ′ (1) which are implicit functions of λ.Inserting v(1) and v ′ (1) into the target Equation ( 22), an implicit nonlinear equation denoted by f a (λ) is obtained as follows: Similarly, for each given value of λ, we can integrate Equation ( 23) with the given initial values v(0) = c 1 and v ′ (0) = c 2 to obtain v(1) and v ′ (1); hence, Equation ( 24) yields the following implicit nonlinear target equation: Equation (29) or Equation ( 30) is an implicit nonlinear equation of λ, which can be solved by the LHL in Equation (28).
No matter the cases in Equations ( 21) and ( 22), or in Equations ( 23) and ( 24), we unify them to where The iterative processes for solving Equations ( 1)-( 3) to determine the eigenvalue are given as follows: 31) by RK4 to obtain v(1; λ 0 ) and v ′ (1; λ 0 ), and obtain f (λ 0 ); compute Equation ( 31) by RK4 to obtain v(1; λ 2 ) and v ′ (1; λ 2 ), and obtain f (λ 2 ); compute Equation ( 31) by RK4 to obtain v(1; λ 1 ) and v ′ (1; λ 1 ), and obtain f (λ 1 ); compute Because we have mathematically transformed the spectral problem of GSLP to the problem for solving a target equation, the convergence is guaranteed by the LHL with three-order convergence.It is known that the RK4 has fourth-order accuracy which together with the LHL, can achieve highly accurate eigenvalue and the corresponding eigenfunction.

Example 1
Consider a generalized Sturm-Liouville problem in [1]: For this problem we solve λ from the following characteristic equation in [1]: which is the explicit form of the characteristic equation and we solve it by using the Newton method.Here w(x) = cos √ λx is the eigenfunction.As shown in [1], the target equation is which can be arranged to Equation (35).For this spectral problem, the target equation is equal to the characteristic equation.We set A 0 = 1, N = 5000, c 1 = c 2 = 0 and determine λ by Equation (29).Table 1 reveals fast convergence with high accuracy for the eigenvalues obtained.NI denotes the number of iterations, the absolute error AE(λ) of eigenvalue, ERBC the absolute error of the right boundary condition, and ME the maximum error of eigenfunction.
The characteristic equation is derived in [3]: No closed-form eigenfunction exists for this spectral problem.We take B 0 = 1, N = 10,000, c 1 = 0 and c 2 = 1, and iteratively determine λ by Equation (30) as shown in Table 2.

Example 3
Taken from [3], we consider the following generalized Sturm-Liouville problem: For this spectral problem we have It follows from Equation ( 10) that Then, in view of Equations ( 17) and ( 18), one has From Equations ( 23) and ( 24) the second-order ordinary differential equation and nonlinear equation are given by where Now, the iterative method listed Section 3.2 can be applied to compute the eigenvalues of Equations ( 39) and (40).

Example 5
We compute the eigenvalues of the following generalized Sturm-Liouville problem: whose eigenvalues are the same as those in Example 2.
We set A 0 = 1, N = 10,000, c 1 = 0 and c 2 = 1.5 and determine λ by Equation (29).As shown in Table 5, the accuracy in the right boundary condition is slightly better than that in Table 2.

First Type Nonlocal Sturm-Liouville Problem
In this section we consider a nonlocal Sturm-Liouville problem (NSLP) [34,35]: where q(x) is a nonlocal potential.To distinguish the NSLP from the GSLP in Equations ( 1)-( 3), we adopt a different notation y(x) for the eigenfunction, instead of w(x).
In terms of a nonlocal linear boundary operator: Equation ( 54) is written as We can impose a normalization condition: where a 0 ̸ = 0 is a given constant.
Equation ( 61) is somewhat similar to the Rayleigh quotient; however, the Rayleigh quotient in the classical Sturm-Liouville problem includes a potential function q(x) in the integral term 1 0 q(x)y 2 (x)dx at the numerator.Equation ( 61) is important for developing a highly precise iterative algorithm in Section 7, where and the second one in Equation (54) constitute coupled nonlinear equations to determine λ.

Boundary Shape Function
Let p 1 (x) and p 2 (x) be We can derive where we suppose that 1 − 1 0 xq(x)dx ̸ = 0.
Theorem 5 is critical such that the nonlocal boundary conditions in Equation (54) can be automatically satisfied for y(x) with the help of two nonlocal shape functions p 1 (x) and p 2 (x) and a new function u(x).
In Sections 6.3 and 6.4 two examples will be given to demonstrate that the solution y(x) is not unique, even the conditions y(0) = 0, y ′ (0) = a 0 and y ′ (1) + 1 0 q(x)y(x)dx = 0 are imposed.This is different from the local Sturm-Liouville problem; the uniqueness of y(x) is guaranteed if y(0) = 0 and y ′ (0) = a 0 are imposed.

Infinite Solutions
We find that for the second NSLP, there may be many solutions with different ended values of the eigenfunctions.Hence, for the uniqueness of the eigenfunction of Equations ( 53) and ( 78), in addition to Equation (57), we can impose an extra condition: where c 0 are some constants.
To resolve this problem by using the BSFM, we can set a suitable initial value of c 1 = u(0), such that the condition (84) can be satisfied.In Equation (69), we insert x = 1 and use Equations ( 71) and (84) to obtain where u(0) = c 1 and u ′ (0) = c 2 were used.Then, we can compute c 1 by For solving Equations ( 53) and (78) by the BSFM to find a particular solution with an ended value y(1) = c 0 , we can modify the iterative algorithm in Section 6.1 to (i) Give q(x) and derive p 1 (x) and p 2 (x); give c 2 , a 0 , λ 0 , α 0 , β 0 , ϵ and N. (ii) For k = 0, 1, 2, . .., calculate integrate Equation ( 73) by RK4 and take We terminate the iterations if (90)

Remark 1.
There are two unconventional phenomena that happened for the second type of the nonlocal Sturm-Liouville problem.First, the eigenvalue may be negative depending on the end value of y(1)y ′ (1).Second, although the conditions y(0) = 0 and y ′ (0) = a 0 are imposed, the eigenfunction corresponding to the same eigenvalue is not unique, due to y(1) appeared in Equation (53).

Precise Eigenvalue for the Nonlocal Sturm-Liouville Problem
As shown in example 6 in Section 5.4, the number of iterations (NI) is 77, and in example 7 in Section 5.5, the relative error of eigenvalue is 4.36 × 10 −3 .Basically, the slow convergence and low accuracy originate from considering one target equation in Equation (54).When Equation ( 61) is also deemed another target equation, we can improve the accuracy and reduce the number of iterations.Because we need to treat two nonlinear equations, the derivative-free iterative algorithm in Section 3.2 is not applicable.
To improve NI and accuracy, we note the derivative-free iterative algorithm in Section 3.2 by taking b = 0 as a one-dimensional fixed-quasi-Newton method, and extend it to a twodimensional fixed-quasi-Newton method (FQNM).It is known that the two-dimensional Newton method (including the quasi-one) is quadratically convergent.If we can enhance the convergence criterion ϵ to a small value, we can achieve a highly accurate eigenvalue within a reasonable value of NI because both Equations ( 54) and (61) are accurately satisfied.

A Fixed-Quasi-Newton Method
For the first type nonlocal Sturm-Liouville problem in Section 5, it follows from Equations ( 54) and (61) and two nonlinear algebraic equations: where A := y(1) and B := λ.Notice that y(1) = A is appeared in the governing Equation (53), which makes y(x) an implicit function of A; hence, f 1 and f 2 are also implicit functions of A.
As an extension of the one-dimensional FQNM in Section 3 with b = 0, we now consider the approximation of a 2 × 2 Jacobian matrix at the root (A * , B * ) of f 1 (A, B) = f 2 (A, B) = 0 by finite differences, which is termed a fixed-quasi-Newton method (FQNM) of two-dimensions.Choose a rectangle with the vertexes (A * 0 , /2] to carry out the approximation of the Jacobian matrix by a constant matrix {a ij }. The iterative algorithm obtained from the two-dimensional FQNM for solving y(x) in Equations ( 53) and ( 54) is summarized as follows.(i) Give A * 0 , A * 2 , B * 0 , B * 2 , ϵ, and ∆x = 1/N.(ii) Compute , (iii) Let A 0 = A * 0 and B 0 = B * 0 and for k = 0, 1, . . ., iterate Let y 1,k (x) = y (k) (x) and y 2,k (x) = y ′ (k) (x).In each iteration, the RK4 is used to integrate .

Conclusions
Because two unknowns with a missing left-end value of the eigenfunction and eigenvalue are involved, the original generalized Sturm-Liouville problem is more difficult to solve than the presented formulation in the paper.We transformed the GSLP to a definite initial value problem.For the uniqueness of the eigenfunction, by giving a normalization condition we reduced two unknowns to an unknown for the eigenvalue, which is then solved from a target equation in terms of the right boundary condition.The proposed method converged very fast and example tests revealed the high precision of the eigenvalues obtained.We have resolved two types of NSLPs with integral boundary conditions specified on the right end.The appearance of nonlocal potential in the ODE led to some interesting phenomena, which are totally different from the SLP.The methods based on the boundary shape function methods and fixed-quasi-Newton methods were developed, which can accurately determine the eigenvalues and compute the eigenfunctions.For the first type of NSLP, an extra slope condition of the solution on the left end guarantees the uniqueness of the solution and the positiveness of the eigenvalue.For the second type of NSLP, two extra conditions of the slope of the solution on the left end and the specification of the ended value on the right end guarantees the uniqueness of the solution.For the second type of NSLP, the eigenvalue may be negative, and many different solutions happened for different end values.Four numerical testing examples demonstrated that the proposed iterative algorithm could achieve high precision eigenvalues.There are two factors to cause the high accuracy of eigenvalue and eigenfunction: the integration of the resulting ODEs by the fourth-order Runge-Kutta method and the boundary conditions being strictly satisfied.
The novelties involved in the paper were as follows: • Mathematically transforming the generalized Sturm-Liouville problem to an initial value problem of a second-order ODE and an implicit nonlinear equation.

•
Obtaining a very precise eigenvalue of the generalized Sturm-Liouville problem by using LHL on the implicit nonlinear equation.

•
Mathematically transforming nonlocal Sturm-Liouville problems to the initial value problems and an implicit nonlinear equation.

•
For nonlocal Sturm-Liouville problems the initial value problems and two implicit nonlinear equations were derived.

•
The two-dimensional fixed-quasi-Newton method is used to solve two implicit nonlinear equations, quickly obtaining highly accurate eigenvalues.

Fig. 1 .
Fig. 1.For example 6 of the first type nonlocal Sturm-Liouville problem, (a) the residuals and (b) numerical and exact solutions and error for the first eigenfunction; (c) the residuals and (d) numerical solution for the second eigenfunction.

Figure 1 .
Figure 1.For example 6 of the first type nonlocal Sturm-Liouville problem, (a) the residuals and (b) numerical and exact solutions and error for the first eigenfunction; (c) the residuals and (d) numerical solution for the second eigenfunction.

Fig. 2 .
Fig. 2. For example 7 of the first type nonlocal Sturm-Liouville problem, (a) the residuals and (b) numerical and exact solutions and error

Figure 2 .
Figure 2.For example 7 of the first type nonlocal Sturm-Liouville problem, (a) the residuals and (b) numerical and exact solutions and errors.

Fig. 3 .
Fig. 3.For example 8 of the second type nonlocal Sturm-Liouville problem, (a) the residuals and (b) numerical and exact solutions and error.

Figure 3 .
Figure 3.For example 8 of the second type nonlocal Sturm-Liouville problem, (a) the residuals and (b) numerical and exact solutions and error.

Fig. 4 .
Fig. 4. For example 9 of the second type nonlocal Sturm-Liouville problem, (a) the residuals and (b) numerical and exact solutions and error.

Figure 4 .
Figure 4.For example 9 of the second type nonlocal Sturm-Liouville problem, (a) the residuals and (b) numerical and exact solutions and error.

Fig. 5 .
Fig. 5.For example 9 comparing three solutions with different ended values.Figure 5.For example 9 comparing three solutions with different ended values.

Figure 5 .
Fig. 5.For example 9 comparing three solutions with different ended values.Figure 5.For example 9 comparing three solutions with different ended values.

Fig. 7 .
Fig. 7.For example 7 solved by the FQNM, the improvements of (a) the residuals and (b) numerical and exact solutions and error

Figure 7 .
Figure 7.For example 7 solved by the FQNM, the improvements of (a) the residuals and (b) numerical and exact solutions and error.

1 Fig. 8 .
Fig. 8.For example 8 comparing four solutions with different ended values.The first solution is obtained by the BSFM, while other solutions are obtained by the BSFM and FQNM.

Figure 8 .
Figure 8.For example 8 comparing four solutions with different ended values.The first solution is obtained by the BSFM, while other solutions are obtained by the BSFM and FQNM.