1 Introduction

We shall consider a linear Schrödinger equation as follows. Let be a bounded rectangular-type domain with a smooth boundary ∂Ω. We find a complex-valued function u(x, t) defined on Ω×[0, T] and satisfying

(1)

where u₀(x) is a given initial complex-valued function, and the trapping potential function V(x) is non-negative bounded and real-valued.

The Schrödinger equation is an important equation in quantum mechanics. There are many numerical methods to solve the Schrödinger equation in the literature, such as the spectral method^[1-2], the finite difference method^[3-5], the finite element method^[6-12], the discontinuous Galerkin method^[13-15], and the local discontinuous Galerkin method^[16-18]. Bao et al.^[1] studied the performance of time-splitting spectral approximations for the general nonlinear Schrödinger equation in the semiclassical regimes. Han et al.^[5] introduced an artificial boundary condition to reduce the one-dimensional time-dependent Schrödinger equation into an initial-boundary value problem in a finite computational domain. Antonopoulou et al.^[7] considered an initial and boundary-value problem for a general Schrödinger-type equation posed on a two space-dimensional noncylindrical domain with mixed boundary conditions. Karakashian and Makridakis^[14] analyzed the convergence of the discontinuous Galerkin method for the nonlinear Schrödinger equation. Guo and Xu^[16] presented a fully discrete scheme by discretizing the space with the local discontinuous Galerkin method and the time with the Crank-Nicholson scheme to simulate the multi-dimensional Schrödinger equation with wave operator.

Superconvergence has been studied for long. Many different numerical methods have been analyzed. It is a powerful tool to improve the approximation accuracy and efficiency. There are numerous studies by many famous scholars^[19-22]. At present, superconvergence results were obtained widely for elliptic, parabolic, Maxwell's equations, and optimal control problems^[23-31]. However, there were not many superconvergence results for the Schrödinger equation^[32-36]. In 1998, Lin and Liu^[32] studied a time-dependent linear Schrödinger equation and analyzed the superconvergence error results. In 2014, Shi et al.^[33] considered a nonlinear Schrödinger equation by the finite element method in the triangular anisotropic meshes and proved the superconvergence result in the semi-discrete scheme. Later, Wang et al.^[35] conducted the superconvergence analysis for a time-dependent Schrödinger equation by using the interpolation operator and obtained the error result in the H¹ norm with O(h^p+1) in the semi-discrete scheme and in the Crank-Nicolson scheme, respectively. Recently, Zhou et al.^[36] studied the superconvergence properties of the local discontinuous Galerkin method for the one-dimensional linear Schrödinger equation.

In this paper, we study a general complex linear Schrödinger equation (1) and extend the previous work^[35]. We analyze the error estimate using the elliptic projection operator. We obtain the error result with O(h^k+1) in the L² norm and the H¹ norm in the semi-discrete finite element scheme. The global superconvergence result is presented by use of the interpolation post-processing technique. Next, we analyze the error estimate in the L² norm with order O(h^k+τ²) in the Crank-Nicolson fully discrete scheme. We extend the idea^[37] and certify that the time-difference of error ηⁿ=Uⁿ-P_huⁿ has a high order error in the L² norm, that is, ||ηⁿ-η^n-1|| ≤ Cτ(h^k+1+τ²), where Uⁿ is the fully discrete solution of Crank-Nicolson scheme. At last, we obtain the superconvergence result in the H¹ norm with O(h^k+τ²) on this basis.

The paper is organized as follows. The notations and the projection operator are given in Section 2. In Section 3, we present a finite element semi-discrete scheme with bi-k-degree rectangular elements. Furthermore, we obtain error results with O(h^k) in the L² norm and the H¹ norm by use of the elliptic projection operator, respectively. In Section 4, we prove the global superconvergence result with O(h^k). In Section 5, we obtain the superconvergence result in the H¹ norm with O(h^k+1τ²) in the Crank-Nicolson fully discrete scheme. In Section 6, numerical examples are given to partly verify the theoretical results.

2 Notation and preliminaries

For an integer m≥ 0 and 1≤ p ≤∞, we shall use W^{m, p} to denote the standard Sobolev space of complex-valued measurable functions defined on Ω with the norm . When p=2, we shall also use the symbol H^m for W^{m, 2}, || · ||_m instead of || · ||_{m, 2}, and || · || instead of || · ||_{0, 2}.

For complex-valued functions ω(x) and υ(x), we define the inner product (ω, υ) with

where υ denotes the complex conjugate of function υ.

Then, we can define the weak solution u(x, t) of problem (1): find a function u(x, t)∈ H₀¹(Ω) such that

(2)

where .

Let Γ_h be a quasi-uniform rectangular partition of Ω with the mesh size h>0, and let e be an arbitrary element of Γ_h. We can define the finite element space of order k as

where

In addition,

Let V₀^{h, k} ⊂ H₀¹(Ω) be the corresponding finite element space of order k. In general given w(x, t) ∈ H₀¹(Ω), the elliptic projection P_hw(x, t) ∈ V₀^{h, k} can be defined by

(3)

Let τ =T/N be the time step of the interval [0, T], time nodes t_j=jτ(j=0, 1, ·, N), , and time elements I_j=[t_j, t_j+1](j=0, 1, ·, N^-1), and set

3 Superconvergence analysis for semi-discrete approximation problem

The semi-discrete finite element solution u_h(x, t) of problem (1) can be defined: find u_h(x, t) ∈ V₀^{h, k} satisfying

(4)

where P_hu₀(x)∈ V₀^{h, k} is the elliptic projection of u₀(x).

Lemma 1^[34]   If for any t ∈ [0, T], the functions u(x, t), u_t(x, t), u_tt(x, t) ∈ H^k+1(Ω), then P_hu(x, t)∈ V₀^{h, k} has the following results:

(5)

(6)

(7)

Lemma 2^[20]   Let u be the solution to the problem (2), and let u_I∈ V₀^{h, k} be the interpolation of u. If u∈ H^k+2(Ω), then

(8)

Theorem 1   If u and u_h are the solutions to the problems (2) and (4), respectively, and u, u_t, u_tt∈ H^k+1(Ω), there hold

(9)

(10)

Proof   It follows from (2) and (4) that

(11)

Let u-u_h=ρ-ξ with

(12)

Then, from (11) and (12), we have

(13)

From (3), we can obtain

(14)

Substituting (14) into (13) yields

(15)

Taking v_h=ξ in (15), we have

(16)

Noticing

and comparing the imaginary parts of (16), we get

(17)

Combining (6) with (17) yields

(18)

Integrating from 0 to t in (18), we have

(19)

It follows from (4) that

(20)

From (19) and (20), we obtain

Therefore, (9) holds. Next, we prove (10). Taking v_h=ε_t(·, 0) in (15) with t=0 and combining (20), we have

Thus,

that is,

(21)

Combining (6) with (21) gives

(22)

Differentiating (15) with respect to t and taking v_h=ξ_t, we can obtain

(23)

Noticing

and comparing the imaginary parts of (23) yield

(24)

From (22) and (25), we get

(25)

Integrating from 0 to t in (25), we can obtain

(26)

It follows from (22) and (26) that

which completes the proof of (10).

Theorem  2   Let u and u_h be the solutions to the problems (2) and (4), respectively, and u, u_t, u_tt∈ H^k+1(Ω). Then, we have

(27)

Proof   Taking v_h=ξ_t in (15), we can get

that is,

(28)

Noticing

and comparing the real parts of (28) give

(29)

From (6), (9), (10), and (29), we can obtain

(30)

Integrating from 0 to t in (30) yields

(31)

Notice

(32)

Substituting (32) into (31) yields

(33)

Therefore, (27) follows from (33).

4 Global superconvergence analysis

Let be a macro element which is the union of four elements e_i ∈ Γ_h (i=1, 2, 3, 4), where the intersection of e_i ∈ Γ_h(i=1, 2, 3, 4) is nonempty (see Fig. 1).

Let the interpolation operator Π_2h² satisfy Π_2h²w ∈ Q₂(), where Q₂ is the space of biquadratic functions, and

(34)

where Z_i(i=1, 2, ·, 9) are the nodes of Γ_h.

When k≥ 2, let Π_2h^2kw ∈ Q_2k() such that

(35)

(36)

(37)

where Z_i(i=1, 2, ·, 9) are the nodes of Γ_h, l_i(i=1, 2, ·, 12) are the edges of Γ_h, e_i(i=1, 2, 3, 4) are the elements of Γ_h, P_k-2 is the set of polynomials of order k-2, and Q_k-2() is the polynomials of order k-2 in x and y.

Lemma 3^{[20, 22]}   The interpolation operator Π_2h^2k is defined in (34)-(37) such that

(38)

(39)

(40)

where w_I∈ V^{h, k} is the interpolant of w.

Lemma 4   Let u and u_h be the solutions to the problems (2) and (4), respectively. If u∈ H^k+2(Ω), and u_t, u_tt∈ H^k+1(Ω), then

(41)

where u_I is the interpolant of u.

Proof   From (14), we can obtain

that is,

(42)

It is easy to check

(43)

Combining (42) with (43) yields

(44)

It follows from (5) that

(45)

and from (8), we have

(46)

In addition,

(47)

Substituting (27) and (45)-(47) into (44), we can get

(48)

By the Poincaré inequality, we can obtain

(49)

Therefore, (48) and (49) show the validity of (41).

Theorem 3   Let u and u_h be the solutions to the problems (2) and (4), respectively. If u∈ H^k+2(Ω), and u_t, u_tt∈ H^k+1(Ω), then

(50)

where Π_2h^2k is the interpolation post-processing operator.

Proof   It follows from (40) and (41) that

(51)

From (38) and (39), we have

(52)

Notice

(53)

Therefore, (50) follows from (51)-(53).

5 Superconvergence analysis in fully discrete scheme

For the function series Uⁿ(x)(n=0, 1, ·), let

Then, the Crank-Nicolson fully discrete finite element solution Uⁿ(x)∈ V₀^{h, k}(n=0, 1, ·, N) to the problem (1) can be defined by

(54)

Theorem 4   Let u(x, t) be the solution to the problem (2), and let the function series Uⁿ(x) be the solution to the problem (54). Then, we have

(55)

Proof   From (2) and (54), we can get

(56)

Let u-U=ρ-η with

(57)

Combining (56) with (57) and (14), we have

(58)

Taking in (58), we can obtain

(59)

Notice

Comparing the imaginary parts of (59) yields

Thus,

(60)

It follows from (5) that

(61)

In addition,

(62)

Substituting (61) and (62) into (60), we have

(63)

Summing up for n in (63) yields

(64)

From (54), we can see

(65)

Therefore, (55) follows from (64) and (65).

Lemma 5   Let u(x, t) be the solution to the problem (2), and let the function series Uⁿ(x) be the solution to the problem (54). Then, the time-difference of error ηⁿ=Uⁿ-P_huⁿ has a high order error

(66)

Proof   It follows from (2) that

(67)

Integrating (67) in I_n by trapezoid and mid-point formulae, respectively, we can obtain

(68)

where

(69)

(70)

From (54), we have

(71)

Combining (68) with (71) yields

(72)

From (57), (14), and (72), we get

(73)

where

(74)

Further, combining (7) and (74) gives

(75)

Substituting n by n-1 in (73), we have

(76)

Let

We can see

(77)

Subtracting (76) from (73) and combining (77) yield

(78)

where

(79)

From (69), (70), (75), and (79), we can obtain

(80)

Taking v=ϵⁿ⁺¹+ϵⁿ in (78), we get

(81)

Comparing the imaginary parts of (81), we have

(82)

Combining (80) with (82) gives

(83)

Without loss of generality, we assume that there is an integer 1≤ K ≤ N such that

(84)

Summing up for n from 2 to K in (83) and combining (84), we have

(85)

Taking n=1 in (64) and combining (65) yield

(86)

Substituting (86) into (85) and using Young's inequality, we can get

(87)

Therefore, (66) follows from (84) and (87).

Theorem  5   Let u(x, t) be the solution to the problem (2), and let the function series Uⁿ(x) be the solution to the problem (54). Then, we have

(88)

Proof   Taking in (58), we have

(89)

Notice

Comparing the real parts of (89), we get

that is,

(90)

Summing up for n in (90) and combining (65), we have

(91)

Substituting (61), (66), (62), and (55) into (91), we can obtain

that is,

which completes the proof.

Similar to Theorem 3, we can obtain the following result.

Theorem  6   Assume that u(x, t) is the solution to the problem (2), and the function series Uⁿ(x) is the solution to the problem (54). Then, we have the global superconvergence estimate

(92)

where Π_2h^2k is the interpolation post-processing operator.

6 Numerical examples

In this section, we carry out some numerical examples with k=1 and k=2 to demonstrate the validity of the theoretical analysis.

Example 1   We consider the following linear Schrödinger equation:

(93)

where Ω =[0, 1]×[0, 1], and let the function f(x, t) be chosen that

is the exact solution.

We have solved the Schrödinger equation on the uniformly rectangular meshes with the mesh size h by the bilinear finite element. First, we calculate the errors with fixing τ =10^-4 by varying h. The error results are presented in Tables 1-4, where Order₁, Order₂, Order₃, and Order₄ denote the convergence orders of || u_I-Uⁿ ||, || u-Uⁿ ||₁, || u_I-Uⁿ||₁, and || u-Π_2h²Uⁿ ||₁, respectively. Moreover, we have shown convergence orders by slopes in Figs. 2-5. Results in all tables show O(h) in || u-Uⁿ ||₁, and O(h²) convergence rate clearly in || u_I-Uⁿ ||, || u_I-Uⁿ ||₁, and || u-Π_2h²Uⁿ ||₁.

To test the convergence rate in terms of τ, we fix the time step τ=h. The error results are shown in Tables 5 and 6. In addition, we also show the convergence orders by slopes in Figs. 6 and 7. Results show the convergence rate O(τ²) clearly in || u_I-Uⁿ ||, || u_I-Uⁿ ||₁, and || u-Π_2h²Uⁿ ||₁.

Example 2   We consider the problem (93) with Ω =[-1, 1]×[-1, 1], and function f(x, t) is chosen corresponding to the exact solution

Similarly, we have solved the Schrödinger equation by the bilinear finite element. We calculate the errors with fixing τ =10^-4 by varying h. The error results at the time level t_n=0.01, 0.1, 0.5, 1.0 are presented in Tables 7-10, respectively. Results in all tables show O(h) in || u-Uⁿ ||₁, and O(h²) convergence rate clearly in || u_I-Uⁿ ||, || u_I-Uⁿ ||₁, and || u-Π_2h²Uⁿ ||₁.

Then, we take the time step τ=h. The error results are listed in Tables 11 and 12. Results show the convergence rate O(τ²) clearly in || u_I-Uⁿ ||, || u_I-Uⁿ ||₁, and || u-Π_2h²Uⁿ ||₁ as well, which are coincident with theoretical results.

The profiles of the exact solution and the numerical solution at t=1.0 on the 64 × 64 mesh grid are plotted in Figs. 8-11.

Example 3   We consider the problem (93) with Ω =[-1, 1]×[-1, 1], and function f(x, t) is chosen corresponding to the same exact solution with Example 2.

The domain Ω is uniformly divided into families Γ_h of quadrilaterals with mesh size h, and V^{h, 2} is the biquadratic rectangular element space defined on Γ_h. The Schrödinger equation is solved by the biquadratic rectangular element. We calculate the errors with fixing τ =10^-3 by varying h. The error results at time level t_n=0.1, 0.2, 0.5, 1.0 are presented in Tables 13-16, respectively. Results in all tables also show O(h²) in || u-Uⁿ ||₁, and O(h³) convergence rate clearly in || u-Uⁿ || and || u_I-Uⁿ ||₁, which are consistent with our theoretical analysis. In addition, the results show O(h⁴) in || u_I-Uⁿ ||. When k≥ 2, there is the superclose property also in the L² norm between the numerical solution with the interpolant of exact solution.

7 Conclusions

In this paper, we consider a two-dimensional time-dependent linear Schrödinger equation with the finite element method. We present the finite element semi-discrete scheme and the Crank-Nicolson fully discrete scheme in the rectangular Lagrange type finite element space of order k. We also obtain the superconvergence result in the H¹ norm by use of the elliptic projection in the semi-discrete scheme and the fully discrete scheme, respectively. Some numerical examples with the order k=1 and k=2 are provided to partly verify our theoretical results. In the future, we shall try to study the problem of superconvergence in the L² norm for the two-dimensional time-dependent Schrödinger equation and the superconvergence in the H¹ norm for the three-dimensional Schrödinger equation with the finite element method.

Acknowledgements We would like to thank anonymous referees for their insightful comments that improved this paper.