Solutions of Sturm-Liouville Problems

: This paper further improves the Lie group method with Magnus expansion proposed in a previous paper by the authors, to solve some types of direct singular Sturm–Liouville problems. Next, a concrete implementation to the inverse Sturm–Liouville problem algorithm proposed by Barcilon (1974) is provided. Furthermore, computational feasibility and applicability of this algorithm to solve inverse Sturm–Liouville problems of higher order (for n = 2, 4) are veriﬁed successfully. It is observed that the method is successful even in the presence of signiﬁcant noise, provided that the assumptions of the algorithm are satisﬁed. In conclusion, this work provides a method that can be adapted successfully for solving a direct (regular/singular) or inverse Sturm–Liouville problem (SLP) of an arbitrary order with arbitrary boundary conditions.

In [39], authors implemented Barcilon [9]'s iterative algorithm to solve inverse SLPs successfully. The present paper extends the results for the case m = 2 by considering the general-order inverse SLP algorithm as described in [36].
Obtain the corresponding boundary conditions for the linear differential equation using the boundary conditions of the m + 1 eigenvalue problems: Assume that (x, λ) and χ(x, λ) satisfy the m boundary conditions common to all the given m + 1 eigenvalue problems.

3.
Solve Equation (4) using above boundary conditions, to obtain the solutions φ n . 4. Solve the adjoint system of equations to the above system, and denote the solutions by η n . 5.
Find the bi-orthogonal set of functions {y n } ∞ 1 to {φ n } ∞ 1 , using the relation y n (x) = N T η n (x). Some remarks on Algorithm 1 follow:  |σ n − σ (k−1) n | < tol (tol: a fixed tolerance) is checked to exit the loop. Line (4) · denotes the integration w.r.t. x from 0 to 1 and e T = [1,0]. y n and φ n are kept fixed at y 0 and φ 0 , which are the solutions by setting p i = 0, ∀i. Solutions y 0 and φ 0 are calculated using Mathematica [41]. Matlab [42] built-in functions: trapz() and griddedInterpolant() with pchip (piecewise cubic Hermite interpolating polynomial) option are used to approximate integrals and to have p as a function, respectively. The latter is required since, p as a function is input to the Magnus method. Again, the eigenvalue sequence σ (k) is calculated using p (k) and the iteration repeats.
You may refer Section 2.6 of [39] for the specific computational details of the Algorithm 1 for the case m = 1 and next section for the case m = 2.
Set: initial guess p (0) , k = 1; 2 Solve: EVPs using p (0) and combine the solutions to {σ and Solve the EVPs using p (k) and set the solutions to {σ

Results for Direct Singular and Inverse SLPs
Here, we present some numerical examples of direct and inverse SLPs of orders 2 and 4. Algorithm 1 was implemented using MATLAB (2014) [42]. The reference solutions of the EVPs are computed using Wolfram Mathematica 11 [41]. Direct problems were solved using the Magnus expansion method [39].
For the numerical calculations, n is the number of subdivisions in the interval [a, b], m is the number of subdivisions in the interval [λ 0 , λ * ]; λ * being the maximum eigenvalue searching, L is the number of multisection steps used to calculate each eigenvalue in the characteristic function and M is the number of inverse algorithm steps.
Error values are stated using max-norm which is defined as ||x|| ∞ := max n (|x n |).
Here, p(x) is symmetric and normalized. In the inverse Algorithm 1 of [39], w  Figure 1 shows the reconstructed and exact potential, and the log absolute errors in the reconstructed potential, respectively. From Figure 1a, it is obvious that the potential is converging towards the exact one. The max-norm of difference between the reconstructed p and the actual p is ≈3.35 × 10 −3 and the max-norm of difference between the reconstructed eigenvalues and the actual ones is ≈1.62 × 10 −3 .
where δ is the noise level and N = 5, m = 20, n = 100, L = 5, M = 6. It is obvious that reconstruction is possible even in the presence of a significant noise.

Discussion
In the first few examples, the Magnus method [39] is further extended to solve some direct singular SLPs on infinite domain, half-axis and for an asymmetric potential. Later, Barcilon's algorithm is extended to solve inverse SLPs in a non-unit domain and for a non-smooth potential. The first two examples, which extend the Magnus method in [39] to infinite interval, illustrates it's applicability to singular SLPs effectively. Furthermore, the present method does not require transforming the singular SLP to a regular one or re-scaling the eigenvalues. It effectively finds the required eigenvalues, just by truncating the interval to a finite one.
Next, the inverse SLP is solved with smooth and non-smooth potential, even in the presence of noise. It is observed that the method is successful even in the presence of significant noise, provided that the assumptions of the algorithm are satisfied.
The last few examples solve the inverse FSLP using Barcilon's algorithm with the initial knowledge of three spectra. A simplified FSLP is solved keeping one of the unknown potential functions zero.

Conclusions
More than 40 years after Barcilon's paper [36], this paper gives a concrete implementation of the inverse SLP algorithm proposed therein, to our knowledge for the first time. Furthermore, computational feasibility and applicability of this algorithm for solving inverse SLPs of higher order is verified successfully in this paper (for n = 2 and n = 4).
In future work, it is possible to extend the results to solve even higher order SLPs.
In conclusion, this work provides a method that can be adapted successfully for solving a direct (regular/singular) or inverse SLP of an arbitrary order with arbitrary BCs.