Introduction

In this thesis, we consider the solution of a nonlinear Schr¨odinger equation [1]:

−∆_qu(x) + α₁|u(x)|²u(x) = λ₁u(x), (1.1) for x ∈ Ω and u : Ω → C with normalized constraint

Z

Ω

|u(x)|²dx = 1, (1.2)

where Ω ⊆ R² is a bounded domain, and the Dirichlet boundary condition

u(x) = 0, x ∈ ∂Ω. (1.3)

According to [1], we define





∇ = (∂x, ∂y)^>,

∇_q = ∇ + iγq,

−∆_q = < ∇_q, ∇_q >,

(1.4)

where q = (q₁, q₂)^> represents a vector on Ω , and γ is the magnitude of magnetic field. From equation (1.4), we can get

−∆_q =< ∇_q, ∇_q >

=< ∇ + iγq, ∇ + iγq >

=< ∇, ∇ > + < iγq, ∇ > + < ∇, iγq > + < iγq, iγq >

= −∆ + iγ∇^∗q − iγq^∗∇ + γ²q^∗q, (1.5) where ∆_q is magnetic laplacian operator[2, 3].

In section 2, we introduce the nonlinear Schr¨odinger equation. We put uc = u1+ iu2 into equation (1.1). We can get the new equation. To consider the energy functional and the eigenvalue, we can find the relation between them. After verifying the relation, we want to transform continuous differ-ential equation to discrete equation in section 3. As a result, we use finite difference method to approximate the solution. In other words, we have to discretize the new equation. We can see details in section 3. Now, we have the discrete nonlinear Schr¨odinger equation. We denote the discrete equation as G(u(s), γ(s)) = 0, where s is an arc length parameter.

We use the fixed point iteration method in subsection 4.1 so as to obtain the initial solution of the discrete equation. Then, we use the continuation method to approximate the solution curve. We will explain the conception of standard continuation method in subsection 4.2. Continuation method

includes predictions and corrections. In the part of predictions, we have to find the tangent vector of possible solution on solution curve. Here, we differentiate the discrete equation with respect to s. Then, we can get the system as follows.

£ G_u G_γ ¤·

˙u(s)

˙γ(s)

¸

= 0.

By solving the system above, we can obtain the tangent vector. We suppose that (u₀, γ₀) is the initial solution. And we denote (u_i, γ_i) as the possible solution. After computing, we can get the new point (u⁽¹⁾_i , γ_i⁽¹⁾) in predictions.

We continue to discuss the part of corrections. The main idea is Newton’s method. We consider the plane N_i whose normal vector is equal to the tangent vector. We hope that the point (u⁽¹⁾_i , γ_i⁽¹⁾) can converge to (u_i, γ_i) on the plane Ni iteratively. There are some problems on the computing process. In order to solve those problems, we have to introduce the idea about the rotation invariance. The solution set contains transformation invariant solution. Hence, we concentrate on solving the quotient solution. As a result, we will develop the quotient invariant transformation continuation method in subsection 4.3. We will add a plane to the system which we have to solve.

And we can obtain the meaningful vectors of predictions. We also introduce a small skill on computing so as to overcome some disadvantages. In subsection 4.4, we show the algorithm of continuation method. Eventually, we can get the approximate solution curve by repeating the steps of continuation method. We can see that the numerical results in section 5. And we have a conclusion in section 6.

2 Nonlinear Schr¨ odinger Equation

We use

q =

· −r cos θ r sin θ

¸ ρ,

where

ρ =









0 if r < ², 1 + cos(^5π_ar + π)

2 if ² ≤ r ≤ a,

1 if r ≥ a.

We substitute equation (1.5) into equation (1.1). Then, we can get [−∆ + iγ(∇^∗q − q^∗∇) + γ²kqk²+ α₁|u|²]u = λ₁u. (2.1)

We let u_c= u₁+ iu₂, where u₁, u₂ : Ω → R. Then we can obtain [−∆ + γ²kqk²+ α₁(u²₁+ u²₂)]u₁− γ(∇^∗q − q^∗∇)u₂ = λ₁u₁,

[−∆ + γ²kqk²+ α₁(u²₁+ u²₂)]u₂+ γ(∇^∗q − q^∗∇)u₁ = λ₁u₂. (2.2) We rewrite equation (2.1) into below form.

−∆u + iγ∂_θ1u + V u + α₁|u|²u = λ₁u. (2.3) The energy functional [7, 8, 9] is

E(φ) = Z

Ω

(|∇φ|²+ iγφ^∗∂θ1φ + V φ²+1

2α1|φ|⁴). (2.4) Let equation (2.3) be multiplied by u^∗. Then we can get equation as following.

−u^∗∆u + iγu^∗∂θ1u + u^∗V u + α1u^∗|u|²u = λ1u^∗u. (2.5) To integrate equation (2.5), it will become

Z

(−u^∗∆u + iγu^∗∂_θ1u + u^∗V u + α₁u^∗|u|²u) = Z

λ₁u^∗u. (2.6) Here,

− Z

u^∗∆u = − µ

u^∗∇u|_∂Ω− Z

∇u · ∇u

¶

= Z

|∇u|². (2.7) We put equation (2.7) into equation (2.6). From equation (1.2), we can get R

u^∗u = 1. Equation (2.6) can be represented as Z

(|∇u|²+ iγu^∗∂_θ1u + u^∗V u + α₁u^∗|u|²u) = λ₁ (2.8) Comparing the energy funtional E(φ) to eigenvalue λ₁ from equations (2.4) and (2.8), we can find that E(φ) = λ₁ −α1

2 Z

|u|⁴.

3 Discretization

In this section, we concern about the process of discretization. In order to transform continuous differential equations into discrete equations, we intro-duce finite difference method [4] to approximate the solutions of equations (2.2), (1.2) and (1.3). First, we divide the solutions into grids of nodes. The domain Ω which we consider is a disc. We cut the disc into n radial parts and m azimuthal parts. We let ∆r = 2

2n + 1, ∆θ = 2π

m. And we define the grid points as below.

r_i = (i − 1

2)∆r, i = 1, 2, · · · , n, θ_j = (j − 1)∆θ, j = 1, 2, · · · , m.

Here, we see an example for n = 3 and m = 10. We show the discrete points on the domain which is a disc in Figure 1.

u1,1

We use the coordinate transformation to convert Cartesian coordinate system (x, y) into the polar coordinate system (r, θ). Therefore, we can obtain

−∆u = −(u_xx + u_yy)

= −(urr+ 1

rur+ 1

r²uθθ). (3.1)

Now, we first consider the approximate solution of −∆u. We consider lapla-cian operator on disc with Dirichlet conditions as following:

½ −∆u = f (r, θ),

u(r, θ) = C(r, θ), (3.2)

where (r, θ) ∈ Ω.

Next, we replace laplacian by central difference formula as following u_rr = u(r_i+1, θ_j) − 2u(r_i, θ_j) + u(r_i−1, θ_j)

∆r² u_r = u(r_i+1, θ_j) − u(r_i−1, θ_j)

2∆r

u_θθ = u(ri, θj+1) − 2u(ri, θj) + u(ri, θj−1)

∆θ² . (3.3)

Then, we substitute equation (3.3) into equations (3.1) and (3.2). We get

− [ 1

∆r²(u(r_i+1, θ_j) − 2u(r_i, θ_j) + u(r_i−1, θ_j)) + 1

2r_i∆r(u(r_i+1, θ_j) − u(r_i−1, θ_j))

+ 1

r_i²∆θ²(u(r_i, θ_j+1) − 2u(r_i, θ_j) + u(r_i, θ_j−1))] = f (r_i, θ_j). (3.4) In order to simplify equation (3.4), we rewrite u(ri, θj) to ui,j and f (ri, θj) to f_i,j. And we substitute r_i = (i − ¹₂)∆r into equation (3.4). Then, we can obtain the following

− 1

∆r²[(ui+1,j− 2ui,j+ ui−1,j) + 1

2(i −¹₂)(ui+1,j− ui−1,j) + 1

(i −¹₂²)∆θ²(ui,j+1− 2ui,j+ ui,j−1)] = fi,j.

Let c_i = 1

2(i −¹₂), d_i = 1

(i − ¹₂)²∆θ². Now, we can get the following:

(i, j) fi,j

(1, 1) − 1

∆r²[(−2 − 2d1)u1,1+ d1u1,2+ d1u1,0+ (1 + c1)u2,1+ (1 − c1)u0,1] = f1,1

(1, 2) − 1

∆r²[(−2 − 2d1)u1,2+ d1u1,3+ d1u1,1+ (1 + c1)u2,2+ (1 − c1)u0,2] = f1,2

(1, 3) − 1

∆r²[(−2 − 2d1)u1,3+ d1u1,4+ d1u1,2+ (1 + c1)u2,3+ (1 − c1)u0,3] = f1,3

(1, 4) − 1

∆r²[(−2 − 2d1)u1,4+ d1u1,5+ d1u1,3+ (1 + c1)u2,4+ (1 − c1)u0,4] = f1,4

... ... (1, m − 1) − 1

∆r²[(−2 − 2d1)u1,m−1+ d1u1,m+ d1u1,m−2+ (1 + c1)u2,m−1+ (1 − c1)u0,m−1] = f1,m−1

(1, m) − 1

∆r²[(−2 − 2d1)u1,m+ d1u1,1+ d1u1,m−1+ (1 + c1)u2,m+ (1 − c1)u0,m] = f1,m

(2, 1) − 1

Then, we want to solve above difference equations related to Dirichlet boundary conditions u(r, θ) = c(r, θ). That is, we solve the linear system of Au = b. We get u_n+1,1 = c_n+1,1 from equation (3.2) and 1 − c₁ = 1 −

b =

We use the central difference formula as following u_θ = u_i,j+1− u_i,j−1

2∆θ .

After computing, we can get

D_θ = 1 matrix of γ²kqk² can be represented as following

γ²· D_q = γ²

where

By using formula for area of sector ¹₂∆r²∆θ, we can get as following

i = 1 1

Hence, we obtain the discrete matrix is C^>· U⁰ = 1.

C^>= ¹₂∆r²∆θ£ ₁

4 · · · ¹₄ 2 · · · 2 · · · 2(i − 1) · · · 2(i − 1)) ¤ ,

U⁰ =



 U¯1

¯...

U_n



, where ¯U_i =





u²₁(ri, θ1) + u²₂(ri, θ1) ...

u²₁(r_i, θ_m) + u²₂(r_i, θ_m)



 with i = 1, 2, · · · , n.

Finally, the discrete nonlinear Schr¨odinger equation relative to equation (2.3) can be represented as

A₁u_c+ iγSu_c+ α₁u¯_c◦ u^°_c² = λ₁u_c,

C^>(¯u_c◦ u_c) = 1. (3.8)

In (3.8), A₁ is the discrete matrix for the operator of −∆ + V . S is the discrete matrix for the operator of ∂_θ1 and S is a skew-symmetric matrix. We denote ◦ as Hadamard product. We can find that C^>(¯uc◦ uc) will be pure imaginary. Therefore, we can get u^°_c² = u_c◦u_c, u^°₁² = u₁◦u₁ and u^°₂² = u₂◦u₂. We also can obtain the energy functional corresponding to (2.4) is

E(φ) = C^>(¯u_c◦ A₁u_c+ iγ ¯u_c◦ Su_c+α₁

2 u¯^°_c² ◦ u^°_c² ).

4 Continuation Methed

4.1 Fixed point iteration method

Now, we want to get the initial solution of equation (3.8). Hence, we use fixed point iteration method [6] so as to converge the solution. We let u_c= u₁+iu₂, where u1,u2 ∈ R^N and N = n · m. And we put it into equation (3.8). Then, we can rewrite equation (3.8) into the form as follows.

(A₁ + α₁U)u₁− γSu₂ = λ₁u₁, (A1 + α1U)u2+ γSu1 = λ1u2,

C^>U⁰ = 1. (4.1)

We let γ = 0. The equations (4.1) will become the form as below.

(A + α₁U)u₁ = λ₁u₁, (A + α₁U)u₂ = λ₁u₂,

C^>U⁰ = 1. (4.2)

Then, we let u₂ = 0. Equations (4.2) can be simplified as

(A + α₁U)u₁ = λ₁u₁, (4.3a)

C^>u^°₁² = 1. (4.3b) Now, we consider the eigenvalue problem of equation (4.3a). First, we let u⁽⁰⁾ ∈ R^N be a random vector. We choose c₁ = C^> · u^(0)°2 such that u⁽⁰⁾ = ^√1_c

1u⁽⁰⁾. u⁽⁰⁾ will satisfy equation (4.3b), that is, C^>u^(0)°2 = C^>1

c₁u^(0)°2 = C^> 1

C^>u^(0)°2 u^(0)°2 = 1.

Next, we substitute equation above into equation (4.3a). To solve the eigenvalue problem, we can obtain eigenvalues and eigenvectors. We choose the smallest eigenvalue λ⁽¹⁾₁ and corresponding eigenvector u⁽¹⁾₁ and let u⁽¹⁾₁ satisfy equation (4.3b). Repeating solving the eigenvalue problem iteratively, we can get λ⁽ⁱ⁾₁ and u⁽ⁱ⁾₁ . Finally, we stop the iterations until the eigenvalue converges. In other words, it means that ku⁽ⁱ⁾₁ − u⁽ⁱ⁺¹⁾₁ k < Tolerance. When the eigenvectors satisfy the condition of convergence, we can get solutions which we need. We suppose the iteration number is k. Hence, we use u^(k)₁ and λ^(k)₁ to denote convergent solutions. After computing above, we obtain a initial solution:

(u1, u2, λ1, γ) = (u^(k)₁ , u2, λ^(k)₁ , 0),

where u^(k)₁ ∈ R^N and u2 is a zero vector ∈ R^N and λ^(k)₁ , 0 ∈ R¹.

• Algorithm Step 1.

u⁽⁰⁾ =rand(N, 1) Let u⁽⁰⁾ = ^√1_c

1u⁽⁰⁾, then C^>u^(0)°2 = 1.

Step 2.

Set i = 0.

Compute (A + α₁u^(i)°2 )u⁽ⁱ⁺¹⁾ = λ_minu⁽ⁱ⁺¹⁾. Let u⁽ⁱ⁺¹⁾ = c₁u⁽ⁱ⁺¹⁾ such that C^>u^(i+1)°2 = 1.

Step 3.

If ku⁽ⁱ⁾− u⁽ⁱ⁺¹⁾k <Tol, then Stop.

Step 4.

Set i = i + 1.

Go to Step 2.

4.2 Standard Continuation Method

In this section, we introduce the process of standard continuation method [5].

We use continuation method approximating to the solution curve. Mainly, we concentrate on solving the system as following

G(u, γ) = 0, (4.4)

Hence, we can find that the Jacobian of G(u, γ) is derived from the G(u, γ). We denote the Jacobian of G(u, γ) as J_G(u, γ) = [G_u, G_γ] ∈ R(2N +1)×(2N +2). We can see that

G_u = the predictions and corrections. We denote (u₀, γ₀)^> as the initial solution, where u0 = (u^>₁₀, u^>₂₀, λ10). And the arc length s is the parameter of the normal vector of the plane is ~p₀. We can find that the point of intersection

of the plane N₁ and solution curve. We denote (u₁, γ₁)^> as the point of intersection. Finally, in order to let (u⁽¹⁾₁ , γ₁⁽¹⁾)^> converge to (u₁, γ₁)^>, we use Newton’s method iteratively. Repeating using the steps above, we can get (u2, γ2)^>, (u3, γ3)^>, · · · . We will show the details of the continuation method as follows.

• predictions

First, we want to find the tangent vector as the predicting vector.

Here, we will use the matrix J_G. We let s be an arc length parameter of solution curve. Therefore, we can rewrite equation (4.4) into the form as following

G(u(s), γ(s)) = 0. (4.5)

Differentiating equation (4.5) with respect to s, we can obtain

G_u˙u(s) + G_γ˙γ(s) = 0. (4.6)

⇒£

G_u G_γ ¤·

˙u(s)

˙γ(s)

¸

= 0. (4.7)

Here, we denote the tangent vector of initial value on solution curve as

~p₀ =

· ˙u₀(s)

˙γ₀(s)

¸ .

Thus, we can obtain

(u⁽¹⁾₁ (s), γ₁⁽¹⁾(s))^> = (u₀(s), γ₀(s))^>+ ¯s~p₀,

=

· u0(s) γ₀(s)

¸ + ¯s

· ˙u0(s)

˙γ₀(s)

¸ .

Next, we consider a plane N₁ whose normal vector is ~p₀. We can see that

N₁ : ˙u^>₀(u − u⁽¹⁾₁ ) + ˙γ₀(γ − γ₁⁽¹⁾) = 0.

⇒ ˙u^>₀u + ˙γ0γ = ˙u^>₀u⁽¹⁾₁ + ˙γ0γ₁⁽¹⁾.

• corrections

Now, we want to find the point of intersection of the plane N1 and the solution curve. We introduce the idea of Newton’s method. By using Newton’s method, we approximate the point of intersection.

– Newton’s method

From Figure 2, we see the graph of a function f : R → R. We can see that f crosses the x-axis at x = c. It shows that c is a root of f (x) = 0, that is, f (c) = 0. We want to approximate root c. First, we choose a point x₁ which is close to c. The point x₂ is the point of intersection of the tangent line at (x₁, f (x₁)) and x-axis. We find x₂ is closer to c than x₁. The tangent line at (x₁, f (x₁)) is y = f (x₁) + f⁰(x₁)(x − x₁). The tangent line intersects the x-axis at y = 0. Therefore, the equation of tangent line will become as follows. 0 = f (x₁) + f⁰(x₁)(x − x₁). Solving the equation, we can get x = x₁− f (x₁)

f⁰(x1), where f⁰(x₁) 6= 0. Here, we let x₂ = x₁ − f (x₁)

f⁰(x1). Next, we use the same steps above to get x₃, x₄, · · · . We can see that the approximations will close to the root c gradually. We denote the tangent line at (xn, f (xn)) as follows.

x_n+1 = x_n− f (x_n)

f⁰(xn), n = 1, 2, · · · .

y

c x x1

x2

( )1

f x

(2) f x

f

x3

Figure 2: Newton’s method

We consider the Newton’s method in Rⁿ, n ≥ 2. We let function f : Rⁿ → Rⁿ. And we choose x which is a solution of f , that is, f (x) = 0. By using the same steps above, we can obtain the formula as follows.

x_n+1 = x_n− f (x_n) · J_f⁻¹(x_n), n = 1, 2, · · · , where J_f is the Jacobian matrix of f .

Now, we have the point (u⁽¹⁾₁ , γ₁⁽¹⁾)^>. We want to get the points which

And we denote the Jacobian matrix of F as J_F =

· G_u G_γ

˙u^>₀ ˙γ0

¸ .

Then, we use the formula of Newton’s method, that is, we solve the system as follows.

Figure 3: predict and correct (Newton’s method)

4.3 Quotient Transformation Invariant Continuation Method

In this subsection, we will introduce the quotient transformation invariant continuation method [7]. And we replace standard continuation method with the quotient transformation invariant continuation method. In the process of computing, some problems will occur. We need to overcome these problems.

First, we will face the problem about the predicting directions. Before discussing the problem, we consider the solution curve of standard continu-ation method. From equcontinu-ation (4.4), we know that there are four variables and three equations. Consequently, we can get the solution which is a two-dimensional manifold. We denote the solution set of equation (4.4) as M.

Hence, we suppose that the solution will be a pipeline.

M

Figure 4: standard continuation method

From Figure 4, we can see that why the standard continuation method will fail. We will find that the predicting directions will move along the manifold. As a result, we can not get meaningful predictions of standard continuation method. Hence, we want to use quotient transformation in-variant continuation method which can improve the disadvantage. For the purpose of letting the solution of equation (4.4) can be a straight line, we add a plane to equation (4.4). Thus, there are four variables and four equations.

M

C

'( ) u

prediction

P

( ) u

θ

θ

Figure 5: quotient transformation invariant continuation method

From Figure 5, we can see that the direction of predictions will be on the plane. And the quotient solution curve C is on the plane P . We will explain the details of the process of transformation.

From equation (4.5), we can rewrite into the forms as follows.

G(y(s)) = 0,

where y(s) = (u(s), γ(s)). And we differentiate G(y(s)) = 0 with respect to s. Thus, we can get

JG(y(s)) ˙y(s) = 0, (4.8)

where J_G(y(s)) is the Jacobian matrix of G. And we denote a tangent vector of C at y(s) as ˙y(s) = ( ˙u(s)^>, ˙γ(s))^>. We let C be the quotient solution curve.

We show that the conception of quotient solution. Let u(θ) = u · e^iθ,

= (u₁+ iu₂) · e^iθ,

= u₁cos θ + iu₁sin θ + iu₂cos θ − u₂sin θ,

= (u₁cos θ − u₂sin θ) + i(u₁sin θ + u₂cos θ).

We define

u(θ) =

· u1cos θ − u2sin θ u₁sin θ + u₂cos θ

¸ .

Let equation (2.3) be multiplied in both right and left sides by e^iθ. We can obtain the equation as follows.

[−∆u + iγ∂_θ1u + V u + α₁|u|²u] · e^iθ = λ₁u · e^iθ.

We know that u will satisfy equation (2.3) if u is a solution. And we suppose that u · e^iθ will be a solution, that is,

[−∆ + iγ∂_θ1 + V + α₁|u · e^iθ|²]u · e^iθ = λ₁u · e^iθ. (4.9) We put |u · e^iθ|² = |u|²|e^iθ|² = |u|² into equation (4.9). Then, we can get

[−∆ + iγ∂_θ1 + V + α₁|u|²]u · e^iθ = λ₁u · e^iθ,

⇒[−∆ + iγ∂_θ1 + V + α₁|u|²]u = λ₁u.

And we consider the normalized constraint.

Z

|u · e^iθ|²dx = Z

|u|²dx = 1.

Hence, we can verify the assumption, that is, if u is a solution, then ue^iθ is also a solution. We can see that the solution M contains transformation invariant solutions u(θ).

We denote the quotient solution as C. Now, we want to get the tangent vector of u(θ). Hence, we will differentiate u(θ) with respect to θ. We can get

u⁰(θ) =

· −u₁sin θ − u₂cos θ u₁cos θ − u₂sin θ

¸

. (4.10)

We choose θ = 0 and put it into equation (4.10). Then we can get u⁰(0) =

∂u

∂θ(0).

∂u

∂θ(0) =

· −u₂ u₁

¸

= (−u^>₂, u^>₁)^>.

Consequently, we can obtain the normal vector of the plane which we add, that is, the tangent vector of u(θ). We denote (^∂u_∂θ(0)^>, 0)^> as the tangent vector.

Then, we begin the further exploration of the quotient transformation invariant continuation method. First, we discuss the part of prediction. We know that the prediction vector ˙y(s) = ( ˙u(s)^>, ˙γ(s))^>. And ( ˙u(s)^>)^> is perpendicular to the tangent vector (^∂u_∂θ(0)^>, 0)^>. Next, we consider the part of correction. We can see that the vector of correction will be on the plane which we add. And the normal vector of the plane is (^∂u_∂θ(0)^>, 0, 0)^>.

• predictions

By using fixed point iteration method, we can get the initial solution (u₀, γ₀)^> ∈ R^{(2N +2)×1}. We let y₀ = (u₀, γ₀)^>. We discuss how we find the vectors of predictions. In the part of predictions, we know that the vector of predictions will be on the plane P . The plane has the normal vector (^∂u_∂θ(0)^>, 0, 0)^>. The vector of predictions ˙y(s) = ( ˙u(s)^>, ˙γ(s))^>

should satisfy the equation (4.8). And ( ˙u(s)^>)^> is perpendicular to the tangent vector (^∂u_∂θ(0)^>, 0)^>. Therefore, we add a plane as follows.

a^>_θ ˙u = 0, (4.11)

where

a_θ = (∂u

∂θ(0)^>, 0)^>

= (((−u^>₂, u^>₁)^>)^>, 0)^>

= (−u^>₂, u^>₁, 0)^>.

Originally, we solve the system of (4.7). We let y(s) = (u(s), γ(s)) and ˙y(s) = ( ˙u(s)^>, ˙γ(s))^>. Hence, we can rewrite it into the form as

follows. £

G_u G_γ ¤ ˙y(s) = 0.

Since we add a plane which is satisfied equation (4.11), we have to solve the system as follows.

· Gu Gγ

a^>_θ 0

¸

˙y =

· 0 0

¸ .

By computing above, we can get the vectors of predictions. Conse-quently, we will get the new solution (u⁽¹⁾₁ , γ₁⁽¹⁾)^> = (u₀, γ₀)^>+ ¯s · ˙y, where k ˙yk₂ = 1.

In the process of predictions, we let y_i = (u_i, γ_i)^> be the approximate solution every time. And we have to solve the system as follows.

· G_u G_γ a^>_θ 0

¸

˙y_i =

· 0 0

¸

. (4.12)

Then, we denote the vector of prediction as ˙y_i = ( ˙u^>_i , ˙γ_i)^>. As a result, we will get the new solution (u⁽¹⁾_i , γ_i⁽¹⁾)^>= (u_i−1, γ_i−1)^>+ ¯s_i· ˙y_i, where i = 1, 2, · · · , ¯s_i is the length of step every time and the two norm of ˙y_iis equivalent to one. After computing the vector of predictions, we hope that the new solution will converge to the solution by using Newton’s method. Hence, we introduce the steps of corrections.

• corrections

We continue discussing the correction part of quotient transforma-tion invariant continuatransforma-tion method. Now, we have the new solutransforma-tion (u⁽¹⁾_i , γ_i⁽¹⁾)^> which is on the plane P_i. Here, P_i represents the plane which has the normal vector (^∂u_∂θ(0)^>, 0)^>. And we know that the vec-tor of corrections will also be on the plane P_i. The solutions which will proceed to converge are also on the plane. We want to get the vectors of correction. As a result, we have to consider the system as follows.





G : G(u, γ) = 0, P : a^>_θ ˙u = 0,

N_i : ˙u^>_i−1u + ˙γ_i−1γ = ˙u^>_i−1u⁽¹⁾_i + ˙γ_i−1γ_i⁽¹⁾

It means that the solution curve is related to the system above. We know that the point which we want to find has to satisfy the system,

that is, (u_i, γ_i)^>. We consider the system as follows.

H =





G : G(u_i, γ_i) = 0, P : a^>_θ ˙u_i−1= 0,

N_i : ˙u^>_i−1u_i+ ˙γ_i−1γ_i = ˙u^>_i−1u⁽¹⁾_i + ˙γ_i−1γ_i⁽¹⁾

For the purpose of getting the vectors of corrections, we still use New-ton’s method. Therefore, we consider the Jacobian matrix of H, that is,

J_H =



 G_u G_γ a^>_θ 0

˙u^>_i−1 ˙γ_i−1



 ∈ R(2N +3)×(2N +2). (4.13)

Then, we use Newton’s method to solve it iteratively. And we put it into the formula as follows.

(u⁽ⁿ⁺¹⁾_i , γ_i⁽ⁿ⁺¹⁾) = (u⁽ⁿ⁾_i , γ_i⁽ⁿ⁾) − H(u⁽ⁿ⁾_i , γ_i⁽ⁿ⁾) · J_H⁻¹(u⁽ⁿ⁾_i , γ_i⁽ⁿ⁾), where n = 1, 2, · · · . Hence, we can get the solutions which will converge gradually. We stop the steps of correction until the solutions converge.

• improvement in computing

Now, we have the system H and J_H. In the process of computing, we observe the last row of G_u. We will find that the norm of the last row is much smaller than other rows. It will result in error in calculation.

The last row will be ignored. To review the system H which we want to solve, we observe G₃, that is, the normalized constraint. We let G3 be multiplied in both sides by accuracy. And we let accuracy = max(k(u, γ)k, 1). We also let the planes P and N_i be multiplied in both sides by accuracy. Now, we have new normalized constraint, in other words, we rewrite G3 to accuracy · (C^>U⁰ − 1). And we rewrite the planes into the forms as following.

P : accuracy · (a^>_θ ˙u_i−1) = 0,

N_i : accuracy · ( ˙u^>_i−1u_i+ ˙γ_i−1γ_i) = accuracy · ( ˙u^>_i−1u⁽¹⁾_i + ˙γ_i−1γ_i⁽¹⁾).

By using the skill above, we can obtain H and J_H accurately. And we have to rewrite equation (4.5) as follows.

G_u =





A + α¯ ₁(3u^°₁² + u^°₂² ) 2α₁(u₁◦ u₂) − γS −u₁ 2α1(u1 ◦ u2) + γS A + α¯ 1(3u^°₂² + u^°₁² ) −u2

accuracy · (2C^>◦ u^>₁) accuracy · (2C^>◦ u^>₂) 0



 ∈ R(2N +1)×(2N +1).

In the part of prediction, system (4.12) has to rewrite into the form as

Furthermore, we consider the part of prediction. We rewrite equation (4.13) into the form as follows.

J_H =

In this section, we will introduce the algorithm of continuation method. We can see that the algorithm consists of predictions and corrections. We show the Algorithm as follows.

Step1. Computing initial value Compute initial value: (u₀, γ₀).

Set the length of step.

Set count = the number of possible solutions.

Step2. Predictions for i = 0, 1, · · · ,count

Compute Jacob_F =

· G_u(u_i, γ_i) G_γ(u_i, γ_i) accuracy · a^>_θ 0

¸ .

Find the desired singular values and singular vectors of Jacob_F.

Let the smallest singular vector be the vector of predictions: P red uγ.

end

Step4. Corrections (Newton’s method) Compute Del uγ = −J_H⁻¹· H(u⁽¹⁾_i+1, γ_i+1⁽¹⁾)

(u⁽²⁾_i+1, γ_i+1⁽²⁾) = (u⁽¹⁾_i+1, γ_i+1⁽¹⁾) + Del uγ. ComputeJ_H, N, G, H

If(kH(u⁽¹⁾_i+1, γ_i+1⁽¹⁾)k/kH(u⁽²⁾_i+1, γ_i+1⁽²⁾)k < 10 and kH(u⁽²⁾_i+1, γ_i+1⁽²⁾)k > 10⁻⁸·accuracy

and kDel uγk > 10⁻⁷)

step = step/2 and go to Step3 else

go to Step5 end

Step5. Corrections (Newton’s method) for j = 2, 3, · · · ,

If (kH(u^(j)_i+1, γ_i+1^(j))k > 10⁻⁸· accuracy) Del uγ = −J_H⁻¹· H(u^(j)_i+1, γ_i+1^(j)) (u^j+1_i+1, γ_i+1^j+1) = (u^j_i+1, γ_i+1^j ) + Del uγ.

ComputeJ_H, N, G, H

if(kH(u^(j)_i+1, γ_i+1^(j))k/kH(u^(j+1)_i+1 , γ_i+1^(j+1))k < 5) if(kDel uγk < 10⁻⁷)

stop else

step = step/2 and go to Step 2.

end end else

(u^(j)_i+1, γ_i+1^(j)) → (u_i+1, γ_i+1) and go to Step6.

end end

Step6. Save

Compute accuracy=p

kJHkk(ui+1, γi+1)k if accuracy < 1

accuracy=1 end

Go to Step 2.

5 Numerical Results

We use quotient transformation invariant continuation method to approxi-mate the solution curve. In this section, we will show that the numerical results as follows. First, we can see that the approximative solution curve in Figure 6 and Figure 7. In Figure 6 and Figure 7, γ represents the value of x-axis. And the values of y-axis are p

u²₁+ u²₂,u₁,u₂ separately. Next, to change the value α₁, we observe the changes of energy. From Figure 8 and Figure 9, it show that the energy E(φ) will increase as α₁ increases. Here,

we let α₁ = 0, 1, 5, 10. Finally, we show that the difference between the eigenvalue λ₁ and the energy functional E(φ) is ^α₂¹ R

|u|⁴. We compute the difference, that is, λ₁ − E(φ) − ^α₂¹ R

|u|⁴. And we obtain the value is 10⁻⁷ approximately. Hence we can verify that the relation between the eigenvalue and the energy functional.

γ

0 0.099 8.62010 17.63620 28.367 30 35.459

2 2

1 2

u+u

Figure 6: Case: α₁ = 5, N = 50 × 50, the value of y-axis: p

u²₁+ u²₂.

γ

00.099 8.62010 17.63620 28.367 30 35.459

u1

γ

00.099 8.62010 17.63620 28.367 30 35.459

u2

Figure 7: Case: α₁ = 5, N = 50 × 50, the value of y-axis: left graph:u₁, right graph: u₂.

0 5 10 15 20 25 30 35

We discuss the process of discretization on a nonlinear Schr¨odinger equation and introduce the standard continuation method. In the process, we face on the difficult computing problem. Hence, we improve the method by using the quotient invariant transformation continuation method. In order to find the solution curve successfully, we add the plane. Hence, we can see the diagrams of approximate solution curve in section 5. We also verify the relation between the eigenvalue and the energy functional.

References

[1] Chia-Feng Tsao. Numerical study of a two-component Bose-Einstein condensate in magnetic field. preprint, 2007.

[2] Volodymyr Sushch On some discrete model of the magnetic Laplacian.

preprint, 2002.

[3] I.Schindler, K.Tintarev A nonlinear Schr¨odinger equation with external magnetic field. preprint.

[4] Ambar K. Mitra Finite Difference Method for the Solution of Laplace Equation. preprint.

[5] H.B. Keller. Lectures on Numerical Methods In Bifurcation problems.

Springer-Verlag, 1987.

[6] Yuen-Cheng Kuo, Wen-Wei Lin, Shih-Feng Shieh, Weichung Wang. A Hyperplane-Constrained ContinuationMethod for Bound States of Cou-pled Nonlinear Schr¨odinger Equations. (Submitted), 2007.

[7] Yuen-Cheng Kuo , Wen-Wei Lin , Shih-Feng Shieh, and Weichung Wang.

Exploring Bistability in Rotating Bose-Einstein Condensates by a Quo-tient Transformation Invariant Continuation Method. (Submitted) [8] Shu-Ming Chang, Yuen-Cheng Kuo, Wen-Wei Lin, Shih-Feng Shieh.

A Continuation BSOR-Lanczos-Galerkin Method for Positive Bound States of A Multi-Component Bose-Einstein Condensate . J. Comput.

Phys., 210(2):439-458, 2005.

[9] W. Z. Bao. Ground states and dynamics of multi-component Bose-Einstein condensates. SIAM Multiscale Modeling and Simulation,

[9] W. Z. Bao. Ground states and dynamics of multi-component Bose-Einstein condensates. SIAM Multiscale Modeling and Simulation,

在文檔中 National University of Kaohsiung Repository System:Item 310360000Q/10839 (頁 9-0)