二、 Formulation
2.4 Time integration
To march the numerical solution in time, we adopt the Crank-Nicolson (CN) method [7] and Implicit-explicit Runge-Kutta (IMEX-RK) method [3, 8]. Denote the time step by ∆t and the n-th time level by tn = n∆t. Let un be the numerical solution at time tn for the following differential equation
du
dt = F (t, u(t)), (39a)
u(0) = u0. (39b)
The Crank-Nicolson algorithm solves Eqs. (39) as un+1− un
∆t = 1
2(F (tn+1, un+1) + F (tn, un)). (40) Then we have the fully-discrete version for Eqs. (14) as
ivjn+1− vnj
∆t = −ρ 2
∂F (xj, tn+1)
∂x +∂F (xj, tn)
∂x
, j = 0, 1, . . . , N, (41a)
vj0 = f (xj), j = 0, 1, . . . , N. (41b)
To analyze the stability of Eqs. (41), we consider the homogeneous boundary conditions.
For convenience, we denote (vjn+1+ vjn)/2 by vn+1/2j . Multiplying −iρ−1ωj
vn+1/2j ∗
and iρ−1ωjvn+1/2j to Eq. (41a) and its complex conjugate, respectively, summing the resultants, and following a similar approach in semi-discrete scheme, we obtain
1 2ρ∆t
N
X
j=0
ωj(|vn+1j |2− |vnj|2) = ir∗−A−r−+ ir∗+A+r+,
where A± are given in Eq. (25) and
r− = r(v0n+1/2, −∂vn+1/20 /∂x), r+ = r(vn+1/2N , ∂vNn+1/2/∂x).
For τ± given in Eq. (26), A± are zero matrices. Thus, we have
N
X
j=0
ωj|vjn+1|2 =
N
X
j=0
ωj|vjn|2 = · · · =
N
X
j=0
ωj|f (xj)|2,
indicating the stability.
In addition, we have the fully-discrete scheme for Eq. (38a) with vnbeing the numerical solution matrix at time tn
vn+1− vn
∆t = 1
2(Lvn+1+ vn+1R + G(tn+1) + Lvn+ vnR + G(tn)). (42) The Eq. (42) can be rewritten as
Avn+1+ vn+1B = F , (43)
where
A = 1
2I(M )− ∆t
2 L, B = 1
2I(N )− ∆t 2 R, F = vn+ ∆t
2 (Lvn+ vnR + +G(tn) + G(tn+1)) .
with I(M ) ∈ R(M +1)×(M +1) and I(N ) ∈ R(N +1)×(N +1) being identity matrices. Assuming that A and B are diagonalizable, we can directly solve vn+1 by the eigenvector decom-position method [9].
For Eqs. (39), the IMEX-RK method splits the flux function F into two parts as F (t, u(t)) = F[im](t, u(t)) + F[ex](t, u(t)),
where F[im]and F[ex] are flux functions to be treated in implicit and explicit ways, respec-tively. For a s-stage IMEX-RK method, we solve Eqs. (39) numerically by the following steps:
u(i) = un+ ∆t
s
X
j=1
a[im]ij F[im](tn+ cj∆t, u(j)) + ∆t
s
X
j=1
a[ex]ij F[ex](tn+ cj∆t, u(j)), 1 ≤ i ≤ s
un+1= un+ ∆t
s
X
i=1
b[im]i F[im](tn+ ci∆t, u(i)) + ∆t
s
X
i=1
b[ex]i F[ex](tn+ ci∆t, u(i)).
The coefficients of a[im]ij , a[ex]ij , b[im]i , b[ex]i , cj can be found in [3]. We consider Lv+vR+G(t) in Eq. (38a) to be the implicit part, and the explicit part is regarded as zero. Then we have a fully-discrete scheme with IMEX-RK
v(i) = vn+ ∆t
s
X
j=1
a[im]ij (Lv(j)+ v(j)R + G(tn+ cj∆t)), 1 ≤ i ≤ s (44a)
vn+1= vn+ ∆t
s
Xb[im]i (Lv(i)+ v(i)R + G(tn+ ci∆t)). (44b)
We can rewrite Eq. (44a) as
Av(i)+ v(i)B = F(i), 1 ≤ i ≤ s, (45)
where
A = 1
2I(M )−∆t
2 L, B = 1
2I(N )− ∆t 2 R,
F(i) =
vn+ ∆t a[im]11 G(tn+ c1∆t), if i = 1, vn+ ∆t
i−1
X
j=1
a[im]ij F(j)+ ∆t a[im]ii G(tn+ ci∆t), if i 6= 1.
Thus, we solve Eqs. (45) directly to obtain v(i), and then vn+1 can be computed by Eqs. (44b).
3 Numerical results
In this section, we illustrate the proposed methods by several examples. Denote ∆x = 2N−1 which is the mean distance for a set of N + 1 grid points. The time step ∆t is computed adaptively as
∆t = CF L∆x,
where CF L is referred to the Courant-Friedrichs-Lewy number. The convergence order is calculated as
q = log(e(N1)/e(N2)) log(N2/N1) ,
where e(N ) = ku − uNk∞is the maximum error with u and uN being the analytic solution and the numerical solution corresponding to polynomial degree N , respectively.
3.1 One-dimensional problem
Example 1. Let I = [−1, 1]. Consider q(x, t) = ei(x−t) satisfying the problem:
i∂q
∂t = −∂2q
∂x2, x ∈ I, t ≥ 0,
q(x, 0) = eix, x ∈ I,
B±q(±1, t) = α±ei(±1−t)± iβ±ei(±1−t), t ≥ 0,
Table 1-3 shows the results subject to different types of boundary conditions imposed at x = ±1 applying Crank-Nicolson method. For each terminal time T , the error decreases as N increases and the rate of convergence is approximately a two. Compare the results from different collocation points, there is little difference between them and the rate of convergence is similar. Thus, we know that the scheme can achieve a convergent result.
Table 1: The rate of convergence for Example 1 with CGL and LGL points at different terminal time.
The numerical results are computed by Crank-Nicolson method. α± = 1, β±= 0, CF L = 0.1
N
T = 10 T = 100
CGL LGL CGL LGL
e(N ) q e(N ) q e(N ) q e(N ) q
8 7.10e-05 - 7.10e-05 - 7.22e-05 - 7.22e-05 -12 3.12e-05 2.03 3.12e-05 2.03 3.21e-05 2.00 3.21e-05 2.00 16 1.76e-05 1.99 1.76e-05 1.99 1.81e-05 1.98 1.81e-05 1.98 20 1.13e-05 1.98 1.13e-05 1.98 1.18e-05 1.94 1.18e-05 1.94
Table 2: The rate of convergence for Example 1 with CGL and LGL points at different terminal time.
The numerical results are computed by Crank-Nicolson method. α± = 0, β±= 1, CF L = 0.1
N
T = 10 T = 100
CGL LGL CGL LGL
e(N ) q e(N ) q e(N ) q e(N ) q
8 1.05e-04 - 1.05e-04 - 6.53e-05 - 6.53e-05 -12 4.64e-05 2.02 4.64e-05 2.02 2.86e-05 2.03 2.86e-05 2.03 16 2.61e-05 2.00 2.61e-05 2.00 1.64e-05 1.94 1.64e-05 1.94 20 1.67e-05 2.00 1.67e-05 2.00 1.06e-05 1.97 1.06e-05 1.97
The IMEX-RK method can be applied to solve Example 1, and the results are shown in Table 4-6. We see the convergence order is approximately a three in each case. So the scheme produced by IMEX-RK is still convergent whether we use CGL or LGL points.
Fig. 1 reveals that the discrete energy is conserved, where ∆E is defined as
∆E =
N
X
j=1
wj|vj|2−
N
X
j=1
wj|f (xj)|2.
Table 3: The rate of convergence for Example 1 with CGL and LGL points at different terminal time.
The numerical results are computed by Crank-Nicolson method. α± = 1.5, β± = 0.5, CF L = 0.1
N
T = 10 T = 100
CGL LGL CGL LGL
e(N ) q e(N ) q e(N ) q e(N ) q
8 2.24e-04 - 2.24e-04 - 2.05e-04 - 2.05e-04 -12 9.97e-05 2.00 9.97e-05 2.00 9.07e-05 2.01 9.07e-05 2.01 16 5.60e-05 2.00 5.60e-05 2.00 5.10e-05 2.00 5.10e-05 2.00 20 3.59e-05 2.00 3.59e-05 2.00 3.26e-05 2.01 3.26e-05 2.01
Table 4: The rate of convergence for Example 1 with CGL and LGL points at different terminal time.
The numerical results are computed by third-order IMEX-RK. α±= 1, β±= 0, CF L = 0.1
N
T = 10 T = 100
CGL LGL CGL LGL
e(N ) q e(N ) q e(N ) q e(N ) q
8 1.43e-06 - 1.43e-06 - 1.33e-06 - 1.34e-06 -12 4.05e-07 3.10 4.07e-07 3.10 4.13e-07 2.88 4.14e-07 2.91 16 1.72e-07 2.98 1.72e-07 2.99 1.75e-07 2.98 1.75e-07 2.99 20 8.84e-08 2.98 8.84e-08 2.98 8.86e-08 3.05 8.86e-08 3.05
Table 5: The rate of convergence for Example 1 with CGL and LGL points at different terminal time.
The numerical results are computed by third-order IMEX-RK. α±= 0, β±= 1, CF L = 0.1
N
T = 10 T = 100
CGL LGL CGL LGL
e(N ) q e(N ) q e(N ) q e(N ) q
8 1.81e-06 - 1.81e-06 - 1.16e-06 - 1.16e-06 -12 2.68e-07 4.71 2.68e-07 4.71 2.79e-07 3.50 2.79e-07 3.50 16 1.14e-07 2.98 1.14e-07 2.98 1.20e-07 2.95 1.20e-07 2.95 20 5.87e-08 2.97 5.87e-08 2.97 6.18e-08 2.96 6.18e-08 2.96
Table 6: The rate of convergence for Example 1 with CGL and LGL points at different terminal time.
The numerical results are computed by third-order IMEX-RK. α±= 1.5, β± = 0.5, CF L = 0.1
N
T = 10 T = 100
CGL LGL CGL LGL
e(N ) q e(N ) q e(N ) q e(N ) q
8 2.42e-06 - 2.45e-06 - 2.17e-06 - 2.17e-06 -12 7.21e-07 2.99 7.21e-07 3.02 6.66e-07 2.92 6.66e-07 2.92 16 3.04e-07 3.00 3.04e-07 3.00 2.82e-07 2.99 2.82e-07 2.99 20 1.56e-07 3.00 1.56e-07 3.00 1.44e-07 3.00 1.44e-07 3.00
0 2 4 6 8 10 0
0.5 1 1.5 2 2.5 3 3.5x 10−5
T
∆ E
0 2 4 6 8 10
−1.5
−1
−0.5 0 0.5 1 1.5x 10−7
T
∆ E
Figure 1: The difference of discrete energy for Example 1 subject to Dirichlet boundary condition. Left:
Time integration by Crank-Nicolson method. Left: Time integration by third-order IMEX-RK.
Example 2. Let I = [−1, 1]. Consider q(x, t) = sin(πx) sin(πt) satisfying the problem:
i∂q
∂t = −∂2q
∂x2 + iπ sin(πx) cos(πt) − π2sin(πx) sin(πt), x ∈ I, t ≥ 0,
q(x, 0) = 0, x ∈ I,
B±q(±1, t) = α±sin(±π) sin(πt) ± β±π cos(±π) sin(πt), t ≥ 0,
In this experiment, the partial differential equation contains a source term. To solve Example 2 with the IMEX-RK method, the source term is treated in the explicit flux. The results are shown in Table 7- 8. We observe an order reduction occurs when the boundary conditions are not Dirichlet-type. However, if we apply the Crank-Nicolson method to solve the same problem, the convergence order is approximately a two for every types of boundary conditions. Since the Dirichlet boundary condition in this problem is exactly zero, we guess that the order reduction occurs if we use IMEX-RK in time to solve a problem which contains a source term and nonzero boundary conditions.
Table 7: Convergence order for Example 2. The numerical solutions are computed by third-order IMEX-RK and Crank-Nicolson at T = 1 with CGL points. CF L = 0.1
N
α±= 1, β± = 0 α±= 0, β± = 1
IMEX-RK CN IMEX-RK CN
e(N ) q e(N ) q e(N ) q e(N ) q
8 1.72e-04 - 8.14e-05 - 1.01e-03 - 5.98e-04 -12 5.38e-05 2.87 3.40e-05 2.16 4.02e-04 2.27 2.65e-04 2.01 16 2.24e-05 3.04 1.92e-05 1.98 1.94e-04 2.54 1.48e-04 2.02 20 1.15e-05 2.97 1.25e-05 1.93 1.11e-04 2.48 9.41e-05 2.03
Table 8: Convergence order for Example 2. The numerical solutions are computed by third-order IMEX-RK and Crank-Nicolson at T = 1 with LGL points. CF L = 0.1
N
α±= 1, β± = 0 α±= 0, β± = 1
IMEX-RK CN IMEX-RK CN
e(N ) q e(N ) q e(N ) q e(N ) q
8 1.77e-04 - 7.85e-05 - 1.00e-03 - 6.08e-04 -12 5.37e-05 2.94 3.39e-05 2.07 4.02e-04 2.25 2.64e-04 2.05 16 2.26e-05 3.01 1.94e-05 1.95 1.94e-04 2.54 1.48e-04 2.02 20 1.15e-05 3.03 1.24e-05 1.99 1.11e-04 2.48 9.41e-05 2.03
Example 3. Let I = [−1, 1]. Consider q(x, t) = ei(x+t) satisfying the problem:
i∂q
∂t = −∂2q
∂x2 − 2ei(x+t), x ∈ I, t ≥ 0,
q(x, 0) = eix, x ∈ I,
B±q(±1, t) = α±ei(±1+t)± iβ±ei(±1+t), t ≥ 0,
We construct Example 3 whose boundary conditions are all nonhomogeneous. Table 9-10 show the results as we solve this problem with IMEX-RK and Crank-Nicolson method, respectively. The order reduction still occurs in the case of IMEX-RK. Meanwhile, the results produced by Crank-Nicolson method present a second-order convergence as we expect.
Table 9: Convergence order for Example 3. The numerical results are computed by third-order IMEX-RK and Crank-Nicolson with CGL points at T = 1. CF L = 0.1
N
α±= 1, β± = 0 α±= 0, β± = 1
IMEX-RK CN IMEX-RK CN
e(N ) q e(N ) q e(N ) q e(N ) q
8 3.25e-04 - 3.33e-05 - 3.27e-05 - 6.37e-05 -12 1.56e-04 1.81 1.48e-05 1.99 1.20e-05 2.46 2.83e-05 2.00 16 8.28e-05 2.20 8.48e-06 1.95 5.47e-06 2.74 1.59e-05 2.00 20 5.16e-05 2.12 5.46e-06 1.97 3.05e-06 2.62 1.02e-05 2.00
Table 10: Convergence order for Example 3. The numerical results are computed by third-order IMEX-RK and Crank-Nicolson with LGL points at T = 1. CF L = 0.1
N
α±= 1, β± = 0 α±= 0, β± = 1
IMEX-RK CN IMEX-RK CN
e(N ) q e(N ) q e(N ) q e(N ) q
8 3.25e-04 - 3.33e-05 - 3.27e-05 - 6.35e-05 -12 1.56e-04 1.81 1.48e-05 1.99 1.20e-05 2.46 2.83e-05 2.00 16 8.28e-05 2.20 8.49e-06 1.94 5.47e-06 2.74 1.59e-05 2.00 20 5.16e-05 2.12 5.46e-06 1.98 3.05e-06 2.62 1.02e-05 2.00
3.2 Two-dimensional problem
Example 4. Let I2 = [−1, 1]×[−1, 1]. Consider q(x, y, t) = e−2iπ2tcos(π(x+y)) satisfying the problem:
i∂q
∂t = −∂2q
∂x2 − ∂2q
∂y2, (x, y) ∈ I2, t ≥ 0,
q(x, y, 0) = cos(π(x + y)), (x, y) ∈ I2,
B(a)q(−1, y, t) = α(a)e−2iπ2tcos(π(y − 1)) + β(a)πe−2iπ2tsin(π(y − 1)), y ∈ I, t ≥ 0, B(b)q(+1, y, t) = α(b)e−2iπ2tcos(π(y + 1)) − β(b)πe−2iπ2tsin(π(y + 1)), y ∈ I, t ≥ 0, B(c)q(x, −1, t) = α(c)e−2iπ2tcos(π(x − 1)) + β(c)πe−2iπ2tsin(π(x − 1)), x ∈ I, t ≥ 0, B(d)q(x, +1, t) = α(d)e−2iπ2tcos(π(x + 1)) − β(d)πe−2iπ2tsin(π(x + 1)), x ∈ I, t ≥ 0.
We consider Example 4 to test the scheme for two-dimension problems. For simpli-fication, the grid revolutions N and M are set to be equal. We apply Crank-Nicolson method in time to solve this problem, and show the convergence results in Table 11-13.
We use different collocation points to compute the numerical solutions. For each terminal time T , The error decreases as N increases with a second-order convergence. The scheme is stable after a long-time computation (T = 100). Moveover, for a fixed N , the error grows linearly in time.
Table 14-16 present the results computed by IMEX-RK in time. For each terminal time T , the the rate of convergence decays rapidly in the beginning, and then goes down to third order. We see the error grows linearly for a fixed grid revolution N . It shows
Table 11: The rate of convergence for Example 4 with CGL and LGL points at different terminal time.
The numerical results are computed by Crank-Nicolson method. For γ ∈ {a, b, c, d}, α(γ) = 1, β(γ) = 0.
CF L = 0.01
N
T = 10 T = 100
CGL LGL CGL LGL
e(N ) q e(N ) q e(N ) q e(N ) q
8 3.58e-02 - 3.37e-02 - 3.50e-01 - 3.32e-01 -12 1.79e-02 1.71 1.79e-02 1.57 1.78e-01 1.67 1.77e-01 1.55 16 9.79e-03 2.10 9.92e-03 2.04 9.71e-02 2.10 9.84e-02 2.05 20 6.30e-03 1.97 6.24e-03 2.08 6.28e-02 1.96 6.21e-02 2.06
Table 12: The rate of convergence for Example 4 with CGL and LGL points at different terminal time.
The numerical results are computed by third-order Crank-Nicolson method. For γ ∈ {a, b, c, d}, α(γ) = 0, β(γ)= 1. CF L = 0.01
N
T = 10 T = 100
CGL LGL CGL LGL
e(N ) q e(N ) q e(N ) q e(N ) q
8 4.81e-02 - 4.84e-02 - 3.95e-01 - 3.95e-01 -12 1.80e-02 2.42 1.80e-02 2.43 1.78e-01 1.97 1.78e-01 1.97 16 1.02e-02 2.00 1.02e-02 2.00 1.00e-01 2.00 1.00e-01 2.00 20 6.50e-03 2.00 6.50e-03 2.00 6.41e-02 2.00 6.41e-02 2.00
Table 13: The rate of convergence for Example 4 with CGL and LGL points at different terminal time.
The numerical results are computed by third-order Crank-Nicolson method. For γ ∈ {a, b, c, d}, α(γ) = 1.5, β(γ)= 0.5. CF L = 0.01
N
T = 10 T = 100
CGL LGL CGL LGL
e(N ) q e(N ) q e(N ) q e(N ) q
8 1.79e-02 - 1.80e-02 - 1.68e-02 - 1.69e-02 -12 4.76e-04 8.95 4.92e-04 8.88 6.58e-04 7.99 6.75e-04 7.94 16 2.76e-04 1.88 2.80e-04 1.95 3.88e-04 1.84 3.91e-04 1.90 20 1.79e-04 1.94 1.80e-04 1.98 2.46e-04 2.04 2.46e-04 2.07
Table 14: The rate of convergence for Example 4 with CGL and LGL points at different terminal time.
The numerical results are computed by third-order IMEX-RK. For γ ∈ {a, b, c, d}, α(γ) = 1, β(γ) = 0.
CF L = 0.01
N
T = 10 T = 100
CGL LGL CGL LGL
e(N ) q e(N ) q e(N ) q e(N ) q
8 4.58e-03 - 4.82e-03 - 2.26e-02 - 2.33e-02 -12 1.84e-04 7.93 1.83e-04 8.06 1.82e-03 6.21 1.82e-03 6.29 16 7.52e-05 3.11 7.61e-05 3.05 7.44e-04 3.11 7.55e-04 3.05 20 3.87e-05 2.98 3.83e-05 3.08 3.85e-04 2.96 3.81e-04 3.06
Table 15: The rate of convergence for Example 4 with CGL and LGL points at different terminal time.
The numerical results are computed by third-order IMEX-RK. For γ ∈ {a, b, c, d}, α(γ) = 0, β(γ) = 1.
CF L = 0.01
N
T = 10 T = 100
CGL LGL CGL LGL
e(N ) q e(N ) q e(N ) q e(N ) q
8 1.05e-02 - 1.06e-02 - 9.26e-03 - 9.18e-03 -12 1.85e-04 9.95 1.84e-04 10.00 1.82e-03 4.02 1.82e-03 3.99 16 7.84e-05 2.99 7.84e-05 2.96 7.68e-04 3.00 7.68e-04 3.00 20 4.02e-05 3.00 4.02e-05 3.00 3.93e-04 3.00 3.93e-04 3.00
Table 16: The rate of convergence for Example 4 with CGL and LGL points at different terminal time.
The numerical results are computed by third-order IMEX-RK. For γ ∈ {a, b, c, d}, α(γ)= 1.5, β(γ)= 0.5.
CF L = 0.01
N
T = 10 T = 100
CGL LGL CGL LGL
e(N ) q e(N ) q e(N ) q e(N ) q
8 1.62e-02 - 1.62e-02 - 1.41e-02 - 1.41e-02 -12 1.21e-05 17.74 1.22e-05 17.74 8.97e-06 18.15 9.13e-06 18.11 16 2.33e-06 5.73 2.34e-06 5.73 3.25e-06 3.53 3.25e-06 3.59 20 1.20e-06 2.97 1.20e-06 3.00 1.66e-06 3.01 1.65e-06 3.04
Example 5. Let I2 = [−1, 1] × [−1, 1]. Consider q(x, y, t) = e−itsin(πx) sin(πy) satisfying the problem:
i∂q
∂t = −∂2q
∂x2 − ∂2q
∂y2 + (1 − 2π2)e−itsin(πx) sin(πy), (x, y) ∈ I2, t ≥ 0,
q(x, y, 0) = sin(πx) sin(πy), (x, y) ∈ I2,
B(a)q(−1, y, t) = α(a)e−itsin(−π) sin(πy) − β(a)πe−itcos(−π) sin(πy), y ∈ I, t ≥ 0, B(b)q(+1, y, t) = α(b)e−itsin(+π) sin(πy) + β(b)πe−itcos(+π) sin(πy), y ∈ I, t ≥ 0, B(c)q(x, −1, t) = α(c)e−itsin(πx) sin(−π) − β(c)πe−itsin(πx) cos(−π), x ∈ I, t ≥ 0, B(d)q(x, +1, t) = α(d)e−itsin(πx) sin(+π) + β(d)πe−itsin(πx) cos(+π), x ∈ I, t ≥ 0.
We can use the scheme to solve Example 5 which contains a source term. Table 17-18 present the results subject to different types of boundary conditions. We see the results computed by the Crank-Nicolson method achieve a convergence of second order indeed.
The order reduction still occurs while we use IMEX-RK to solve this problem with the nonhomogeneous boundary conditions.
It has been observed that if we use pseudospectral penalty method to construct a scheme, the order reduction may occurs for Runge-Kutta methods. The paper [2] provide us a procedure to deal with explicit Runge-Kutta. However, there is no way for implicit Runge-Kutta up to now.
Table 17: Convergence order for Example 5 at T = 1. The numerical solutions are computed by third-order IMEX-RK and Crank-Nicolson with LGL points. CF L = 0.1
N
α(γ)= 1, β(γ)= 0 α(γ) = 0, β(γ) = 1
IMEX-RK CN IMEX-RK CN
e(N ) q e(N ) q e(N ) q e(N ) q
8 1.12e-03 - 1.07e-03 - 1.68e-03 - 3.20e-03 -12 1.69e-05 10.36 4.59e-07 19.13 2.51e-04 4.69 6.98e-06 15.11 16 6.96e-06 3.07 1.43e-07 4.05 1.24e-04 2.44 3.76e-06 2.15 20 3.51e-06 3.07 7.50e-08 2.89 6.98e-05 2.59 2.45e-06 1.92
Table 18: Convergence order for Example 5 at T = 1. The numerical solutions are computed by third-order IMEX-RK and Crank-Nicolson with CGL points. CF L = 0.1
N
α(γ)= 1, β(γ)= 0 α(γ) = 0, β(γ) = 1
IMEX-RK CN IMEX-RK CN
e(N ) q e(N ) q e(N ) q e(N ) q
8 1.13e-03 - 1.27e-03 - 1.62e-03 - 3.17e-03 -12 1.69e-05 10.36 4.66e-07 19.52 2.55e-04 4.56 6.99e-06 15.08 16 6.86e-06 3.13 1.41e-07 4.15 1.25e-04 2.47 3.76e-06 2.15 20 3.54e-06 2.96 7.57e-08 2.79 6.99e-05 2.61 2.45e-06 1.92
4 Concluding remarks
In this study, we proposed a numerical scheme for solving the Sch¨ordinger equation based on pseudospectral penalty method. For stable computations, we determine the penalty parameters subject to different types of boundary conditions by conducting the energy estimate. Although we establish the stability through LGL points, the scheme still works when we employ CGL points. We apply Crank-Nicolson method and IMEX-RK for time-discretization. Several numerical experiments are shown to validate the scheme, and we observe the expected convergence rate in most cases. An order reduction occurs only when we use IMEX-RK to solve the problems that contain a source term and nonhomogeneous boundary conditions.
The present method well solves the problems in a simple domain with scalar boundary parameters. In the future, we can conduct some numerical experiments for the problems whose boundary conditions are parameterized by functions. The method may be gener-alized to solve the non-linear Sch¨ordinger equations defined on complicated domains.
References
[1] E. Schr¨odinger, An Undulatory Theory of the Mechanics of Atoms and Molecules.
Phys. Rev. Vol. 28, No. 6 (1926) pp. 1049-1070
[2] W.S. Don, D. Gottlieb, The Chebyshev-Legendre method: Implementing Legendre methods on Chebyshev points. SIAM J. Numer. Anal., 31 (1994), pp. 1519-1534.
[3] C.A. Kennedy, M.H. Carpenter, Additive Runge-Kutta schemes for convection-diffusion-reaction equations. Appl. Numer. Math., 44 (2003), pp. 139-181.
[4] P. D. Lax, R. D. Richtmyer, Survey of the stability of linear finite difference equations.
Communications on Pure and Applied Mathematics, 9 (1956) 267-293.
[5] J.S. Hesthaven, S. Gottlieb, D. Gottlieb, Spectral Methods for Time-Dependent Problems. Cambridge Universuty Press. Prentice Hall, (2007).
[6] Anna Nissen, Gunilla Kreiss, Margot Gerritsen, High Order Stable Finite Difference Methods for the Schr¨odinger Equation. J. Sci. Comput. 55(1): 173-199 (2013) [7] J. Crank, P. Nicolson, A practical method for numerical evaluation of solutions of
partial differential equations of the heat conduction type. Proc. Camb. Phil. Soc., 43, 1 (1947) 50-67.
[8] A. Kanevsky, M.H. Carpenter, D. Gottlieb, J.S. Hesthaven, Application of implic-itVexplicit high order RungeVKutta methods to discontinuous-Galerkin schemes. J.
Comput. Phys., 225 (2007), pp. 1753-1781.
[9] Robert E. Lynch, John R. Rice, Donald H. Thomas, Direct solution of partial differ-ential equations by tensor product methods. Numer. Math., 6 (1964) 185-199.
[10] D. Funaro, D. Gottlieb, A new method of imposing boundary conditions in pseu-dospectral approximations of hyperbolic equations. Math. Comput. 51 (1988) 599-613.
[11] Delfour, M, Fortin, M, Payr, G, Finite-difference solutions of a non-linear Schr¨odinger equation. Journal of Computational Physics, 1981, Vol.44(2), pp.277-288