# Symplectic Schemes for Linear Stochastic Schrödinger Equations with Variable Coefficients

Authors: Xiuling Yin; Yanqin Liu
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

This paper proposes a kind of symplectic schemes for linear Schrödinger equations with
variable coefficients and a stochastic perturbation term by using compact schemes in space. The
numerical stability property of the schemes is analyzed. The schemes preserve a discrete charge conservation
law. They also follow a discrete energy transforming formula. The numerical experiments
verify our analysis.


---

## Body

## 1. Introduction

Differential equations (DEs) are important models in sciences and engineering. By theoretical and numerical analysis of DEs, we can yield some mathematical explanation of many phenomena in applied sciences [1–5]. Time-dependent Schrödinger equations arise in quantum physics, optics, and many other fields [6, 7]. Some numerical methods for such equations, such as symplectic scheme and multisymplectic schemes, have been proposed in [8–14]. The schemes possess good numerical stability. Compact schemes are popular recently due to high accuracy, compactness, and economic resource in scientific computation [15–17]. In this paper, applying compact operators, we construct symplectic methods to the initial boundary problems of the linear Schrödinger equation with a variable coefficient and a stochastic perturbation term (denoted by LSES):
(1)iut+ux4+f(x)u=ϵu∘χ˙,x∈[0,L],u(x,0)=u0(x),t∈[0,T],u(0,t)=u(L,t),
where i2=-1, f(x) is a real differential function, u0(x) is a differential function, ϵ is a small real number, and ∘ means Stratonovich product. χ˙ is a real-valued white noise which is delta correlated in time and either smooth or delta correlated in space. For an integer m, uxm and utm mean the m-order partial derivatives of u with respect to x and t, respectively. The system (1) with ϵ=0 is a deterministic system. When ϵ is small, we can think that (1) is perturbed by the stochastic term.By multiplying (1) by u- or u-t and then integrating it with respect to t and x, it is easy to verify the following result.Proposition 1.
Under the periodic boundary condition,
(a)
the solution of (1) satisfies the charge conservation law
(2)
𝒬
(
t
)
=
∫
0
L
|
u
(
x
,
t
)
|
2
d
x
=
𝒬
(
0
)
,
(b)
the corresponding deterministic system (
ϵ
=
0
) of (1) possesses the energy conservation law
(3)
ℰ
(
t
)
=
∫
0
L
‍
|
u
x
x
|
2
+
f
(
x
)
|
u
|
2
d
x
=
ℰ
(
0
)
.The paper is organized as follows. In Section2, we give a symplectic structure of the LSES. In Section 3, we present the new symplectic methods to the LSES. First, we use a kind of compact schemes in discretization of spatial derivative. Then, in temporal discretization, we adopt the symplectic midpoint method. The new methods are denoted by LSC schemes. We also analyze the numerical stability of LSC schemes. We give two numerical examples to support our theory in Section 4. At last, we make some conclusions.


## 2. Symplectic Structure of the LSES

Letu=p+iq. The LSES (1) can be written in
(4)pt+qx4+f(x)q=ϵq∘χ˙,-qt+px4+f(x)p=ϵp∘χ˙.Introducing the variablez=(p,q)T, (4) reads in stochastic symplectic context
(5)zt=J-1∇zH(z)+ϵJ-1∇zS(z)∘χ˙,
where
(6)H(z)=12(pxx2+qxx2)+f(x)2(p2+q2),S(z)=-12(p2+q2),J=(01-10).
The system satisfies the symplectic conservation law [7, 12, 18]:
(7)ωt=0,ω=dp∧dq.
Numerical methods which preserve the discrete symplectic conservation law are called symplectic methods. Symplectic methods have good numerical stability.


## 3. LSC Schemes

### 3.1. Compact Scheme

Introduce the following uniform mesh grids:(8)xk=kh,k=0,1,…,N;tn=nτ,n=0,1,…,
where h=L/N and τ are spatial and temporal step sizes, respectively. Denote the numerical values of u(xk,tn) at the nodes (xk,tn) by ukn. The symbols un and uk mean the numerical solution vectors at t=tn and x=xk with components ukn, respectively. Furthermore, we denote
(9)ukn+(1/2):=ukn+1+ukn2,δtukn+(1/2):=ukn+1-uknτ.Introducing the following linear operators(10)𝒜uk=αuk-1+uk+αuk+1,ℬuk=buk+3-9uk+1+16uk-9uk-1+uk-36h4+auk+2-4uk+1+6uk-4uk-1+uk-2h4,
we adopt formula [19]
(11)δx4uk=𝒜-1ℬuk
to approximate ux4, which means that
(12)𝒜δx4uk=ℬuk.
By Taylor expansion, we can derive a family of fourth-order schemes with
(13)a=2(1-α),b=4α-1.
The leading term of the truncation error of the method is ((7-26α)/240)(ux8)kh4. If b=0, we get a scheme with smaller stencil. A sixth-order scheme is obtained with
(14)α=726,a=1913,b=113.
Denote two symmetric and cyclic matrices by(15)A=[1α0⋯⋯0αα1α0⋯⋯00α1α0⋯0⋮⋱⋱⋱⋱⋱⋮⋮⋮⋱α1α000⋯0α1αα0⋯⋯0α1]N×N,B=16h4[c0c16ab0⋯⋯0b6ac1c1c0c16ab0⋯⋯0b6a6ac1c0c16ab0⋯⋯0bb6ac1c0c16ab0⋯⋯00b6ac1c0c16ab0⋯0⋮⋱⋱⋱⋱⋱⋱⋱⋱⋱⋮0⋯0b6ac1c0c16ab00⋯⋯0b6ac1c0c16abb0⋯⋯0b6ac1c0c16a6ab0⋯⋯0b6ac1c0c1c16ab0⋯⋯0b6ac1c0]N×N,where c0=16b+36a and c1=-9b-24a. Then the matrix form of (11) is
(16)δx4un=A-1Bun.


### 3.2. Discretization of the LSES

Applying the approximation (11) to linear system (4), we obtain the following semidiscretization stochastic Hamiltonian system:
(17)(zk)t=J-1∇zH-(zk)+ϵJ-1∇zS-(zk)∘χ˙k,
where
(18)H-(z)=f(xk)2(pk2+qk2)+12δx4(pk2+qk2),S-(z)=-12(pk2+qk2).
In temporal discretization of (17), we apply the symplectic midpoint method
(19)δtzkn+(1/2)=J-1∇zH-(zkn+(1/2))+ϵJ-1∇zH-(zkn+(1/2))∘χ˙kn+(1/2).
Its componentwise formulation is
(20)pkn+1-pknτ=-f(xk)qkn+(1/2)-δx4qkn+(1/2)+ϵqkn+(1/2)∘χ˙kn+(1/2),qkn+1-qknτ=f(xk)pkn+(1/2)+δx4pkn+(1/2)-ϵpkn+(1/2)∘χ˙kn+(1/2).According to the Fourier analysis, the LSC schemes (19) are unconditionally stable. In fact, we can derive
(21)δx4zkn=μzkn,
where
(22)μ=4(ω-1)2(3a+2b+bω)3h4(1+2αω),ω=cos⁡⁡βh.
Then, with (19) and (21), we can obtain zn+1=Gzn with (23)(1+c2τ24)G=(1-c2τ24-cτcτ1-c2τ24),
where c=μ+f(xk)-ϵ∘χ˙kn+(1/2). By direct computation, we can derive that the spectral radius of the matrix G is 1 and ∥G∥2=1. Therefore, the scheme (19) is unconditionally stable. Moreover, by symmetry, they are nondissipative.Theorem 2.
Let∥zn∥2=h∑kzknzkn¯. Then, ∥zn∥ is the discrete charge invariant of the LSC schemes (19), which implies the discrete charge conservation law of (2).Scheme (19) can be rewritten in compact form
(24)i𝒜ukn+1-uknτ+ℬukn+(1/2)+𝒜f(xk)ukn+(1/2)=ϵ𝒜ukn+(1/2)∘χ˙kn+(1/2).
Multiplying (24) by (ukn+1¯-ukn¯) and summing over k, we obtain
(25)iτ2∑k‍[ℬukn+1ukn+1¯-ℬuknukn¯]+iτ2∑k‍𝒜f(xk)[ukn+1ukn+1¯-uknukn¯]-iτϵ2∑k‍[𝒜ukn+1∘χ˙kn+(1/2)ukn+1¯-𝒜ukn∘χ˙kn+(1/2)ukn¯]+iτ2∑k‍[ℬuknukn+1¯-ℬukn+1ukn¯]+iτ2∑k𝒜f(xk)[uknukn+1¯-ukn+1ukn¯]‍-iτϵ2∑k‍[𝒜ukn∘χ˙kn+(1/2)ukn+1¯-𝒜ukn+1∘χ˙kn+(1/2)ukn¯]=∑k‍[𝒜(ukn+1-ukn)(ukn+1¯-ukn¯)].
Since 𝒜 and ℬ are symmetric, the first three summation terms in the above equality are purely imaginary, while the last four summation terms are real. Denote
(26)En=h∑k‍[ℬuknukn¯]+h∑k‍[𝒜f(xk)uknukn¯]-h2ϵ∑k‍[𝒜ukn∘χ˙knukn¯].
Now, taking the imaginary parts of (25), we can get that
(27)En+1-En=h2ϵ∑k‍[𝒜ukn+1∘χ˙knukn+1¯-𝒜ukn∘χ˙kn+1ukn¯].
Denote
(28)vkn=b6uk+3n-uknh2,v~kn=3b2uk+1n-uknh2,wkn=auk+2n-uknh2,w~kn=2auk+1n-uknh2,ykn=𝒜f(xk)ukn,y~kn=𝒜ukn∘χ˙kn,〈un,yn〉=h∑k‍uknykn¯,∥un∥2=〈un,un〉.
According to the Green formula, we obtain that
(29)h∑k‍[ℬuknukn¯]=∥v~n∥2-∥vn∥2+∥w~n∥2-∥wn∥2.
Then,
(30)En=∥v~n∥2-∥vn∥2+∥w~n∥2-∥wn∥2+〈yn,un〉-ϵ2〈y~n,un〉.
Therefore, from the above analysis, we give the following result.Theorem 3.
Under the periodic boundary condition, the LSC schemes (19) satisfy the discrete energy transforming law (27).


## 3.1. Compact Scheme

Introduce the following uniform mesh grids:(8)xk=kh,k=0,1,…,N;tn=nτ,n=0,1,…,
where h=L/N and τ are spatial and temporal step sizes, respectively. Denote the numerical values of u(xk,tn) at the nodes (xk,tn) by ukn. The symbols un and uk mean the numerical solution vectors at t=tn and x=xk with components ukn, respectively. Furthermore, we denote
(9)ukn+(1/2):=ukn+1+ukn2,δtukn+(1/2):=ukn+1-uknτ.Introducing the following linear operators(10)𝒜uk=αuk-1+uk+αuk+1,ℬuk=buk+3-9uk+1+16uk-9uk-1+uk-36h4+auk+2-4uk+1+6uk-4uk-1+uk-2h4,
we adopt formula [19]
(11)δx4uk=𝒜-1ℬuk
to approximate ux4, which means that
(12)𝒜δx4uk=ℬuk.
By Taylor expansion, we can derive a family of fourth-order schemes with
(13)a=2(1-α),b=4α-1.
The leading term of the truncation error of the method is ((7-26α)/240)(ux8)kh4. If b=0, we get a scheme with smaller stencil. A sixth-order scheme is obtained with
(14)α=726,a=1913,b=113.
Denote two symmetric and cyclic matrices by(15)A=[1α0⋯⋯0αα1α0⋯⋯00α1α0⋯0⋮⋱⋱⋱⋱⋱⋮⋮⋮⋱α1α000⋯0α1αα0⋯⋯0α1]N×N,B=16h4[c0c16ab0⋯⋯0b6ac1c1c0c16ab0⋯⋯0b6a6ac1c0c16ab0⋯⋯0bb6ac1c0c16ab0⋯⋯00b6ac1c0c16ab0⋯0⋮⋱⋱⋱⋱⋱⋱⋱⋱⋱⋮0⋯0b6ac1c0c16ab00⋯⋯0b6ac1c0c16abb0⋯⋯0b6ac1c0c16a6ab0⋯⋯0b6ac1c0c1c16ab0⋯⋯0b6ac1c0]N×N,where c0=16b+36a and c1=-9b-24a. Then the matrix form of (11) is
(16)δx4un=A-1Bun.


## 3.2. Discretization of the LSES

Applying the approximation (11) to linear system (4), we obtain the following semidiscretization stochastic Hamiltonian system:
(17)(zk)t=J-1∇zH-(zk)+ϵJ-1∇zS-(zk)∘χ˙k,
where
(18)H-(z)=f(xk)2(pk2+qk2)+12δx4(pk2+qk2),S-(z)=-12(pk2+qk2).
In temporal discretization of (17), we apply the symplectic midpoint method
(19)δtzkn+(1/2)=J-1∇zH-(zkn+(1/2))+ϵJ-1∇zH-(zkn+(1/2))∘χ˙kn+(1/2).
Its componentwise formulation is
(20)pkn+1-pknτ=-f(xk)qkn+(1/2)-δx4qkn+(1/2)+ϵqkn+(1/2)∘χ˙kn+(1/2),qkn+1-qknτ=f(xk)pkn+(1/2)+δx4pkn+(1/2)-ϵpkn+(1/2)∘χ˙kn+(1/2).According to the Fourier analysis, the LSC schemes (19) are unconditionally stable. In fact, we can derive
(21)δx4zkn=μzkn,
where
(22)μ=4(ω-1)2(3a+2b+bω)3h4(1+2αω),ω=cos⁡⁡βh.
Then, with (19) and (21), we can obtain zn+1=Gzn with (23)(1+c2τ24)G=(1-c2τ24-cτcτ1-c2τ24),
where c=μ+f(xk)-ϵ∘χ˙kn+(1/2). By direct computation, we can derive that the spectral radius of the matrix G is 1 and ∥G∥2=1. Therefore, the scheme (19) is unconditionally stable. Moreover, by symmetry, they are nondissipative.Theorem 2.
Let∥zn∥2=h∑kzknzkn¯. Then, ∥zn∥ is the discrete charge invariant of the LSC schemes (19), which implies the discrete charge conservation law of (2).Scheme (19) can be rewritten in compact form
(24)i𝒜ukn+1-uknτ+ℬukn+(1/2)+𝒜f(xk)ukn+(1/2)=ϵ𝒜ukn+(1/2)∘χ˙kn+(1/2).
Multiplying (24) by (ukn+1¯-ukn¯) and summing over k, we obtain
(25)iτ2∑k‍[ℬukn+1ukn+1¯-ℬuknukn¯]+iτ2∑k‍𝒜f(xk)[ukn+1ukn+1¯-uknukn¯]-iτϵ2∑k‍[𝒜ukn+1∘χ˙kn+(1/2)ukn+1¯-𝒜ukn∘χ˙kn+(1/2)ukn¯]+iτ2∑k‍[ℬuknukn+1¯-ℬukn+1ukn¯]+iτ2∑k𝒜f(xk)[uknukn+1¯-ukn+1ukn¯]‍-iτϵ2∑k‍[𝒜ukn∘χ˙kn+(1/2)ukn+1¯-𝒜ukn+1∘χ˙kn+(1/2)ukn¯]=∑k‍[𝒜(ukn+1-ukn)(ukn+1¯-ukn¯)].
Since 𝒜 and ℬ are symmetric, the first three summation terms in the above equality are purely imaginary, while the last four summation terms are real. Denote
(26)En=h∑k‍[ℬuknukn¯]+h∑k‍[𝒜f(xk)uknukn¯]-h2ϵ∑k‍[𝒜ukn∘χ˙knukn¯].
Now, taking the imaginary parts of (25), we can get that
(27)En+1-En=h2ϵ∑k‍[𝒜ukn+1∘χ˙knukn+1¯-𝒜ukn∘χ˙kn+1ukn¯].
Denote
(28)vkn=b6uk+3n-uknh2,v~kn=3b2uk+1n-uknh2,wkn=auk+2n-uknh2,w~kn=2auk+1n-uknh2,ykn=𝒜f(xk)ukn,y~kn=𝒜ukn∘χ˙kn,〈un,yn〉=h∑k‍uknykn¯,∥un∥2=〈un,un〉.
According to the Green formula, we obtain that
(29)h∑k‍[ℬuknukn¯]=∥v~n∥2-∥vn∥2+∥w~n∥2-∥wn∥2.
Then,
(30)En=∥v~n∥2-∥vn∥2+∥w~n∥2-∥wn∥2+〈yn,un〉-ϵ2〈y~n,un〉.
Therefore, from the above analysis, we give the following result.Theorem 3.
Under the periodic boundary condition, the LSC schemes (19) satisfy the discrete energy transforming law (27).


## 4. Numerical Examples

We use the LSC scheme to solve the LSESs and investigate its numerical behavior. According to the precise mathematical definition of the white noise [13, 14], we can simulate the noise as χ˙kn+(1/2)=(1/hτ)χkn+(1/2), where χkn+(1/2), k=0,1,…,N is a sequence of independent random variables with normal law 𝒩(0,1) at each time increment. Denote
(31)∥un∥∞=max⁡1≤k≤N|ukn|,eQn=Qn-Q0,eHn=Hn-Hn-1.
The numerical residuals of 𝒬(t) and ℋ(t) are measured by eQn and eHn, respectively. For numerical computation, we take τ=0.01, h=π/50, and ϵ=0.05,0.1.Example 1.
LSE with constant coefficients and periodic boundary condition.
Consider(32)iut+ux4-15u=ϵu∘χ˙,x∈[0,π],t>0,u(x,0)=exp⁡(iπ3)cos⁡⁡2x.
The exact solution of its deterministic system is u(x,t)=exp⁡[i(t+(π/3))]cos⁡⁡2x. The right side in the above system can be seen as a stochastic perturbation term.
Figure1 plots the amplitude |ukn| for one trajectory. Figure 2 shows the evolution of ∥u∥∞ for one trajectory and the average norm over 50 trajectories. We see that the white noise produces stochastic perturbation on the solitary wave and the size of perturbation depends on the size of noise. Figure 2 plots the residuals of discrete charge of one trajectory, which verifies the conservation of discrete charge of the compact schemes. Figure 3 plots the residuals of discrete energy for one trajectory and the average norm over 50 trajectories. The figure tells us that the stochastic noise makes residuals of discrete energy increase linearly over time.|
u
k
n
| for one trajectory.
(a)
ϵ
=
0.05
(b)
ϵ
=
0.1∥
u
∥
∞ for one trajectory and the average norm over 50 trajectories (a). Residuals of discrete charge for one trajectory (b).
(a)
(b)Residuals of discrete energy for one trajectory and the average energy over 50 trajectories.
(a)
e
H
n
(b)
H
n
-
H
0Example 2.
LSE with a variable coefficient and periodic boundary condition.
Consider(33)iut+ux4+(8*cot2x-7)u=ϵu∘χ˙,x∈(0,π),t>0,u(x,0)=sin2x.
The exact solution of its deterministic system is u(x,y,t)=eitsin2x.
Figure4 plots the amplitude |ukn| for one trajectory. Figure 5 shows the evolution of ∥u∥∞ for one trajectory and the average norm over 50 trajectories. We see that the white noise produces stochastic perturbation on the solitary wave and the size of perturbation depends on the size of noise. Figure 5 plots the residuals of discrete charge of one trajectory, which verifies the conservation of discrete charge of the compact schemes. Figure 6 plots the residuals of discrete energy for one trajectory and the average norm over 50 trajectories. The figure tells us that the stochastic noise makes residuals of discrete energy increase linearly over time.|
u
k
n
| for one trajectory.
(a)
ϵ
=
0.05
(b)
ϵ
=
0.1∥
u
∥
∞ for one trajectory and the average norm over 50 trajectories (a). Residuals of discrete charge for one trajectory (b).
(a)
(b)Residuals of discrete energy for one trajectory and the average energy over 50 trajectories.
(a)
e
H
n
(b)
H
n
-
H
0


## 5. Conclusion

In this paper, we apply a symplectic scheme in time and a kind of compact difference schemes in space to solve the LSES. The methods are unconditionally stable. Under periodic boundary conditions, they preserve a discrete charge invariant and satisfy a discrete energy transforming law.



---

*Source: 427023-2014-04-22.xml*