An Evolution of the Splitting from Cosine to Eikonal Schemes

An Evolution of the Splitting from Cosine to Eikonal Schemes

Qin Sheng

Department of Mathematics

Center for Astrophysics, Space Physics & Engineering Research Baylor University

Texas, USA

http://bearspace.baylor.edu/Qin Sheng/www/

... A few words about Baylor...

Founded in 1845, Baylor University is one of the biggest private universities in the United States. It consists of 11 academic colleges and has approximately 14,400 undergradu- ate students and 2,000 graduate students. The university locates in the heart of Texas, between Austin and Dallas. A research university, Baylor is one of the Big 12 Confer- ence universities in America.

Department of Mathematics is one of 27 departments in the College of Arts & Sci- ences. Research groups include pure, applied mathematics and mathematical educa-

1. The Idea of Splitting.

Although modern splitting methods can be really sophisticated, such as the cosine split- ting and eikonal splitting which are nonlinear transformation based splitting methods for wave equations, basic ideas of splitting is very simple. They can be traced back to the Baker-Campbell-Hausdorff formula (1897, 1906),

log (exp{ γ ^A }exp{ γ ^B }) = γ ^A + γ ^B + γ

²

2 [A,B] + γ

³

12 [A,[A,B]] − γ

³

12 [B,[A,B]]

− γ

⁴

48 [B,[A,[A,B]]] − γ

⁴

48 [A,[B,[A,B]]] + ⋅⋅⋅ ,

where

A , B

are from any Lie algebra

℘

^and

[A,B]

is the commutator of them.

The most popular classical applications of splitting may be the Alternating Direction Implicit (ADI) method:

u (x,y, ˜t) − u(x,y,t

k

) = γ ⁽ δ

x²

u (x,y, ˜t) + δ

y²

(x,y,t

k

) ) , u (x,y,t

k+1

) − u(x,y, ˜t) = γ ⁽ δ

x²

u (x,y, ˜t) + δ

y²

(x,y,t

k+1

) )

and Locally One Dimensional (LOD) method:

u (x,y, ˜t) − u(x,y,t

k

) = γ ⁽ δ

x²

u (x,y, ˜t) + δ

x²

(x,y,t

k

) )

,

^(1.1)

u (x,y,t

k+1

) − u(x,y, ˜t) = γ ⁽ δ

_y²

^u (x,y, ˜t) + δ

_y²

(x,y,t

k+1

) )

.

^(1.2)

They are for the solution of the linear convection-diffusion equation

u

_t

= a(u

xx

+ u

yy

), a > 0,

^(1.3)

together with proper IBVs.

Splitting methods, associated with finite differences, finite elements, hybrid multi- scale settings or adaptations, have been widely used for solving all the major classes of different differential equations. Latest challenging research issues in the area range from splitting for higher accuracy and flexibility in different parallel environments, split- ting for nonlinear partial differential equations, splitting for singular differential equations and inverse problems, nonlinear stability and convergence of splitting schemes, iterative and adaptive splitting strategies, geometric integration and domain decompositions, to quantum splitting computations in modern bio-chemistry and defense applications.

↬

Operator splitting

↬

Dimensional splitting

↬

Asymptotic formulas

↬

Adaptive splitting

↬

Cosine splitting

↬

Eikonal splitting...

A more commonly used splitting concept starts from a flow operator:

E (t) = exp{tX} = exp{t(X

1

+ X

2

+ X

3

+ ⋅⋅⋅ + X

n

)},

where each

X

_k

, 1 ≤ k ≤ n,

is simpler than the original vector field

X

in either of the following two ways:

1. The types of

X

_k are simpler. An example is that a conservation law may contain fast and slow wave terms which can be treated separately.

2. The

X

_k’s are easier to discuss numerically. For instance, dimensional splitting for multidimensional diffusion equations may be decomposed dimension-wisely. An- other example is the split-step Fourier method for the linear Schr¨odinger equation

iu = u

xx

+V (x)u,

where each term is linear and Hamiltonian, but the first term can be integrated more quickly if a splitting is employed.

ϕ (t) = ∑

^K

k=1

β

k

E

_k

(t),

^(1.4)

where

E

_k

(t)

are finite products of the exponentials

exp

{ α

k, j

tX

_{r(k, j)}

(t) }

, 1 ≤ r(k, j) ≤ n ,

^and

α

k, j

, β

k are nonnegative. Let

∥ ⋅ ∥

be a suitable norm. We want

∥ ϕ (t) − E(t)∥ = O (

t

^p+1

) .

Theorem 1.1 (Sheng–Suzuki, 1989). Let the real parts of the eigenvalues of

X

_k

, k = 1 ,2,...,n,

be nonpositive. If all

t , α

k, j

, β

k

≥ 0

^then

p ⩽ 2.

^{That is,} the highest order of a stable splitting is two.

Theorem 1.2 (Blanes and Casas, 2004). Let the real parts of the eigenvalues of

X

_k

, k = 1,2,...,n,

be nonpositive. If all

t , α

k, j

, β

k

∈ ℂ

⁺ then there is no bound

.

Frequently used splitting formulas include

ϕ

1

(t) = exp(tX

1

)exp(tX

2

)⋅⋅⋅exp(tX

n

), ϕ

2

(t) = 1

2 [exp(tX

1

)exp(tX

2

)⋅⋅⋅exp(tX

n

) + exp(tX

n

)exp(tX

n−1

)⋅⋅⋅exp(tX

1

)], ϕ

3

(t) = exp(tX

1

/2)exp(tX

2

/2)⋅⋅⋅exp(tX

n−1

/2)exp(tX

n

)

×exp(tX

n−1

/2)exp(tX

n−2

/2)⋅⋅⋅exp(tX

1

/2)).

It is easy to show that

∥ ϕ

1

(t) − E(t)∥ = O ( t

²

) ,

∥ ϕ

2

(t) − E(t)∥ = O ( t

³

) ,

∥ ϕ

3

(t) − E(t)∥ = O ( t

³

) .

Cosine splitting and eikonal splitting are extensions of the basic splitting concept for nonlinear wave equations. They are extremely useful methods built on top of nonlinear transformations. Such a splitting process often consists of FIVE stages:

1. transforming the original PDE problem into an integrable system with a phase space

M

and a vector field

X ,

u

^′

= X(u), u ∈ M.

^(1.5)

2. selecting a set of vector fields

X

_k

, k = 1,2,...,n,

^{such that}

X = ∑

ⁿk=1

X

_k

;

3. integrating either exactly or approximately each

X

_k

;

4. combining these solutions to yield an integrator for approximating

X ;

and 5. estimating, analyzing the numerical stability, conservation and error involved.

2. Cosine Splitting

For a typical nonlinear wave, we consider the sine-Gordon equation

u

_tt

= u

xx

+ u

yy

− φ (x,y)sinu, t > 0,

^(2.1)

in the region

Ω = {(x,y)∣−a < x < a, −b < y < b}.

The BCs associated with

(

^2.1

)

impose a zero gradient along the boundary

Γ

^of

Ω,

u

_x

= 0, x = ±a; u

y

= 0, y = ±b,

^(2.2)

while the initial conditions are

u (x,y,0) = f (x,y), u

t

(x,y,0) = g(x,y).

^(2.3)

The function

φ

can be interpreted as a Josephson current density, and

f

and

g

are wave modes or kinks and their velocity, respectively.

Let

u

_{i+1, j,k}

− 2u

i, j,k

+ u

i−1, j,k

h

²_x

, u

_{i, j+1,k}

− 2u

i, j,k

+ u

i, j−1,k

h

²_y

⇒ u

xx

(x

i

,y

j

,t

k

), u

yy

(x

i

,y

j

,t

k

), u

_{i+1, j,k}

− u

i−1, j,k

2h

_x

, u

_{i, j+1,k}

− u

i, j−1,k

2h

_y

⇒ u

x

(x

i

,y

j

,t

k

), u

y

(x

i

,y

j

,t

k

).

Denote

B = h

⁻²_x

B

₁

+ h

⁻²_y

B

₂

∈ ℝ

^mn×mn

,

where

B

₁

= I ⊗ T, B

2

= T ⊗ I

and

T =

⎡

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎢ ⎢

⎣

−2 2

1 −2 1

⋅⋅⋅

⋅⋅⋅

1 −2 1 2 −2

⎤

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎥ ⎥

⎦

is a

n ×n

tridiagonal matrix and

I ∈ ℝ

^m×m is the identity matrix. We obtain the second- order semidiscretized system from

(

^2.1

)

^,

(

^2.2

)

^:

u

^′′_k

= Bu

k

+ r(u

k

).

^(2.4)

Now, from

(

^2.4

)

we acquire that

u

_k+1

− 2u

k

+ u

k−1

= τ

²

ψ ( τ

²

^A )Bu

k

+ τ

²

ψ ( τ

²

^A )r(u

k

)

^(2.5)

together with

(

^2.3

)

^{, where}

τ > 0

^and

A = ∂

∂ ^u (Bu + r(u))  

u=uk

= B + r

u

(u

k

), ψ (z) = cos √

−z − 1 z /2 , r

_u

(u

k

) =

^diag

[

φ

1,1

cos u

_1,1,k

, φ

1,2

cos u

_1,2,k

, ...,

φ

1,n

cos u

_1,n,k

, φ

2,1

cos u

_2,1,k

, ..., φ

m,n

cos u

_m,n,k

] .

Recall that

cos √

−z =

^cosh

√ z .

From

(

^2.5

)

, we obtain the cosine scheme

u

_k+1

− 2

^cosh

( τ ^{√ A} )u

k

+ u

k−1

= τ

²

ψ ( τ

²

^A )[r(u

k

) − r

u

(u

k

)u

k

].

^(2.6)

But how to link

(

^2.6

)

to splitting?

Note that

cosh

( τ ^{√ A} ) = e

^{τ22 A}

+ O( τ

⁴

)

when the real parts of the eigenvalues of

A

are negative. Thus,

(

^2.6

)

can be approxi- mated by

u

_k+1

− 2e

^{τ22 A}

u

_k

+ u

k−1

= τ

²

(r(u

k

) − r

u

(u

k

)u

k

)

^(2.7)

incurring a local error

O ( τ

⁴

).

Therefore we don’t need to evaluate the matrix exponential of

√

A

directly!

Set

A

₁

= 1

h

²_x

B

₁

− 1

2 r

_u

(u

k

), A

2

= 1

h

²_y

B

₂

− 1

2 r

_u

(u

k

).

Thus,

A = A

1

+ A

2

.

Recall Strang splitting

ϕ

³

( τ ),

^{we have}

e

^{τ22 A}

= e

^{τ24 A1}

e

^{τ22 A2}

e

^{τ24 A1}

+ O( τ

⁶

).

A substitution of it to

(

^2.7

)

^yields

u

_k+1

− 2e

^{τ24 A1}

e

^{τ22 A2}

e

^{τ24 A1}

u

_k

+ u

k−1

= τ

²

(r(u

k

) − r

u

(u

k

)u

k

).

Now, replace the matrix exponential functions in the above equation by

[0/1], [1/1]

and

[1/0]

Pad ´e appromants, respectively. We obtain the following two-stage cosine splitting scheme:

(

I − τ

²

4 A

₁

)

v

_k

= (

I + τ

²

4 A

₂

)(

I − τ

²

4 A

₂

)

₋₁

(

I + τ

²

4 A

₁

)

u

_k

,

^(2.8)

u

_k+1

= 2v

k

+ τ

²

(r(u

k

) − r

^u

(u

k

)u

k

) − u

k−1

, k = 1,2,...

^(2.9)

Lemma 3.1. All eigenvalues of

B

₁

, B

2

,

^and

A

₁

, A

2 are real.

Lemma 3.2. For

i = 1,2,...,mn,

^{we have}

− 1

2 (r

u

(u

k

))

_i

− 4

h

²_x

≤ λ

_i^(A1)

≤ − 1

2 (r

u

(u

k

))

_i

,

− 1

2 (r

u

(u

k

))

_i

− 4

h

²_y

≤ λ

_i^(A2)

≤ − 1

2 (r

u

(u

k

))

_i

.

Further, based on the Bauer-Fike bound,

Lemma 3.3. Let

λ

^{(B1) and}

λ

^(B2) be eigenvalues of ¹

h²_x

B

₁ and ¹

h²_y

B

₂

,

respectively. We have

min

i

{  λ

_i^(A1)

− λ

^(B1)

  }

≤ κ (Q

1

V

₁

)max

j

 (r

^u

(u

k

))

_j

 , min {  λ

^(A2)

− λ

^(B2)

  }

≤ κ (Q

2

V

₂

)max  (r

^u

(u

k

))  ,

where

κ (⋅)

is the condition number,

V =

^diag

(1,2,...,2,1) ∈ ℝ

^n×n

, V

₁²

= I

m

⊗V, V

₂²

= V ⊗ I

m

, I

_m

∈ ℝ

^{m×m and}

Q

₁

, Q

2 are two orthogonal matrices.

Rewrite

(

^2.8

)

^,

(

^2.9

)

in the embedded form:

w

_k

= Pu

k

,

^(3.1)

u

_k+1

= 2w

k

− u

k−1

+ τ

²

(r(u

k

) − r

u

(u

k

)u

k

), k = 1,2,...,

^(3.2)

where

P = (

I − τ

²

4 A

₁

)

₋₁

(

I + τ

²

4 A

₂

)(

I − τ

²

4 A

₂

)

₋₁

(

I + τ

²

4 A

₁

)

is the Peaceman-Rachford splitting operator.

Consequently, the perturbed equation for the numerical error is

v

_k+1

= Qv

k

,

where

v

_k

=

⎛

⎝ u

_k

− ˜u

k

u

_k−1

− ˜u

k−1

⎞

⎠, Q =

⎛

⎝ 2P −I I 0

⎞

⎠

and

u ˜

_k is the perturbed solution.

Denote

κ

¹

= τ

²

4h

²_x

[

4 + h

²_x

max

i

{(r

u

(u

k

))

_i

} ]

, κ

2

= τ

²

4h

²_y

[

4 + h

²_y

max

i

{(r

u

(u

k

))

_i

} ]

.

Theorem 3.4. Let

(r

u

(u

k

))

_i

≥ 0, i = 1,2,...,mn.

^If

κ

1

, κ

2

≤ 2,

then the two-stage cosine splitting method

(

^3.1

)

^,

(

^3.2

)

is stable in the weak von Neumann sense.

Further, denote

σ

1

= max

i

  

2 + τ

²

λ

_i^(A1)

2 − τ

²

λ

_i^(A1)

 



^and

σ

2

= max

i

  

2 + τ

²

λ

_i^(A2)

2 − τ

²

λ

_i^(A2)

   .

Theorem 3.5. If

2 σ

1

σ

2

≤ 1,

then the two-stage cosine splitting scheme

(

^3.1

)

^,

(

^3.2

)

^is

stable in the von Neumann sense.

Theorem 3.6. Let

λ

_i^{(A1) and}

λ

_i^(A2)

, i = 1,2,...,mn,

be eigenvalues of

A

₁ and

A

₂

,

respectively. If

−2 (√

2 + 1 )

₂

≤ τ

²

λ

_i^(Aℓ)

≤ −2 (√

2 − 1 )

₂

, ℓ = 1,2, i = 1,2,...,nm,

^(3.3)

then

(

^3.1

)

^,

(

^3.2

)

is stable in the von Neumann sense.

Theorem 3.7. Let

λ

_i^{(A1) and}

λ

_i^(A2)

, i = 1,2,...,mn,

be eigenvalues of

A

₁ and

A

₂

,

respectively. A necessary condition for

(

^3.3

)

^{to hold is}

τ

²

h

²_x

, τ

²

h

²_y

≤ 2 √

2 .

The total energy involved in the solitary waves can be computed via the integral

∫ ∫

[u

t

u

_tt

− u

t

(u

xx

+ u

yy

) + u

t

φ (x,y)sinu]dxdy = 0.

Integrating by parts, we obtain

1 2

∂

∂ ^t

∫ ∫

[u

xx

+ u

yy

+ u

tt

+ 2(1 − cosu)]dxdy = 0.

This leads to the undamped energy formula

E = 1 2

∫ ∫

[u

xx

+ u

yy

+ u

tt

+ 2(1 − cosu)]dxdy.

^(3.4)

Via an appropriate numerical quadrature, we may show that the two-stage cosine split- ting method

(

^2.8

)

^,

(

^2.9

)

is numerically conservative.

The strategy can be extended for solving other nonlinear wave equations including Schr ¨odinger and the KdV equations.

4. Eikonal Splitting

Think in optics. According to Maxwell’s field equations, a slowly varying envelope ap- proximation of the light beam can be ideally modeled by a three-dimensional paraxial Helmholtz equation,

2i κ ^u

^z

= u

xx

+ u

yy

, 0 ≤ x,y ≤ ℓ, z ≥ z

0

,

^(4.1)

where

i = √

−1, κ = 2 π / λ

is the wave number,

λ

is the wavelength,

z

is the beam propagation direction and

x ,y

represent the spatial Cartesian coordinates perpendicu- lar to the light. Further,

u

is the complex envelope of the wave function and

ℓ

is a finite boundary location needed by computations. Since

κ ≥ 10

⁵ in most important applica- tions,

u

is highly oscillatory. As a consequence,

(

^4.1

)

is difficult to solve effectively via any conventional methods.

0.8 0.9 1 1.1 1.2 1.3 1.4 1.5

−120

−100

−80

−60

−40

−20 0 20 40 60

z

v = real(u)

real part of the solution

FIGURE 4.1. A typical highly oscillatory nonlinear wave in optical beam propagation studies.

Take the eikonal transformation in geometric optics,

u (x,y,z) = φ (x,y,z)exp{i κψ (x,y,z)},

^(4.2)

where

φ , ψ

are nonoscillatory real functions and

φ ∕= 0.

^Surfaces

φ =

c are often called wavefronts of the disturbance. Substitute

(

^4.2

)

^into

(

^4.1

)

and separate real and imaginary parts of resulted equation. We acquire readily the following coupled ray

equations

φ

z

= α ( ψ

xx

+ ψ

yy

) + f

1

,

^(4.3)

ψ

^z

= β ( φ

^xx

+ φ

^yy

) + f

2

,

^(4.4)

where

α = φ

2 , β = − 1 2 κ

²

φ ^,

f

₁

= φ

x

ψ

x

+ φ

y

ψ

y

, f

2

= 1 2

( ( ψ

x

)

²

+ ( ψ

y

)

²

) .

Note that solutions of the coupled nonlinear equations

(

^4.3

)

^,

(

^4.4

)

are oscillation-free even for large values of

κ .

Equations

(

^4.3

)

^,

(

^4.4

)

are not easy to solve, however. Many explicit or semi-explicit schemes, such as those suggested by Richtmyer and Morton (1994, 2003), perform poorly. This cries out for splitting!

w = ⎝ φ

ψ ^{⎠, f = ⎝} f

₁

f

₂

⎠, M = ⎝ 0 α

β ^{0 ⎠.}

Lemma 4.1. The matrix

M

is similar to a skew symmetric matrix and thus its eigen- values are pure imaginary. Further, for any

φ ∕= 0,

the matrix has a pair of constant eigenvalues

λ

^M

= ± i/(2 κ )

and thus a spectral radius

ρ (M) = 1/(2 κ ).

Rewrite

(

^4.3

)

^,

(

^4.4

)

in a matrix form

w

_z

= Mw

xx

+ Mw

yy

+ f.

^(4.5)

Different BCs may be introduced based on the assumption that light beams are con- centrated symmetrically around its geometric center and the light intensity is negligible near the boundaries in both

x

and

y

directions. Examples:

w (x,y,z

0

) = g

0

(x,y),

^(4.6)

w

_x

(0,y,z) = w

x

(ℓ,y,z) = 0, w

y

(x,0,z) = w

y

(x,ℓ,z) = 0.

^(4.8)

Recall the LOD scheme. We obtain from

(

^4.5

)

^-

(

^4.7

) w

^r+1/2

− w

^r

= h

_z

2 (

Mw

^r_xx^+1/2

+ Mw

^r_xx

)

+ h

_z

2 f

^r

, w

^r+1

− w

^r+1/2

= h

_z

2 (

Mw

^r_yy^+1/2

+ Mw

^r_yy⁺¹

)

+ h

_z

2 f

^r+1/2

.

The above can be formulated to

(I − μ ^M

^r

)w

^r+1/2

= (I + μ ^M

^r

)w

^r

+ h

_z

2 f

^r

,

^(4.9)

(

I − η ^M

^r+1/2

^{) v}

^r+1

= (

I + η ^M

^r+1/2

^{) v}

^r+1/2

+ h

_z

2 g

^r+1/2

,

^(4.10)

M

^σ

=

^diag

(M

₁^σ

,M

₂^σ

,...,M

_n^σ

),

M

^σ_j

=

⎛

⎜ ⎜

⎜ ⎜

⎜ ⎜

⎜ ⎜

⎜ ⎜

⎜ ⎜

⎝

−2M

₁^σ_{, j}

M

₁^σ_{, j}

M

₂^σ_{, j}

−2M

₂^σ_{, j}

M

₂^σ_{, j}

M

₃^σ_{, j}

−2M

₃^σ_{, j}

M

₃^σ_{, j}

⋅⋅⋅ ⋅⋅⋅ ⋅⋅⋅

M

_n^σ_{−1, j}

−2M

_n^σ_{−1, j}

M

_n^σ_{−1, j}

M

_n^σ_{, j}

−2M

_n^σ_{, j}

⎞

⎟ ⎟

⎟ ⎟

⎟ ⎟

⎟ ⎟

⎟ ⎟

⎟ ⎟

⎠ ,

j = 1,2,...,n.

Note that

v

^σ

= Pw

^σ

, g

^σ

= P f

^σ

,

^(4.11)

where

P

is a permutation matrix. Therefore

(

^4.10

)

can be replaced by

( − η

^r+1/2

⁾

^r+1

= (

+ η

^r+1/2

⁾

^r+1/2

+ h

_{z r+1/2}

,

where

N

^σ

= P

⁻¹

M

^σ

P .

Lemma 4.2. The eigenvalues of matrices

M

^σ

, N

^σ are pure imaginary.

Let

S

^σ

=

^diag

(N

₁^σ

,N

₂^σ

,...,N

_n^σ

) ∈ ℝ

^2n2×2n2

,

where

N

^σ_j

=

^diag

(

M

₁^σ_{, j}

,M

₂^σ_{, j}

,...,M

_n^σ_{, j}

)

, j = 1,2,...,n,

and

T

₀

=

⎛

⎜ ⎜

⎜ ⎜

⎜ ⎜

⎜ ⎜

⎜ ⎜

⎜ ⎜

⎝

−2I

2

I

₂

I

₂

−2I

2

I

₂

I

₂

−2I

2

I

₂

⋅⋅⋅ ⋅⋅⋅ ⋅⋅⋅

I

₂

−2I

2

I

₂

I

₂

−2I

2

⎞

⎟ ⎟

⎟ ⎟

⎟ ⎟

⎟ ⎟

⎟ ⎟

⎟ ⎟

⎠

∈ ℝ

^2n×2n

,

T

₁

=

^diag

(T

0

,T

0

,...,T

0

) ∈ ℝ

^2n2×2n2

.

Lemma 4.3. We have

ρ (S

^σ

) = O(1/ κ ),

^where

κ

is the wave number.

Definition 4.4. Let

A

be the amplifying matrix of a finite difference scheme, such as

(

^4.9

)

^,

(

^4.10

)

, for solving an highly oscillatory problem. We say that the scheme is asymptotically stable if there exists

c > 0

^{such that}

ρ (A) = 1 + O

( 1 κ

^c

)

,

^(4.12)

where

κ

is the wave number. Further,

c

is called the asymptotical stability index of the scheme. If

(

^4.12

)

holds for all mesh steps

h

_x

,h

y and

h

_z then we say that the scheme is unconditionally asymptotically stable.

Theorem 4.5. The eikonal splitting scheme

(

^4.9

)

^,

(

^4.10

)

is unconditionally asymptoti- cally stable with a index 2.

↬

Similar results can be obtained with the Neumann BCs

(

^4.8

)

^.

↬

Discussions can be carried out for the eikonal ADI method:

w

^r_s^+1/2_,t

− w

^r_s,t

= h

_z

2 M

_s^r_,t

( δ

x²

w

^r_s^+1/2_,t

+ δ

y²

w

^r_s,t

) + h

_z

2 f

_s^r_,t

, w

^r_s⁺¹_,t

− w

^r_s,t^+1/2

= h

_z

2 M

_s^r_,t^+1/2

( δ

x²

w

^r_s^+1/2_,t

+ δ

y²

w

^r_s⁺¹_,t

) + h

_z

2 f

_s^r_,t^+1/2

.

However, analysis of the numerical stability of the above splitting method can be sig- nificantly different due to the fact that structures of amplification matrices utilized are different. In fact, in our latest efforts, investigations are carried out via nested Kronecker products of the matrices involved.

We may comprise our linear system

(

^4.9

)

^,

(

^4.10

)

in a different way, that is,

I φ

^r+1/2

− μ ^D

_α^r

^S ψ

^r+1/2

= g

^r₁

,

^(5.1)

I ψ

^r+1/2

− μ ^D

_β^r

^S φ

^r+1/2

= g

^r₂

,

^(5.2)

where

I ∈ ℝ

^n2×n2 is the identity matrix.

D

_α^σ

=

^diag

( α

₁^σ_,1

, α

₂^σ_,1

,..., α

_n^σ_,1

, α

₁^σ_,2

, α

₂^σ_,2

,..., α

_n^σ_,2

,..., α

₁^σ_,n

, α

₂^σ_,n

,..., α

_n^σ_,n

), D

_βσ

=

^diag

( β

₁^σ_,1

, β

₂^σ_,1

,..., β

_n^σ_,1

, β

₁^σ_,2

, β

₂^σ_,2

,..., β

_n^σ_,2

,..., β

₁^σ_,n

, β

₂^σ_,n

,..., β

_n^σ_,n

), S =

^diag

(S

_θ

,S

_θ

,...,S

_θ

) ∈ ℝ

^{n2×n2 in which}

S

_θ

=

^tridiag

(1, θ

k

,1) ∈ ℝ

^{n×n with}

θ

1

= θ

n

= −1; θ

k

= −2, k = 2,3,...,n − 1,

^and

φ

^σ

= (

φ

₁^σ_,1

, φ

₂^σ_,1

,..., φ

_n^σ_,1

, φ

₁^σ_,2

, φ

₂^σ_,2

,..., φ

_n^σ_,2

,..., φ

₁^σ_,n

, φ

₂^σ_,n

,..., φ

_n^σ_,n

⁾

^T

,

^(5.3)

ψ

^σ

= (

ψ

₁^σ_,1

, ψ

₂^σ_,1

,..., ψ

_n^σ_,1

, ψ

₁^σ_,2

, ψ

₂^σ_,2

,..., ψ

_n^σ_,2

,..., ψ

₁^σ_,n

, ψ

₂^σ_,n

,..., ψ

_n^σ_,n

⁾

^(5.4)T

On the other hand,

g

^σ₁

, g

^σ₂ are corresponding right-hand-sides in

(

^??

)

^and

(

^??

)

^,

that is,

g

^σ₁

= I φ

^σ

+ μ ^D

_α^σ

^S ψ

^σ

+ h

_z

2 f

₁^σ

, g

^σ₂

= I ψ

^σ

+ μ ^D

_β^σ

^S φ

^σ

+ h

_z

2 f

₂^σ

,

with

f

₁^σ

= (( f

1

)

^σ_1,1

,( f

1

)

^σ_2,1

,...,( f

1

)

^σ_n,1

,( f

1

)

^σ_1,2

,( f

1

)

^σ_2,2

,...,( f

1

)

^σ_n,2

,...,( f

1

)

^σ_1,n

,( f

1

)

^σ_2,n

, ...,( f

1

)

^σ_n,n

)

^{T and}

f

₂^σ

= (( f

2

)

^σ_1,1

,( f

2

)

^σ_2,1

,...,( f

2

)

^σ_n,1

,( f

2

)

^σ_1,2

,( f

2

)

^σ_2,2

,...,( f

2

)

^σ_n,2

,..., ( f

2

)

^σ_1,n

, ( f

2

)

^σ_2,n

,...,( f

2

)

^σ_n,n

)

^T

.

Note that

φ

^r+1/2

= μ ^D

α^r

S ψ

^r+1/2

+ g

^r₁

.

A

_r

ψ

^r+1/2

= μ ^D

_β^r

^Sg

^r₁

+ g

^r₂

,

^(5.5)

φ

^r+1/2

= μ ^D

α^r

S ψ

^r+1/2

+ g

^r₁

,

^(5.6)

B

_r+1/2

ψ ˜

^r+1

= η ^D

_β^r+1/2

^Sg

^r₃^+1/2

+ g

^r₄^+1/2

,

^(5.7)

φ ˜

^r+1

= η ^D

_α^r+1/2

^{S ˜} ψ

^r+1

+ g

^r₃^+1/2

,

^(5.8)

r = 0,1,2,...,

where the coefficient matrices

A

_r

= I − μ

²

^D

_β^r

^SD

_α^r

^S , B

r+1/2

= I − η

²

^D

_β^r+1/2

^SD

_α^r+1/2

^S

are quintic diagonal and

φ ˜

^σ

= (

φ

₁^σ_,1

, φ

₁^σ_,2

,..., φ

₁^σ_,n

, φ

₂^σ_,1

, φ

₂^σ_,2

,..., φ

₂^σ_,n

,..., φ

_n^σ_,1

, φ

_n^σ_,2

,..., φ

n^σ,n

)

_T

, ψ ˜

^σ

= (

ψ

₁^σ_,1

, ψ

₁^σ_,2

,..., ψ

₁^σ_,n

, ψ

₂^σ_,1

, ψ

₂^σ_,2

,..., ψ

₂^σ_,n

,..., ψ

_n^σ_,1

, ψ

_n^σ_,2

,..., ψ

n^σ,n

)

_T

.

Theorem 5.1. There exist reasonable values of

μ , η

^{such that}

^A

^{r and}

^B

r+1/2 are nonsingular. Therefore the solution of

(

^5.5

)

^-

(

^5.8

)

exists and is unique.

↬

Parallelization

↬

Iterative procedures

↬

Acceleration

↬

^HPC...

6. Simulation Experiments

6.1. Solitary waves via sine-Gordon equations.

−10

−5 0

5 10

−10

−5 0 5 10 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

x t = 0

y

sin(u/2)

−10

−5 0

5 10

−10

−5 0 5 10 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

x t = 4

y

sin(u/2)

−10

−5 0

5 10

−10

−5 0 5 10

−0.5

−0.4

−0.3

−0.2

−0.1 0 0.1 0.2 0.3 0.4

x t = 8

y

sin(u/2)

−10

−5 0

5 10

−10

−5 0 5 10 0 0.2 0.4 0.6 0.8

x t = 11.5

y

sin(u/2)

−10

−5 0

5 10

−10

−5 0 5 10

−0.8

−0.6

−0.4

−0.2 0 0.2

x t = 15

y

sin(u/2)

FIGURE6.1. Circular and elliptic ring solitons: the functionsin(u/2) (from left to right:

t = 0, 4, 8, 11.5 ^and 15).