Electronic Journal of Differential Equations, Vol. 2024 (2024), No. 53, pp. 1–17. ISSN: 1072-6691. URL: https://ejde.math.txstate.edu, https://ejde.math.unt.edu DOI: 10.58997/ejde.2024.53 DEEP LEARNING METHOD FOR FINDING EIGENPAIRS IN STURM-LIOUVILLE EIGENVALUE PROBLEMS SEN ZHANG, JIAN ZU, JINGQI ZHANG Abstract. Solving the eigenvalue problem for differential equations in inho- mogeneous media poses a significant challenge across diverse scientific fields. While classical finite difference methods and finite element methods have pro- duced numerous outcomes, they heavily rely on discretizing the computational domain, which can introduce complexities and limitations. In this study, we present an unsupervised neural network approach tailored for finding eigen- pairs in Sturm-Liouville eigenvalue problems within inhomogeneous media. Our method introduces eigenvalues as trainable parameters, crafts a novel cost function, incorporates an adaptive hyper-parameter tuning strategy, and se- quentially trains the eigenpairs. The simplicity, accuracy, and interpretability of our approach significantly expand its applicability across various domains. The method we present in this paper can easily tackle boundary value condi- tions with derivatives, resulting in orthogonal eigenfunctions. This is a very important advantage of deep learning methods that has not yet been noticed. Quantitative estimation of eigenpairs is given for the Sturm-Liouville eigen- value problems. Furthermore, we extend the proposed methodology to tackle two-dimensional cases, periodic scenarios, demonstrating its versatility and broad potential. 1. Introduction The Sturm-Liouville eigenvalue problem holds substantial significance across a diverse range of fields. It can be demonstrated that there are an infinite number of eigenvalues, each with a unique eigenfunction, and that these eigenfunctions form an orthonormal basis within a specific Hilbert space of functions. The adaptability of the Sturm-Liouville equation to various boundary value conditions allows for the precise representation of specific physical phenomena. For example, in the con- text of heat equations, Dirichlet conditions maintain a consistent temperature at the rod’s extremities, while Neumann conditions simulate insulating environments with no heat transfer and negligible temperature gradients at the ends. Further- more, more comprehensive boundary value conditions, such as the Robin condition, accommodate partially insulated boundaries. In practical applications, it is also crucial to compute multiple orthogonal eigen- vectors of the Sturm-Liouville eigenvalue problem. For instance, in the construction 2020 Mathematics Subject Classification. 34L16, 65F18, 65L15, 65N25, 68T07. Key words and phrases. Eigenvalue problems; Sturm-Liouville problems; neural network; self-adjoint operator; inhomogeneous media. ©2024. This work is licensed under a CC BY 4.0 license. Submitted December 2, 2023. Published September 12, 2024. 1 2 S. ZHANG, J. ZU, J. ZHANG EJDE-2024/53 of a suitable Lyapunov function, there are instances where the need for several mu- tually orthogonal minimal eigenvectors arises, see [4]. Therefore, the quest for mul- tiple minimal eigenpairs (i.e., eigenvalues and eigenvectors) for the Sturm-Liouville eigenvalue problem remains an enduring research pursuit. Classical numerical methods employed in the quest for eigenpairs include finite difference methods [1, 17] and finite element methods [2]. The finite difference method, characterized by its simplicity and ease of comprehension, is a popular choice. However, it also exhibits inherent limitations. For instance, in the case of Sturm-Liouville eigenvalue problems with derivative boundary value conditions, such as Neumann boundary and Robin boundary, the finite difference approxima- tion of the differential operator may not preserve self-adjointness. Consequently, the resulting eigenvectors are not orthogonal to each other, which contradicts the nature of Sturm-Liouville differential operators being self-adjoint and having orthogonal eigenvectors. The finite element method offers high computational accuracy and wide applicability. However, it requires manual specification of appropriate finite element units prior to computation, and its effectiveness is often reliant on indi- vidual expertise. It is a more complex method and may not be user-friendly for non-experts. Furthermore, it is worth noting that both finite difference methods and finite element methods yield numerical solutions, which lack analytical expressions and cannot be readily used for operations such as differentiation. Their applicabil- ity is relatively limited in situations requiring analytical expressions. The search for effective analytical approximations remains a significant scientific challenge. In recent years, deep learning methods have rapidly advanced and integrated into diverse aspects of people’s lives. Within the study of differential equations, their exceptional strengths and capabilities stand out. Chen et al. [6] proposed Neural ODE, which approximate the right-hand side function of differential equations us- ing deep neural networks, achieving remarkable results. Currently, the most widely used approach is to approximate the solutions of differential equations using deep neural networks, with advantages like increased applicability, ease of understanding, and the ability to perform differentiable operations. Raissi et al. [19] introduced a deep learning framework known as ’Physics-informed neural networks’ (PINN). Within this innovative framework, the physical equations themselves are seamlessly integrated into the network as constraints, thereby enabling the solution of partial differential equations (PDEs) without the necessity of a dataset. Pang et al. [18] proposed a fractional PINN (fPINN), addressing the challenge of automatic differen- tiation being unsuitable for fractional operators. Yang et al. [20] further advanced the field by developing a novel class of physics-informed generative adversarial net- works, which effectively tackle forward, inverse, and mixed stochastic problems in a unified manner, relying on a limited number of scattered measurements. Lu et al. [16] have presented a Python library, named DeepXDE, specifically designed for PINNs. This library serves dual purposes: As an educational tool for classroom use and as a research tool for addressing computational science and engineering challenges. For specific mathematical problems, such as identifying Hamiltonian systems, Mattheakis et al. [7] developed a Hamiltonian neural network (HNN) with the pri- mary aim of conserving the symplectic structure of the Hamiltonian system. This allows for more accurate learning and prediction of the dynamic behavior of the EJDE-2024/53 DEEP LEARNING METHOD FOR FINDING EIGENPAIRS 3 system. Jin et al. [13] have developed symplectic networks (SympNets) to effec- tively learn the symplectic flow within Hamiltonian systems. Moreover, they [12] have also proposed Poisson neural networks (PNNs) which are specifically designed to learn Poisson systems and accurately capture the trajectories of autonomous sys- tems, all derived directly from data. Lu et al. [15] proposed DeepONet, which has demonstrated its efficacy as a powerful tool for handling nonlinear operators using supervised data-driven methods. What’s even more thrilling is the potential that emerges from merging DeepONet with the physics encoded by PINNs. This union opens up the possibility of achieving precise, real-time predictions in many fields with extrapolation capabilities. Zang et al. [21] presented an adversarial neural net- work rooted in partial differential equation weak solutions, offering a pathway for tackling high-dimensional partial differential equations. For further advancements related to physics-informed machine learning, we highly recommend delving into the comprehensive review article [14] for a deeper understanding. The application of deep learning techniques in addressing eigenvalue problems has witnessed significant advancements. PINN, owing to its ease of implementation and ability to maintain the self-adjointness of the operator throughout the compu- tational process, has been employed by scholars to identify orthogonal eigenpairs. Jin et al. [9] used PINN to solve the quantum problems related to finite wells, mul- tiple finite wells, and hydrogen atom eigenvalues, and introduced a strong bound- ary value condition. However, their investigation was limited to one-dimensional homogeneous media under the Dirichlet boundary value condition. Holliday et al. [8] extended Jin’s algorithm to quantum billiard eigenvalue problems, success- fully identifying the eigenvalues and eigenfunctions of the stationary Schrödinger differential equation with Dirichlet boundary value conditions in two-dimensional domains. Ben-Shuaul et al. [5] have effectively employed PINN techniques to tackle the problem of finding the smallest eigenpairs for one-dimensional Sturm-Liouville boundary value problems. The strength of their approach lies in its ability to simul- taneously identify multiple eigenpairs. However, their method, which incorporates the Rayleigh quotient for eigenvalue calculation, can be challenging for non-experts and may yield lower precision results. Additionally, their focus is specifically on scenarios involving the Laplace operator as the differential operator, disregarding in- homogeneous cases. Furthermore, their framework is limited to Dirichlet boundary value condition and does not account for derivative boundary situations. Unlike the Dirichlet boundary, the discrete matrix corresponding to the Sturm-Liouville differ- ential operator under derivative boundaries is not self-adjoint when using classical numerical methods. This aspect underscores the advantages of leveraging neural networks for addressing such problems. In this article, we consider how to find the smallest eigenpairs using deep learn- ing methods. We introduce eigenvalues as trainable parameters, construct a novel cost function, add an adaptive hyper-parameter tuning strategy, and enhance the speed and accuracy of training by sequentially training eigenpairs. This setup al- lows us to compute a more extensive set of smallest orthogonal eigenpairs. By utilizing a weighted inner product space and adaptive hyper-parameters, we tackle the Sturm-Liouville eigenvalue problem under varying boundary value conditions in inhomogeneous media. The simplicity, accuracy, and interpretability of our method 4 S. ZHANG, J. ZU, J. ZHANG EJDE-2024/53 significantly broaden the applicability of the proposed approach. Quantitative esti- mation of eigenpairs is given for the Sturm-Liouville eigenvalue problems. Further- more, we extend the proposed method to two-dimensional cases, periodic scenarios. The organization of this paper is as follows. In Section 2, we introduce the prob- lem model, construct a novel cost function, and present our algorithm. In Section 3, we give a quantitative estimation of eigenpairs. In Section 4, we compare our ap- proach with previous algorithms, including state-of-the-art deep learning methods and classical numerical techniques. Section 5 presents extensive numerical exper- iments that demonstrate the efficiency and accuracy of our method in addressing the eigenvalue problem under various boundary conditions. Finally, in Section 6, we conclude the paper by summarizing our key findings and contributions. 2. Model problem In this article, we aim to find the M smallest eigenvalues and eigenfunctions of the general Sturm-Liouville eigenvalue problems (p(x)φ′(x))′ + λσ(x)φ(x) = 0, (2.1) with homogeneous boundary value conditions a1φ ′(0)− b1φ(0) = 0, a2φ ′(π) + b2φ(π) = 0, a2i + b2i ̸= 0, i = 1, 2, (2.2) where p ∈ C1((0, π),R+) and σ ∈ C((0, π),R+), λ is the eigenvalue and φ ∈ C2((0, π),R) is the corresponding eigenfunction. In some instances, p and σ are elastic coefficent and rock density coefficient respectively, see [3]. By a change of variable x→ z given by z = ∫ x 0 √ σ(s) p(s) ds, equation (2.1) leads to (ρ(z)φ′(z))′ + λρ(z)φ(z) = 0, (2.3) where ρ(z) = √ σ(z)p(z) is the impedance function. The mathematical theoretical results for the Sturm-Liouville eigenvalue problems (2.1) or (2.3) are relatively abundant. Let Ω = (0, π). Define the Sturm-Liouville operator in inhomogeneous media by L[φ] := − (p(x)φ′(x))′ σ(x) = λφ, in Ω. (2.4) Define the weighted inner product by (φ,ψ)σ := ∫ Ω φψ∗σdx, φ, ψ ∈ L2(Ω, σdx), which induces the L2 norm ∥ · ∥σ. Note that when σ ≡ 1 on Ω, ∥ · ∥σ reduces to be the standard L2 norm, which is denoted by ∥ · ∥. It is well known that L is a self- adjoint operator, the eigenvectors of L construct an orthonormal basis (φk)k∈N+ of L2((0, π), σdx): L[φk] = λkφk, k = 1, 2, . . . , (2.5) where (λk)k∈N is increasingly approach to +∞ as k → +∞. Finding analytical expressions for the eigenpairs of Sturm-Liouville problems in inhomogeneous media is highly challenging, so the deep learning solution we de- signed can only be compared with the classical numerical method. Using deep learning methods to solve eigenvalue problems offers numerous advantages. One such advantage is that the resulting neural network solutions are semi-analytical, meaning they are differentiable. The method is also simple, practical, and easily EJDE-2024/53 DEEP LEARNING METHOD FOR FINDING EIGENPAIRS 5 extensible, making it highly accessible to non-experts. Personally, one of the great- est benefits I have found is that for derivative boundary value cases, the neural network solutions obtained through this dimension-free approach are self-adjoint, aligning with the original problem. In contrast, classical numerical methods, such as finite difference method, often yield approximate solutions that do not satisfy self-adjointness, leading to numerical eigenfunctions that are not orthogonal, see [1]. Using neural networks to find the eigenpairs of {λk, φk} in inhomogeneous me- dia only requires starting with the problem itself. We gradually learn eigenpairs, starting from smaller values of k and progressing to larger ones. Inspired by the PINN, we integrate the equation equalities, boundary value conditions, normaliza- tion of eigenfunctions, and orthogonality among distinct eigenfunctions into the loss function. This approach ensures the feasibility and practical implementation of our method. Specifically, we introduce a parameterized form of the k-th eigenfunction, denoted by φθk . In this context, θk functions as the primary training parame- ter. Moreover, within our comprehensive framework, we incorporate the eigenvalue λk as an additional parameter that undergoes training alongside θk. For each k, when φθk learns the expected target through training, λk will also be automati- cally updated with the network to the eigenvalue corresponding to the eigenfunction network φθk . This design greatly reduces the complexity of the network and the difficulty of calculation, so that we only need to pay attention to the training of the eigenfunction network, and do not need to design an additional network to train the eigenvalue. To optimize the neural network φθk with respect to its parameter θk, we define the loss function J(θk, λk) = αJeq(θk, λk) + βJbdry(θk) + γJnor(θk) + δkJorth(θk), (2.6) where α, β, γ, δk are hyper-parameters that control the relative importance of each term. The individual terms are as follows: Jeq(θk, λk) = ∥L[φθk ]− λkφθk∥2σ (2.7) enforces the constraints imposed by the equation. Jbdry(θk) = ∣∣B[φθk ] ∣∣2 (2.8) accounts for the constraints imposed by the boundary value conditions. Jnor(θk) = ( ∥φθk∥2σ − 1 )2 ∥φθk∥2σ , (2.9) is responsible for the normalization of the eigenfunction in the weight space. The denominator is included to prevent the norm of the eigenfunction from vanishing during the computation of higher eigenvalues. Jorth(θk) = k−1∑ j=1 ∣∣(φθk , φθj )σ ∣∣2, (2.10) ensures orthogonality of the current eigenfunction with respect to previously com- puted eigenfunctions in the weight space. Under this framework, the choice of δk is crucial for accurately computing higher- order eigenvalues λk. When the epoch is less than half of the total number of epochs, we set δk equal to 10 · λk−1 to ensure that the eigenvalues can correctly 6 S. ZHANG, J. ZU, J. ZHANG EJDE-2024/53 approach the desired values. However, when the epoch exceeds half of the total, we set δk to 1 to reduce the weight of the orthogonality term, thereby facilitating better computation of the loss associated with the equation term. Note that both the divisor of the third term and the addition of the δk are specially designed for obtaining the higher eigenpairs. For details of the process, see Algorithm 1. Algorithm 1 Deep learning method for finding M eigenpairs 1: for k = 1 to M do 2: Set total number of epochs Nepochs 3: Define hyper-parameters α, β, γ, λ0 4: Initialize neural network φθk with parameters θk and λk = λk−1 5: for epoch = 1 to Nepochs do 6: if epoch < Nepochs 2 then 7: Set δk = 10 · λk−1 8: else 9: Set δk = 1 10: end if 11: Compute loss function J(θk, λk) using Eq. (2.6) 12: Compute gradients of J(θk, λk) with respect to θk and λk 13: Update θk and λk using gradient descent optimizer 14: if convergence criterion is met then 15: Break the loop 16: end if 17: end for 18: Output optimized parameters θk and eigenvalue λk 19: end for In the individual terms of cost function, Jeq(θk, λk) is calculated by ∥L[φθk ]− λkφθk∥2σ = 1 Nk Nk∑ j=1 σ ( x(j)r )( L[φθk ](x (j) r )− λkφθk(x (j) r ) )2 , (2.11) where {x(j)r } are sets of inner points sampled uniformly in the entire domain Ω. The weight norm ∥ · ∥σ and weight inner product in Jnor and Jorth are calculated by the mid-point integral method. Indeed, our framework is applicable not only to 1-dimensional cases but also to 2-dimensional cases, as well as periodic scenarios. For the 1D case, the weight inner product is calculated by the mid-point integral method: (φθk , φθj )σ = h N−1∑ m=0 σ(x (m) h )φθk(x (m) h )φθj (x (m) h ) (2.12) where {x(m) h |x(m) h = h/2 + m · h} are sets of integral points with h = π/N . The boundary value conditions for the eigenfunction φ are defined by B[φ] = [a1φ ′(0)− b1φ(0), a2φ ′(π) + b2φ(π)]. (2.13) In the experiments, the eigenfunction network φθk is set as a fully connected neural network with 4 hidden layers, each consisting of 32 neurons. Silu is used as the activation function. During training, we utilize the Adam optimizer to update the parameters of the neural network, with the initial learning rate of 0.01, EJDE-2024/53 DEEP LEARNING METHOD FOR FINDING EIGENPAIRS 7 which drops by 10% every 1000 epoch. A total of 50, 000 epochs are set, ensuring sufficient time for network convergence. This computation is performed using the backpropagation algorithm implemented in PyTorch, a deep learning framework that facilitates automatic gradient computation. 3. Quantitative estimation of eigenpairs We consider the Sturm-Liouville eigenvalue problem: (ρ(x)φ′(x))′ + λρ(x)φ(x) = 0, (3.1) with homogeneous boundary value conditions a1φ ′(0)− b1φ(0) = 0, a2φ ′(π) + b2φ(π) = 0, a2i + b2i ̸= 0, i = 1, 2. (3.2) The boundary value conditions are divided into 6 cases: (1) Dirichlet boundary value condition: a1 = 0, b1 > 0,a2 = 0, b2 > 0; (2) Neumann boundary value condition: a1 > 0, b1 = 0, a2 > 0, b2 = 0; (3) Dirichlet-Neumann boundary value condition: a1 = 0, b1 > 0, a2 > 0, b2 = 0 or a1 > 0, b1 = 0, a2 = 0, b2 > 0; (4) Mixed boundary value conditions I: a1 = 0, b1 > 0, a2 > 0, b2 ̸= 0 or a1 > 0, b1 ̸= 0, a2 = 0, b2 > 0; (5) Mixed boundary value conditions II: a1 > 0, b1 = 0, a2 > 0, b2 ̸= 0 or a1 > 0, b1 ̸= 0, a2 > 0, b2 = 0; (6) General boundary value condition: a1 > 0, b1 ̸= 0, a2 > 0, b2 ̸= 0. Assume that ρ(x) ∈ H2((0, π),R+). Let z(x) = √ ρ(x)φ(x), then the above equa- tion can be transformed to z′′(x) + (λ− ηρ(x))z(x) = 0, where ηρ(x) = 1 2 ρ′′ ρ − 1 4 ( ρ′ ρ ) 2 is assumed to be a real function in L2(0, π). The boundary value condition is transformed to a new one: a1z ′(0)− (a1 2 ρ′(0) ρ(0) + b1 ) z(0) = 0, a2z ′(π) + ( b2 − a2 2 ρ′(π) ρ(π) ) z(π) = 0. Let (Λk, φk) be the exact eigenpairs of (3.1)-(3.2). Let (λk, φθk) be approximate eigenpairs calculated by our Neural Networks. Denote ρ0 := ess inf ηρ(x) and ρ1 := ρ0 + 2 π ∫ π 0 {ηρ(x)− ρ0}+ dx, where {h(x)}+ = max{h(x), 0}. Now, we give the estimation of Λk for Sturm- Liouville eigenvalue probems. The eigenvalue estimation for Dirichlet boundary value conditions has been well studied. Lemma 3.1 ([22]). The Sturm-Liouville problems (3.1) with Case (1) have the following inequalities: k2 + ρ0 ≤ Λk ≤ k2 + ρ1, k ∈ N+. 8 S. ZHANG, J. ZU, J. ZHANG EJDE-2024/53 The proof of above lemma can be seen in [22]. Notice that we have not assumed that ρ0 > 0, so the case ρ(x) ≡ constant is included. Denote α1 := ρ′(0) 2ρ(0) + b1 a1 and α2 := − ρ′(π) 2ρ(π) + b2 a2 . For other boundary value problems, we have the following result. Lemma 3.2. Let k ∈ N+. If a2 = 0 and α1 ≥ 0 (resp. a1 = 0 and α2 ≥ 0), the Sturm-Liouville problems (3.1) with Case (3) and (4) satisfy( k − 1 2 )2 + ρ0 ≤ Λk ≤ ( k − 1 2 )2 + ρ1 + 2α1 π ( resp. ≤ ( k − 1 2 )2 + ρ1 + 2α2 π ) . Proof. The case a2 = 0 and α1 = 0 can be found in [11]. Now we prove that it still holds for the case a2 = 0 and α1 > 0. It is obviously that Λk > ρ0, k ∈ N+. We introduce the Prüfer transformation: zk = r sin θ, z′k = √ Λk − ρ0r cos θ, where r(x) > 0. Assume that zk(x) has (k − 1) zeros in the interval (0, π) denoted as 0 = κ0 < κ1 < κ2 < · · · < κk−1 < κk = π. Define θ(κi) = iπ, i = 1, 2, . . . , k − 1. Using the boundary conditions, we obtain θ(0) = arctan ( 1 α1 √ Λk − ρ0 ) , θ(π) = kπ. Firstly, we estimate the lower bound of Λk. It is not difficult to know that θ′ = √ Λk − ρ0 − (ηu(x)− ρ0) sin 2 θ√ Λk − ρ0 ≤ √ Λk − ρ0. Integrating both sides of the inequality over (0, π), we obtain θ(π)− θ(0) ≤ √ Λk − ρ0π. Since α1 > 0, utilizing the properties of arctan(x), we know that θ(0) = arctan ( 1 α1 √ Λk − ρ0 ) ∈ (0, π 2 ), i.e., θ(π)− θ(0) > ( k − 1 2 ) π. Thus, we conclude that Λk ≥ ( k − 1 2 )2 + ρ0 Secondly, we prove the upper bound of Λk. Let pk = √ Λk − ρ0. Then, as an increasing function of k, pk approaches infinity and satisfies pk > 1/2. Notice that lim k→+∞ pk ( arctan ( pk α1 ) − π 2 ) = −α1. Define f(p) := p ( arctan ( p α1 ) − π 2 ) . Taking the derivative of f(p), we obtain f ′(p) = arctan ( p α1 ) − π 2 + α1p p2 + α2 1 . and limp→+∞ f ′(p) = 0. Since f ′′(p) = 2α1 p2 + α2 1 − 2α1p 2 (p2 + α2 1) 2 = 2α3 1 (p2 + α2 1) 2 > 0, EJDE-2024/53 DEEP LEARNING METHOD FOR FINDING EIGENPAIRS 9 we obtain f ′(p) < 0, p ∈ [ 12 ,+∞). It implies that√ Λk − ρ0 ( arctan (√Λk − ρ0 α1 ) − π 2 ) ≥ −α1. Hence θ(0) ≥ π 2 − α1√ Λk − ρ0 . By the Prüfer transformation, we obtain θ′ = √ Λk − ρ0 cos 2 θ − (ηu(x)− ρ0) sin 2 θ√ Λk − ρ0 ≥ √ Λk − ρ0 − {ηu(x)− ρ0}+ sin2 θ√ Λk − ρ0 ≥ √ Λk − ρ0 − {ηu(x)− ρ0}+√ Λk − ρ0 . By integrating both sides over the interval (0, π), we obtain: θ(π)− θ(0) ≥ √ Λk − ρ0π − ∫ π 0 {ηu − ρ0}+dx√ Λk − ρ0 . Therefore, kπ − π 2 + α1√ Λk − ρ0 ≥ √ Λk − ρ0π − ∫ π 0 {ηu − ρ0}+dx√ Λk − ρ0 . Defining D := Λk − ρ0, B := k − 1 2 , C := 1 π ( ∫ π 0 {ηu(x) − ρ0}+dx + α1), we obtain D ≤ B √ D + C, which implies √ D ≤ ( √ B2 + 4C +B)/2. By inequalities √ B2 + 4C = B √ 1 + 4C B2 ≤ B(1 + 2C B2 ) = B + 2C B , we obtain the upper bound D ≤ B2 + 2B √ B2 + 4C +B2 + 4C 4 ≤ B2 + 2C. Thus, we conclude that Λk ≤ ( k − 1 2 )2 + ρ1 + 2α1 π . The aforementioned conclusion can be easily extend to the case where a1 = 0 and α2 ≥ 0. We omit the proof for brevity. □ Lemma 3.3. Assume that α1 ≥ 0 and α2 ≥ 0. Then the Sturm-Liouville problems (3.1) with Case (2), (5) and (6) satisfy (k − 1)2 + ρ0 ≤ Λk ≤ (k − 1)2 + ρ1 + 2α1 π + 2α2 π , k ∈ N+. The above lemma can be proved by a similar argument to that of Lemma 3.2. Now, we give a rough estimation of λk obtained by the neural network method for the Sturm-Liouville eigenvalue problems. 10 S. ZHANG, J. ZU, J. ZHANG EJDE-2024/53 Proposition 3.4. Assume that ρ1− ρ0 ≤ 1. If λk ∈ [k2+ ρ0, k 2 + ρ1], k ∈ N+, the Sturm-Liouville eigenvalue problems with Case (1) imply that |λk − Λk| ≤ 1 ≤ min j ̸=k {|λk − Λj |}. Proof. By Lemma 3.1, we have |λk − Λk| ≤ |ρ1 − ρ0| ≤ 1. For j > k ≥ 1, we have |λk − Λj | ≥ |Λj − Λk| − |λk − Λk| ≥ |j2 + ρ0 − k2 − ρ1| − |k2 + ρ1 − k2 − ρ0| ≥ |j2 − k2| − 2|ρ1 − ρ0| ≥ 1. The second inequality follows from Lemma 3.1, and the fourth inequality follows from the fact that j2 − k2 ≥ j + k ≥ 3. Similarly, for k > j ≥ 1, we also conclude that |λ− Λj | ≥ 1. Thus, for j ̸= k ∈ N+, we have |λk − Λk| ≤ |ρ1 − ρ0| ≤ 1 ≤ min j {|λk − Λj |}. □ Proposition 3.5. Assume that a2 = 0, α1 ≥ 0 and ρ1−ρ0+ 2α1 π ≤ 2 3 (resp. a1 = 0, α2 ≥ 0 and ρ1 − ρ0 + 2α2 π ≤ 2 3). If λk ∈ [ ( k − 1 2 )2 + ρ0, ( k − 1 2 )2 + ρ1 + 2α1 π ] (resp. λk ∈ [ ( k − 1 2 )2 + ρ0, ( k − 1 2 )2 + ρ1 + 2α2 π ]), k ∈ N+, the Sturm-Liouville eigenvalue problems with Case (3) and (4) imply that |λk − Λk| ≤ 2 3 ≤ min j ̸=k {|λk − Λj |}. By a similar argument to that of Proposition 3.4, we can prove Proposition 3.5 based on Lemma 3.2. Proposition 3.6. Assume that α1 ≥ 0, α2 ≥ 0 and ρ1 − ρ0 + 2α1 π + 2α2 π ≤ 1 3 . If λk ∈ [(k−1)2+ρ0, (k−1)2+ρ1+ 2α1 π + 2α2 π ], k ∈ N+, the Sturm-Liouville eigenvalue problems with Case (2), (5) and (6) imply that |λk − Λk| ≤ 1 3 ≤ min j ̸=k {|λk − Λj |}. By a similar argument to that of Proposition 3.4, we can prove Proposition 3.6 based on Lemma 3.3. Denote Ek(x) = L[φθk ](x)−λkφθk(x), we can estimate the error of the eigenvalue by ∥Ek(x)∥σ and ∥φθk∥σ. Theorem 3.7 (Eigenvalue estimation). If |λk−Λk| ≤ minj ̸=k{|λk−Λj |}, k ∈ N+, we have |λk − Λk| ≤ ∥Ek(x)∥σ ∥φθk∥σ . Proof. Note that ∥L[φθk ]− λkφθk∥2σ = (L[φθk ]− λkφθk , L[φθk ]− λkφθk)σ = ∞∑ j=1 (Λj − λk) 2|(φθk , φj)σ|2 EJDE-2024/53 DEEP LEARNING METHOD FOR FINDING EIGENPAIRS 11 ≥ (Λk − λk) 2 ∞∑ j=1 |(φθk , φj)σ|2 Thus, we conclude that |λk − Λk| ≤ ∥Ek(x)∥σ ∥φθk∥σ . □ We define Pkφ = (φ,φk)σφk and P⊥ k φ = ∑∞ j ̸=k(φ, φj)σφj . Obviously, φθk = Pkφθk + P⊥ k φθk . Next, we use ∥P⊥ k φθk | to evaluate the eigenfunction of our algo- rithm. Theorem 3.8 (Eigenfunction estimation). Assume that ρ1 − ρ0 < 1. If λk ∈ [k2 + ρ0, k 2 + ρ1], k ∈ N+, the Sturm-Liouville eigenvalue problems with Case (1) satisfy ∥P⊥ 1 φθ1∥σ ≤ ∥E1(x)∥σ |Λ2 − λ1| , or ∥P⊥ k φθk∥σ ≤ ∥Ek(x)∥σ min{|Λk+1 − λk|, |Λk−1 − λk|} , k ≥ 2. Proof. By Proposition 3.4, we have |λk − Λj | ≥ 1, k ̸= j. For k = 1, we obtain ∥L[φθ1 ]− λ1φθ1∥2σ = (L[φθ1 ]− λ1φθ1 , L[φθ1 ]− λ1φθ1)σ = ∞∑ j=1 (Λj − λ1) 2|(φθ1 , φj)σ|2 ≥ (Λ2 − λ1) 2 ∞∑ j=2 |(φθ1 , φj)σ|2 + (Λ1 − λ1) 2|(φθ1 , φ1)σ|2 ≥ (Λ2 − λ1) 2 ∞∑ j=2 |(φθ1 , φj)σ|2 Thus, we conclude that ∥P⊥ 1 φθ1∥σ ≤ ∥E1(x)∥σ |Λ2 − λ1| For k ≥ 2, we obtain ∥L[φθk ]− λkφθk∥2σ = (L[φθk ]− λkφθk , L[φθk ]− λkφθk)σ = ∞∑ j=1 (Λj − λk) 2|(φθk , φj)σ|2 ≥ min{(Λk+1 − λk) 2, (Λk−1 − λk) 2} ∞∑ j=1,j ̸=k |(φθk , φj)σ|2 Thus, we conclude that ∥P⊥ k φθk∥σ ≤ ∥Ek(x)∥σ min{|Λk+1 − λk|, |Λk−1 − λk|} □ Remark 3.9. By Proposition3.5 and 3.6, Theorem 3.8 can be easily extended to: 12 S. ZHANG, J. ZU, J. ZHANG EJDE-2024/53 (1) Sturm-Liouville eigenvalue problems with Case (3) and (4), if a2 = 0, α1 ≥ 0, ρ1−ρ0+ 2α1 π ≤ 2 3 and λk ∈ [ ( k − 1 2 )2 +ρ0, ( k − 1 2 )2 +ρ1+ 2α1 π ] (resp. a1 = 0, α2 ≥ 0, ρ1−ρ0+ 2α2 π ≤ 2 3 and λk ∈ [ ( k − 1 2 )2 +ρ0, ( k − 1 2 )2 +ρ1+ 2α2 π ]). (2) Sturm-Liouville eigenvalue problems with Case (2), (5) and (6), if α1 ≥ 0, α2 ≥ 0, ρ1 − ρ0 + 2α1 π + 2α2 π ≤ 1 3 and λk ∈ [(k − 1)2 + ρ0, (k − 1)2 + ρ1 + 2α1 π + 2α2 π ]. 4. Comparison with previous algorithms 4.1. Comparison with the state-of-the-art deep learning methods. Con- sider the smallest eigenpairs problems φ′′(x) + λφ(x) = 0, (4.1) with Dirichlet boundary value conditions φ(0) = φ(π) = 0. (4.2) Ben-Shaul et al. [5] utilized a more intricate cost function to simultaneously iden- tify the four smallest eigenpairs. In Table 4.1, Ben-Shaul1 denotes the simultaneous training of six smallest eigenpairs using 60000 epochs with their approach. Ben- Shaul2 represents the concurrent training of seven smallest eigenpairs using 70000 epochs. We find that their algorithm has a large error when n = 6, and they ob- tained incorrect results when n = 7. In contrast, our objective differs significantly. Our aim is to devise an easy-to-understand algorithm and adapt it for handling larger eigenvalues. Our method involves sequentially training 10 smallest eigen- pairs, with each pair trained using 10000 epochs. Note that it is a fair comparison, the total number of epochs used in our method is equivalent to that of the other methods. We set the integral point to 500 and lr = 0.01. The hyper-parameters are set as in [5]. We found that the runtime of our method is less than both Ben- Shaul1 and Ben-Shaul2, but it obviously has more accuracy, see Table 4.1. For our method, if we set 30000 epochs for each, 1000 inner points and 200 integral points, we will get more accurate eigenpairs. Table 4.1. Comparison with Ben-Shaul et al.’s method. Exact Value Ben-Shaul1 Ben-Shaul2 Ours (Runtime) λ1 = 1 0.9859 0.8639 1.0000 (39s) λ2 = 4 3.9042 1.0964 4.0000 (81s) λ3 = 9 8.8453 4.6499 8.9943 (149s) λ4 = 16 15.4377 8.5765 16.0017 (270s) λ5 = 25 24.4220 16.2587 24.9451 (338s) λ6 = 36 34.8055(1485s) 23.5614 35.9928 (404s) λ7 = 49 N/A 35.5628(2128s) 48.9332 (475s) λ8 = 64 N/A N/A 63.9604 (552s) λ9 = 81 N/A N/A 80.8455 (636s) λ10 = 100 N/A N/A 99.8336 (725s) EJDE-2024/53 DEEP LEARNING METHOD FOR FINDING EIGENPAIRS 13 4.2. Comparison with classical numerical methods. The following finite dif- ference method’s discretization scheme is derived from reference [1]. When a1 = a2 = 0, we have Dirichlet boundary value condition. φ(0) = φ(π) = 0. (4.3) If, at the internal grid points of the uniform grid G = {xi;xi = ih, i = 0, 1, . . . , n, n+ 1, h = π/(n+ 1)}, (4.4) the operator Lφ is approximated using central differences, with L =  2 h2 + ηρ(x1) − 1 h2 − 1 h2 2 h2 + ηρ(x2) − 1 h2 . . . . . . . . . − 1 h2 2 h2 + ηρ(xn−1) − 1 h2 − 1 h2 2 h2 + ηρ(xn)  . (4.5) Such a matrix is symmetric, ensuring the orthogonality of different eigenfunctions. The error of the k-th eigenvalue λk is about O(k4h2). However, for other boundary value conditions, especially those involving derivatives(a1 ̸= 0 or a2 ̸= 0), loss of symmetry brings additional difficulties when using classical numerical methods. When a1 ̸= 0 or a2 ̸= 0, if we use finite difference method, we found that it is not symmetric, resulting in the eigenfunctions not being orthogonal, while the original differential operator is self-adjoint, ensuring the orthogonality of the real eigenfunctions. Specifically, the central difference eigenvalues λk, k = 1, . . . , n + 1 of the above equation satisfy on the augmented grid G∗ = {xj ;xj = jh, j = −1, 0, 1, . . . , n, n+ 1, h = π/n}, (4.6) the following eigenvalue problem of matrix L =  2 h2 + ηρ(x0) − α h2 − 1 h2 2 h2 + ηρ(x1) − 1 h2 . . . . . . . . . − 1 h2 2 h2 + ηρ(xn−1) − 1 h2 − β h2 2 h2 + ηρ(xn)  (4.7) with α = 2a1 a1 + h ( b1 + a1 2 ρ′(0) ρ(0) ) and β = 2a2 a2 + h ( b2 − a2 2 ρ′(π) ρ(π) ) . Obviously, it is not a diagonal matrix. Although the eigenvalues we calculated are correct, the corresponding eigenfunctions are not orthogonal, which limits the application of such a method. Despite the fact that the study of the Dirichlet eigenvalue problem using deep learning methods is a hot topic, to our knowledge, there is no work on the boundary value condition involving derivatives using the deep learning method. Indeed, the method we present in this paper can easily tackle boundary value conditions with derivatives, resulting in orthogonal eigenfunctions. This is a very important advan- tage of deep learning methods that has not yet been noticed. Additionally, we don’t require ρ(x) to have a second-order derivative, which is an assumption required by the above method. Unlike classical numerical methods, our neural network solution provides an approximate analytical solution, in a certain sense, which allows us to obtain derivatives using automatic differentiation. 14 S. ZHANG, J. ZU, J. ZHANG EJDE-2024/53 5. Numerical experiments In this section, we will utilize the designed neural network algorithm to identify the smallest eigenpairs under various boundary value conditions. The numerical results obtained from the neural network algorithm will be compared to those from the classical numerical method, demonstrating the efficiency of our approach. More- over, we illustrate that our proposed method can be seamlessly extended to handle two-dimensional cases and periodic scenarios. Example 5.1. The Sturm-Liouville problem with Neumann boundary value condi- tions: −(ρ(x)φ′(x))′ = λρ(x)φ(x), x ∈ (0, π), φ′(0) = φ′(π) = 0, (5.1) where ρ(x) = 1 + 0.3 sinx. Table 5.1 shows the eigenvalues of the Sturm-Liouville Problem (5.1). λ̃n rep- resents the numerical results obtained via the finite difference method (4.5) with h = π 1000 . λ̂n represents the numerical results obtained from the deep learning algorithm. We get the smallest eigenvalues and eigenfunctions simultaneously, see Figure 5.1. Table 5.1. Eigenvalues of Sturm-Liouville Problem (5.1) λn λ1 λ2 λ3 λ4 λ5 λ6 λ7 λ8 λ9 λ10 λ̃n 7.56 · 10−6 1.1129 4.1098 9.1099 16.1119 25.1151 36.1193 49.1243 64.1302 81.1368 λ̂n 6.62 · 10−6 1.1126 4.1082 9.1077 16.1057 25.1079 36.0993 49.1056 64.0550 81.0975 Figure 5.1. Eigenfunctions of the Sturm-Liouville Problem (5.1) Example 5.1 shows that our neural network method is well adapted to the in- homogeneous media case, especially to the derivative boundary value case. The resulting eigenfunction is orthogonal and has broad prospects for use in other math- ematical problems. Example 5.2. Consider the periodic eigenpair problems: −φ′′(x) = λφ(x), x ∈ (0, π), φ(0) = φ(π), φ′(0) = φ′(π). (5.2) EJDE-2024/53 DEEP LEARNING METHOD FOR FINDING EIGENPAIRS 15 Periodic phenomena are ubiquitous in physics. Regarding the periodic eigen- value problem, where each eigenvalue corresponds to a pair of mutually orthogonal eigenvectors, to our knowledge, no scholar has yet studied how to use deep learning methods to find multiple smallest eigenpairs for such problems. In Table 5.2, λ̂n represents the numerical results obtained from the deep learning algorithm, and λn represents the exact value of the eigenvalue. It shows that our algorithm can be easily applied to periodic boundary value cases with enough accuracy. Table 5.2. Eigenvalues of periodic boundary value problems (5.2) λn λ1 λ2 λ3 λ4 λ5 λ6 λ7 λ8 λ9 λ10 λn 0 4 4 16 16 36 36 64 64 100 λ̂n 0.0000 4.0000 3.9999 15.9998 16.0000 36.0002 36.0001 64.0007 63.9998 99.9998 Example 5.3. 2D Sturm-Liouville problem with Robin boundary value conditions [10]: −∇(ρ(x, y)∇φ(x, y)) = λρ(x, y)φ(x, y), (x, y) ∈ Ω φ(0, y) = φx(π, y) = φ(x, 0) = φy(x, π) = 0, (5.3) where ρ(x, y) = (1 + 0.3 sinx)(1 + 0.3 sin y). In this example, the finite element method is used as the baseline to compare our Neural network method. The number of finite elements is set as 47905 by using the software of Freefem++. Usually, the finite elements method has very high accuracy. In Table 5.3, λ∗n represents the numerical eigenvalues obtained by finite element method and λ̂n represents the numerical results obtained from our deep learning algorithm. In Figure 5.2, we show the mutual orthogonal eigenfunctions that are learned simultaneously with eigenvalues. Such eigenfunctions are hard to learn for a non-expert using classical numerical methods. But our method is beginner-friendly. Table 5.3. Eigenvalues of 2D Sturm-Liouville problem (5.3) λn λ1 λ2 λ3 λ4 λ5 λ6 λ7 λ8 λ9 λ10 λ∗ n 0.4959 2.5067 2.5067 4.5175 6.5069 6.5069 8.5177 8.5177 12.5069 12.5069 λ̂n 0.4959 2.5067 2.5067 4.5174 6.5068 6.5068 8.5176 8.5176 12.5075 12.5071 The final example demonstrates that our neural network possesses considerably good accuracy. Our neural network can be easily utilized, whereas the finite ele- ment method poses significant difficulties for non-experts. How to further improve accuracy of neural network methods remains an important question for future stud- ies. 6. Conclusion Deep learning techniques have shown significant advancements in addressing eigenvalue problems, particularly with the use of PINN (Physics-Informed Neural Networks) due to its ease of implementation and ability to maintain operator self- adjointness. PINN has been applied to solve quantum problems and Sturm-Liouville 16 S. ZHANG, J. ZU, J. ZHANG EJDE-2024/53 Figure 5.2. Eigenfunctions of 2D Sturm-Liouville problem with (5.3) boundary value problems, but with limitations such as one-dimensional homoge- neous media or specific boundary value conditions. This paper introduces a novel deep learning method to find the smallest eigenpairs by constructing a new cost function, incorporating adaptive hyper-parameter tuning, and sequentially training eigenpairs. The proposed approach extends to varying boundary value conditions in inhomogeneous media, two-dimensional cases, and periodic scenarios, demon- strating improved simplicity, accuracy, and interpretability. Quantitative estima- tion of eigenpairs is given for the Sturm-Liouville eigenvalue problems. Numerical experiments validate the method’s efficiency and accuracy compared to previous algorithms [5]. References [1] R. Anderssen, F. De Hoog; On the correction of finite difference eigenvalue approximations for Sturm-Liouville problems with general boundary conditions, BIT Numerical Mathematics, 24 (1984), 401–412. [2] A. Andrew, J. Paine; Correction of finite element estimates for Sturm-Liouville eigenvalues, Numerische Mathematik, 50 (1986), 205–215. [3] V. Barbu, N. Pavel; Periodic solutions to nonlinear one dimensional wave equation with x-dependent coefficients, Transactions of the American Mathematical Society, 349 (1997), 2035–2048. [4] K. Beauchard, M. Mirrahimi; Practical stabilization of a quantum particle in a one- dimensional infinite square potential well, SIAM Journal on Control and Optimization, 48 (2009), 1179–1205. [5] I. Ben-Shaul, L. Bar, D. Fishelov, N. Sochen; Deep learning solution of the eigenvalue problem for differential operators, Neural Computation, 35 (2023), 1100–1134. [6] R. Chen, Y. Rubanova, J. Bettencourt, D. Duvenaud; Neural ordinary differential equations, NeurIPS 2018, 31 (2018). [7] S. Greydanus, M. Dzamba, J. Yosinski; Hamiltonian neural networks, NeurIPS 2019, 32 (2019). [8] E. Holliday, J. Lindner, W. Ditto; Solving quantum billiard eigenvalue problems with physics- informed machine learning, AIP Advances, 13 (2023), 085013. [9] H. Jin, M. Mattheakis, P. Protopapas; Physics-informed neural networks for quantum eigen- value problems, In IJCNN at IEEE World Congress on Computational Intelligence (2022). [10] S. Ji; Periodic solutions of two-dimensional wave equations with x-dependent coefficients and Sturm-Liouville boundary conditions, Nonlinearity, 35 (2022), 5033. [11] S. Ji, Y. Li; Periodic solutions to one-dimensional wave equation with x-dependent coeffi- cients, Journal of Differential Equations, 229 (2006), 466–493. EJDE-2024/53 DEEP LEARNING METHOD FOR FINDING EIGENPAIRS 17 [12] P. Jin, Z. Zhang, I. Kevrekidis, G. Karniadakis; Learning Poisson systems and trajectories of autonomous systems via Poisson neural networks, IEEE Transactions on Neural Networks and Learning Systems, 34 (2022), 8271–8283. [13] P. Jin, Z. Zhang, A. Zhu, Y. Tang, G. Karniadakis; SympNets: Intrinsic structure-preserving symplectic networks for identifying Hamiltonian systems, Neural Networks, 132 (2020), 166– 179. [14] G. Karniadakis, I. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, L. Yang; Physics-informed machine learning, Nature Reviews Physics 3 (2021), 422–440. [15] L. Lu, P. Jin, G. Pang, Z. Zhang, G. Karniadakis; Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators, Nature Machine Intelligence, 3 (2021), 218–229. [16] L. Lu, X. Meng, Z. Mao, G. Karniadakis; DeepXDE: A deep learning library for solving differential equations, SIAM Review 63 (2021), 208–228. [17] J. Paine, R. Anderssen, F. De Hoog; On the correction of finite difference eigenvalue approx- imations for sturm-liouville problems, Computing, 26 (1981), 123–139. [18] G. Pang, L. Lu, G. Karniadakis; fPINNs: Fractional physics-informed neural networks, SIAM Journal on Scientific Computing, 41 (2019), A2603–A2626. [19] M. Raissi, P. Perdikaris, G. Karniadakis; Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, Journal of Computational physics, 378 (2019), 686–707. [20] L. Yang, D. Zhang, G. Karniadakis; Physics-informed generative adversarial networks for stochastic differential equations, SIAM Journal on Scientific Computing, 42 (2020), A292– A317. [21] Y. Zang, G. Bao, X. Ye, H. Zhou; Weak adversarial networks for high-dimensional partial differential equations, Journal of Computational Physics, 411 (2020), 109409. [22] J. Zu; Approximate stabilization of one-dimensional Schrödinger equations in inhomogeneous media, Journal of Optimization Theory and Applications, 153 (2012), 758–768. Sen Zhang Center for Mathematics and Interdisciplinary Sciences, and School of Mathematics and Statistics, Northeast Normal University, Changchun, 130024, China Email address: zhangs832@nenu.edu.cn Jian Zu (corresponding author) Center for Mathematics and Interdisciplinary Sciences, and School of Mathematics and Statistics, Northeast Normal University, Changchun, 130024, China Email address: zuj100@nenu.edu.cn Jingqi Zhang Center for Mathematics and Interdisciplinary Sciences, and School of Mathematics and Statistics, Northeast Normal University, Changchun, 130024, China Email address: zhangjq906@nenu.edu.cn 1. Introduction 2. Model problem 3. Quantitative estimation of eigenpairs 4. Comparison with previous algorithms 4.1. Comparison with the state-of-the-art deep learning methods 4.2. Comparison with classical numerical methods 5. Numerical experiments 6. Conclusion References