# An Efficient Finite Element Method and Error Analysis for Schrödinger Equation with Inverse Square Singular Potential

Authors: Hui Zhang; Fubiao Lin; Junying Cao
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

We provide in this study an effective finite element method of the Schrödinger equation with inverse square singular potential on circular domain. By introducing proper polar condition and weighted Sobolev space, we overcome the difficulty of singularity caused by polar coordinates’ transformation and singular potential, and the weak form and the corresponding discrete scheme based on the dimension reduction scheme are established. Then, using the approximation properties of the interpolation operator, we prove the error estimates of approximation solutions. Finally, we give a large number of numerical examples, and the numerical results show the effectiveness of the algorithm and the correctness of the theoretical results.


---

## Body

## 1. Introduction

Schrödinger equation with the inverse square or centrifugal potential plays an important role in quantum mechanics, quantum cosmology, nuclear physics, molecular physics, and so on [1–8]. The potential has the same differential order as the Laplacian operator near the origin, which usually leads to strong singularities and cannot be treated as a lower-order perturbation term [9–14]. Li et al. [15] proposed an efficient finite element method to discuss the numerical solution of time-fractional Schrö dinger equations. Thus, we need to develop some new numerical methods to solve Schrödinger equation with inverse square singular potential.In recent years, more and more attention has been paid to the numerical methods of the schrödinger equations with similar singular potential [1, 16–21]. However, many numerical methods are based on low-order finite element methods. If we solve these problems directly in two-dimensional domain, it will cost a lot of computing time and memory capacity to obtain high-precision numerical solutions [22–24]. In practice, we usually need to solve the Schrödinger equation with inverse square singular potential on circular domain. As far as we know, there are few reports on an effective numerical method for the Schrödinger equation with inverse square potential in circular domain. Thus, the purpose of this paper is to propose an effective finite element method of the Schrödinger equation with inverse square singular potential on circular domain. By introducing proper polar condition and weighted Sobolev space, we overcome the difficulty of singularity caused by polar coordinates transformation and singular potential and establish the weak form and corresponding discrete scheme based on the dimension reduction format. Then, using the approximation properties of interpolation operator, we prove the error estimates of approximation solutions. Finally, we give a large number of numerical examples, and the numerical results show the effectiveness of the algorithm and the correctness of the theoretical results [25, 26].The rest of this paper is organized as follows. InSection 2, we derive an equivalent scheme based on variable separation. In Section 3, we prove the existence and uniqueness of the solution. In Section 4, we prove the error estimation of approximation solutions. In Section 5, we describe the details for an efficient implementation of the algorithm. In Section 6, we provide some numerical experiments to show the accuracy and efficiency of our algorithm. Finally, in Section 7, we give some concluding remarks.


## 2. An Equivalent Scheme Based on Variable Separation

We are interested in studying the following Schrödinger equation with inverse square singular potential:(1)−Δψx,y+βrr2ψx,y=fx,y,inD,(2)ψx,y=0,on∂D,where 0<β∗≤βr≤β∗ and D=x∈ℝ2:0≤r<R with r=x2+y2. Let x=rcosθ,y=rsinθ, Ur,θ=ψx,y, and Fr,θ=fx,y. Then, the Laplace operator in polar coordinates is as follows:(3)ℒU=1r∂∂rr∂U∂r+1r2Uθθ.We can rewrite (1) and (2) as follows:(4)−rℒUr,θ+βrrUr,θ=rFr,θ,r,θ∈0,R×0,2π,(5)UR,θ=0.SinceUr,θ and Fr,θ are 2π periodic in θ, then we have(6)Ur,θ=∑|m|=0∞umreimθ,Fr,θ=∑|m|=0∞fmreimθ,where um and fm are the Fourier coefficients of Ur,θ and Fr,θ, respectively. We can derive from (3) and (6) that(7)rℒUr,θ=∑|m|=0∞∂∂rr∂umr∂r−m2rumreimθ,(8)Ur,θr=∑|m|=0∞umrreimθ.To make (7) and (8) meaningful, we need introduce the following essential pole conditions:(9)m2um0=0,um0=0.The pole condition (9) can be further reduced to(10)um0=0.Using the orthogonal properties of Fourier basis functions and polar condition (10), we can reduce (4) and (5) to a series of equivalent one-dimensional Schrödinger equations as follows:(11)−rℒmumr+βrrumr=rfmr,r∈0,R,(12)um0=umR=0,where(13)ℒumr=1r∂∂rr∂umr∂r−m2r2umr.


## 3. Existence and Uniqueness of the Solution

For convenience, we use the expressiona≲b to mean that a≤Cb, where C is a positive constant. In order to derive the weak form and corresponding discrete scheme of equations (11) and (12), we need to introduce the usual weighted Sobolev space:(14)Lω20,R≔u:∫0Rωu2dr<∞,with the corresponding inner product and norm,(15)u,vω=∫0Rωuvdr,uw=∫0Rωu2dr1/2,and the nonuniformly weighted Sobolev space:(16)ℋ0,ω10,R≔um:∂rkum∈Lω2k−120,R,k=0,1,um0=umR=0,with the corresponding inner product and norm,(17)um,vm1,ω=∫0Rr∂rum∂rvm+1rumvmdr,um1,ω=um,um1,ω1/2,where ω=r is a weight function. Then, the weak form of equations (11) and (12) is to find um∈ℋ0,ω10,R, such that(18)amum,vm=Fmvm,∀vm∈ℋ0,ω10,R,where(19)amum,vm=∫0Rr∂rum∂rvm+m2+βrrumvmdr,Fmvm=∫0Rrfmvmdr.Define approximation spaceXh=Ph∩ℋ0,ω10,R, where Ph is a piecewise linear interpolation polynomial space. Then, the corresponding discrete scheme of (18) is to find umh∈Xh, such that(20)amumh,vmh=Fmvmh,∀vmh∈Xh.Lemma 1.
amum,vm is a continuous and coercive bilinear functional on ℋ0,ω10,R×ℋ0,ω10,R, i.e.,(21)amum,vm≲um1,ωvm1,ω,amum,um>rsimum1,ω2.Proof.
We derive from the Cauchy–Schwarz inequality that(22)amum,vm=∫0Rr∂rum∂rvm+m2+βrrumvmdr≤∫0Rr∂rum∂rvm+m2+β∗rumvmdr≲∫0Rr∂rum∂rvm+1rumvmdr≲∫0Rr∂rum2dr1/2∫0Rr∂rvm2dr1/2+∫0R1rum2dr1/2∫0R1rvm2dr1/2≲∫0Rr∂rum2+1rum2dr1/2∫0Rr∂rvm2+1rvm2dr1/2≲um1,ωvm1,ω,amum,um=∫0Rr∂rum2+m2+βrrum2dr≥∫0Rr∂rum2+m2+β∗rum2dr≲∫0Rr∂rum2+1rum2dr≲um1,ω2.Lemma 2.
Iffm∈Lω20,R, then Fmvm is a bounded linear functional on ℋ0,ω10,R, i.e.,(23)Fmvm≲um1,ω.Proof.
From Cauchy–Schwarz inequality, we have(24)Fmvm=∫0Rrfmvmdr⩽∫0Rrfm2dr1/2∫0Rrvm2dr1/2⩽fmω∫0Rr∂rum2+R2rvm2dr1/2≲um1,ω.This finishes our proof.Theorem 1.
Iffm∈Lω20,R, then problems (18) and (20) have unique solutions um and umh, respectively.Proof.
From Lemma1, we have(25)amum,vm≲um1,ωvm1,ω,amum,um≲um1,ω2,for ∀um,vm∈ℋ0,ω10,R. That is, amum,vm is a bounded and positive definite bilinear functional defined on ℋ0,ω10,R×ℋ0,ω10,R. In addition, from Lemma 2, we have(26)Fmvm≲um1,ω,which means Fmvm is a bounded linear functional defined on ℋ0,ω10,R. Then, from Lax–Milgram lemma, we know that equations (18) and (20) have unique solutions um and umh, respectively.


## 4. Error Estimation of Approximation Solutions

In this section, we will present the error estimates of approximate solutions. Define the piecewise linear interpolation operatorℐh:ℋ0,ω10,R⟶Xh by(27)ℐhumr=pmir,r∈Ii,where pmir is the linear interpolation polynomial of um on interval Ii=ri−1,ri. Let(28)umir=umr,r∈Ii.Then, from error formula of linear interpolating remainder term, we have(29)umir−pmir=∂r2umiξi2!r−ri−1r−ri,where ξi∈Ii is a function depending on r.Theorem 2.
LetGmir=∂r2umiξi/2!. Suppose that um is smooth enough such that ∂rkGmir≤Mk=0,1, where M is a large enough constant(0≤m≤M) in Fourier transformation. Then, it holds(30)∫0R∂rℐhumr−umr2dr≲h2,where h=max1≤i≤nhi, hi=ri−ri−1.Proof.
Since(31)umir−ℐhumir=umir−pmir=Gmirr−ri−1r−ri,then we derive that(32)∂rumir−ℐhumir=∂rGmirr−ri−1r−ri+Gmir∂rr−ri−1r−ri.
Then, we obtain(33)∂rumir−ℐhumir2≲r−ri−1r−ri2+∂rr−ri−1r−ri2=r−ri−1r−ri2+r−ri+r−ri−12≤hi24+2hi2≲hi2.
Thus, we derive that(34)∫0R∂rℐhumr−umr2dr=∑i=1n∫ri−1ri∂rℐhumir−umir2dr≲∑i=1nhi3≲h2.The proof of Theorem2 is complete.Lemma 3.
For anyumr∈ℋ0,ω10,R, the following inequality holds:(35)∫0R1rum2rdr≤R∫0R∂rumr2dr.Proof.
Sinceumr∈ℋ0,ω10,R, then we have(36)1rum2r=1r∫0r∂rumtdt2≤1r∫0r12dt∫0r∂rumt2dt=∫0r∂rumt2dt≤∫0R∂rumr2dr.
Thus, we derive that(37)∫0R1rum2rdr≤R∫0R∂rumr2dr.Theorem 3.
Letum and umh be the solutions of problems (18) and (20), respectively. Then, the following inequality holds:(38)um−umh1,ω≲h.Proof.
We derive from (18) and (20) that(39)amum,vmh=Fmvmh,∀vmh∈Xh,amumh,vmh=Fmvmh,∀vmh∈Xh.
Then, we have(40)amum−umh,vmh=0,∀vmh∈Xh.
We derive from Lemma1 and (40) that(41)um−umh1,ω2≲amum−umh,um−umh=amum−umh,um−vmh+vmh−umh=amum−umh,um−vmh≲umr−umh1,ωum−vmh1,ω.
Then, we obtain(42)um−umh1,ω≲um−vmh1,ω,∀vmh∈Xh.
Then, we have(43)um−umh1,ω2≲infvmh∈XNum−vmh1,ω2≲um−ℐhum1,ω2=∫0Rr∂rum−ℐhum2+1rum−ℐhum2dr.
We derive from Lemma3 that(44)um−umh1,ω2≲∫0Rr∂rum−ℐhum2+1rum−ℐhum2dr≲∫0RR∂rum−ℐhum2dr≲∫0R∂rum−ℐhum2dr.
Combining with Theorem2, we can obtain the desired result.


## 5. Implementation of the Algorithm

To solve the discrete scheme (20), we need to construct a set of basis functions of approximation space. Let(45)ψir=r−ri−1hi,ri−1≤r≤ri,−r−ri+1hi+1,ri≤r≤ri+1,0,others,where i=1,…,N−1. It is clear that(46)Xh=spanψ1r,…,ψN−1r.Set(47)aij=∫0Rrψjr′ψir′dr,bij=∫0R1rψjrψirdr,cij=∫0Rβrrψjrψirdr,fj=∫0Rrfmrψjrdr.Let(48)umh=∑k=1N−1umkψkr.Substituting expression (48) into (20) and taking vmh through all the basis functions in Xh, we obtain the following linear system:(49)A+m2B+ℂU=F,where(50)A=aij,B=bij,ℂ=cij,F=f1,…,fN−1T,U=u1,…,uN−1T.From the properties of the basis functions, we know that the stiff matrices and mass matrix in (49) are all tridiagonal sparse matrices.


## 6. Numerical Experiments

We will perform some numerical tests in order to show the convergence and the effectiveness of our algorithm. We operate our programs in MATLAB 2015b.LetψMhx,y be the approximation solution of exact solution ψx,y.By using the following polar coordinate transformation,(51)x=rcosθ,y=rsinθ,0≤r≤R,0≤θ<2π,we obtain(52)ψx,y=Ur,θ=∑|m|=0∞umreimθ,ψMhx,y=UMhr,θ=∑|m|=0Mumhreimθ.We define the error between the exact solutionψx,y and the approximation solution ψMhx,y as follows:(53)eψx,y,ψMhx,y=ψx,y−ψMhx,yL∞D=Ur,θ−UMhr,θL∞D.Example 1.
We takeβr=r2+4,R=1, ψx,y=xyx2+y2−1. It is clear that ψx,y satisfies the boundary condition (2). The fx,y can be obtained by substituting ψx,y into equation (1). Next, we solve problems (1) and (2) by using the algorithm proposed in this paper. We list the errors between exact solution and approximation solutions in Table 1 for different M and h. In order to further show the accuracy and convergence of our algorithm, we present the figures and their error figures of exact solution and approximation solution in Figures 1 and 2, respectively.
We observe from Table1 that the error eψx,y,ψMhx,y achieve about 10−4 with h≤1/32 and M=6. In addition, we can see from Figures 1 and 2 that the numerical solution converges to exact solution with the decrease of h.Table 1
The erroreψx,y,ψMhx,y for different M and h.
hM=6M=8M=10M=121/80.05170.05170.07540.07541/164.4437e − 044.4437e − 044.4437e − 044.4437e − 041/321.1380e − 041.1380e − 041.1380e − 041.1380e − 041/642.8623e − 052.8623e − 052.8623e − 052.8623e − 05Figure 1
Figures of the exact solution (a) and the numerical solution (b) withh=1/64,m=12.
(a)(b)Figure 2
The error figures of exact solution and numerical solution withh=1/32,m=10 (a) and h=1/64,m=12 (b).
(a)(b)Example 2.
We takeβr=2/3,R=1, and u=x2+y2−1sinx2+y2. It is obvious that ψx,y satisfies the boundary condition (2). Similarly, fx,y can be obtained by substituting ψx,y into equation (1). We list the errors between exact solution and approximation solutions in Table 2 for different M and h. In order to further show the accuracy and convergence of our algorithm, we present the figures and their error figures of exact solution and approximation solution in Figures 3 and 4, respectively.
We observe from Table2 that the error eψx,y,ψMhx,y achieves about 10−4 with h≤1/32 and M=6. In addition, we can see from Figures 3 and 4 that the numerical solution converges to exact solution with the decrease of h.Table 2
The erroreψx,y,ψMhx,y for different M and h.
hM=6M=8M=10M=121/80.00270.09780.09780.09781/166.5774e − 046.5774e − 046.5774e − 046.5774e − 041/321.6360e − 041.6360e − 041.6360e − 041.6360e − 041/644.0798e − 054.0798e − 054.0798e − 054.0798e − 05Figure 3
Figures of the exact solution (a) and the numerical solution (b) withh=1/64,m=12.
(a)(b)Figure 4
The error figures of numerical solution and exact solution withh=1/32,m=10 (a) and h=1/64,m=12 (b).
(a)(b)Example 3.
We takeβ=1/2,R=1, and ψx,y=x2+y2−1xex+y. We list the errors between exact solution and approximation solutions in Table 3 for different M and h. We also present the figures and their error figures of exact solution and approximation solution in Figures 5 and 6, respectively.
We observe from Table3 that the error eψx,y,ψMhx,y achieves about 10−4 with h≤1/64 and M=6. In addition, we can see from Figures 5 and 6 that the numerical solution converges to exact solution with the decrease of h.Table 3
The erroreψx,y,ψMhx,y for different M and h.
hM=6M=8M=10M=121/160.00220.00220.00220.00501/325.4324e − 045.4324e − 045.4324e − 045.4324e − 041/641.3583e − 041.3583e − 041.3583e − 041.3583e − 041/1284.1754e − 053.3967e − 053.3944e − 053.3944e − 05Figure 5
Figures of the exact solution (a) and the numerical solution (b) withh=1/128,m=12.
(a)(b)Figure 6
The error figures of numerical solution and exact solution withh=1/64,m=10 (a) and h=1/128,m=12 (b).
(a)(b)


## 7. Conclusions

We present in this paper an efficient finite element method for the Schrödinger equation with the inverse square potential on the circular domain. By using polar coordinate transformation and Fourier basis function expansion, we reduce the original problem into a series of equivalent one-dimensional problems. By introducing polar conditions, we overcome not only the difficulty brought by the singular potential but also the degree of freedom which is greatly reduced by dimension reduction. Thus, we only spend less computing time and memory capacity to obtain high-precision numerical solutions. Numerical results show that our algorithm is very effective. We mainly focus on, in this paper, the Schrödinger equation with the inverse square potential on the circular domain. In fact, we can extend our method to the Schrödinger equation with more complex potentials.



---

*Source: 8740634-2021-06-14.xml*