Efficient simulation of Schrödinger equation with piecewise constant positive potential
EEFFICIENT SIMULATION OF SCHR ¨ODINGER EQUATIONWITH PIECEWISE CONSTANT POSITIVE POTENTIAL
XUXIN YANG, ANTTI RASILA, AND TOMMI SOTTINEN
Abstract.
In this paper we introduce a new method for the simulationof a weak solution of the Schr¨odinger-type equation where the potentialis piecewise constant and positive. The method, called killing walk onspheres algorithm, combines the classical walk of spheres algorithm withkilling that can be determined by using panharmonic measures. Introduction
We investigate the simulation of the following special case of the time-independent Schr¨odinger equation(1.1) 12 ∆ u ( x ) − q ( x ) u ( x ) = 0 , where D ⊂ R n , n ≥ u : D → R is a piecewise twotimes continuously differentiable function and q is of the form(1.2) q ( x ) = M (cid:88) m =1 λ j D m ( x ) , where, for each m = 1 , . . . , M , λ m is a non-negative constant, D m ⊂ D isa domain such that D m ∩ D (cid:96) = ∅ for m (cid:54) = (cid:96) , and D = D ∩ (cid:0) M (cid:91) m =1 D m (cid:1) . This class of functions q given by (1.2) is motivated by the phenomenoncalled quantum tunneling in particle physics. Also, the form (1.2) can bemotivated by being an approximation of general positive q . Mathematics Subject Classification.
Key words and phrases.
Brownian motion; killing walk on spheres; harmonic measure;numerical algorithm; Yukawa equation; Schr¨odinger equation;This work is supported by the NNSF of China (No. 11571088), Hunan ProvincialNatural Science Foundation of China (14JJ7083), a Key Project Supported by ScientificResearch Fund of Hunan Provincial Education Department (14A028), Scientific ResearchFund of Hunan Provincial Education Department (14C0253), the Aid Program for Scienceand Technology Innovative Research Team in Higher Educational Institutions of HunanProvince. A. Rasila was partially supported by Academy of Finland (No. 289576). T.Sottinen was partially funded by the Finnish Cultural Foundation (National Foundations’Professor Pool). This work was done during T. Sottinen’s research visit to Hunan FirstNormal University, Changsha, PRC. T. Sottinen wishes to thank them for their hospitality. a r X i v : . [ m a t h . NA ] D ec XUXIN YANG, ANTTI RASILA, AND TOMMI SOTTINEN
Since we are considering the case where q is discontinuous, we cannotobtain classical u ∈ C ( D ) solution. Instead, we consider weak solutions,i.e., u is a solution of (1.1) if (cid:90) D u ( x ) 12 ∆ ϕ ( x ) d x = (cid:90) D q ( x ) u ( x ) ϕ ( x ) d x for all ϕ ∈ C ∞ c ( D ). 2. Panharmonic measures
Let W = ( W t ) t ≥ be standard n -dimensional Brownian motion. Let P x be the probability measure such that, P x -almost surely, W = x . Let E x be the expectation associated with the measure P x . Let τ D = inf { t > W t (cid:54)∈ D } , i.e., τ D is the first exit time of the Brownian motion from the domain D .We assume that the domain D , and all the sub-domains D m , are (Wiener)regular, i.e., P y [ τ D = 0] = 1 for all y ∈ ∂D. Then, by [2, Theorem 4.7], a bounded weak solution to the boundary valueproblem(2.1) (cid:26) ∆ u ( x ) − q ( x ) u ( x ) = 0 on D,u ( x ) = f ( x ) on ∂D, where f is a bounded Borel function, is given by(2.2) u ( x ) = E x (cid:104) e − (cid:82) τD q ( W s ) d s f ( W τ D ) (cid:105) . Moreover, if f is continuous, then (2.2) is the unique solution satisfying u ∈ C ( ¯ D ).In the case where q is of the form (1.2) the formula (2.2) takes the form(2.3) u ( x ) = E x (cid:104) e − (cid:80) Mm =1 λ m T m f ( W τ D ) (cid:105) , where T m is the total time the Brownian motion spends in the sub-domain D m before exiting the domain D at time τ D . Since the Brownian motionhas the strong Markov property and the “discounting term” e − λ m T m canbe considered as exponential killing, we can analyze the terms in the sumindependently. This means that on the stochastic set { t ≥ W t ∈ D m } we can consider the corresponding Yukawa equation, or Brownian motionwith exponential killing. This can be analyzed by using the panharmonicmeasures as in [8]. Now we recall some facts of panharmonic measures.The harmonic kernel is(2.4) h x ( D ; d y, t ) = P x [ W τ D ∈ d y | τ D = t ] d P x d t [ τ D ≤ t ] . In [8] it was shown that the harmonic kernel exists for all regular domains D and that the harmonic measure can be written as H x ( D ; d y ) = (cid:90) ∞ t =0 h x ( D ; d y, t ) d t. IMULATION OF SCHR ¨ODINGER EQUATION 3
The λ -panharmonic measure is(2.5) H xλ ( D ; d y ) = (cid:90) ∞ t =0 e − λt h x ( D ; d y, t ) d t. It was shown in [8] that all bounded solutions to boundary value problemto the Yukawa equation, i.e. (2.1) with constant q ( x ) ≡ λ >
0, can berepresented via the λ -panharmonic measure as(2.6) u ( x ) = (cid:90) y ∈ ∂D f ( y ) H xλ ( D ; d y ) . It was also shown in [8] that all the panharmonic measures H xλ ( D ; · ), λ ≥ Z xλ ( D ; · ) such that H xλ ( D ; d y ) = Z xλ ( D ; y ) H x ( D ; d y ) . The harmonic measure corresponding to λ = 0 is a probability measure.Indeed, H x ( D ; d y ) = P x [ W τ D ∈ d y ] . The panharmonic measure corresponding to case λ > Z xλ ( D ; y ) = E x (cid:104) e − λτ D (cid:12)(cid:12)(cid:12) W τ D = y (cid:105) ≤ . This Radon–Nikodym derivative and the classical harmonic measure play acentral part in constructing the killing walk on spheres (KWOS) simulationalgorithm in the next section. The algorithm is an extension of the classicalwalk of spheres algorithm (WOS) due to Muller [7].3.
Construction of the algorithm
The general formula (2.2) for the Schr¨odinger equation provides an obvi-ous way to solve the boundary value problem (2.1) by simulating Brownianparticles. Indeed, let w kx,δt , k = 1 , . . . , K be independent simulated Brown-ian trajectories starting from x ∈ D with time-mesh δt . So, each w kx,δt is anindependent varying-length vector ( w kx,δt (0) , w kx,δt (1) , . . . , w kx,δt ( T kD )), where w kx,δt (0) = x , w kx,δt ( j ) = w kx,δt ( j −
1) + √ δt ξ kj ,ξ kj ’s are independent standard n -dimensional normal random variables, and T kD is the first index j such that w kx,δt ( j ) (cid:54)∈ D . Set(3.1) ˆ u K,δt ( x ) = 1 K K (cid:88) k =1 exp − T k (cid:88) j =1 q (cid:16) w kx,δt ( j ) (cid:17) δt f (cid:16) w kx,δt ( T k ) (cid:17) . Then, by the strong law of large numbers and by the dominated convergencetheorem, u K,δt ( x ) → u ( x ) for all x ∈ D with P x -probability 1 as δt → K → ∞ .The obvious problem in the general simulation method described aboveis that it is computationally very heavy. However, since the solution of XUXIN YANG, ANTTI RASILA, AND TOMMI SOTTINEN the Schr¨odinger equation (2.2) depends on the entire path of the Brow-nian motion, there seems to be no other way than take the δt to be verysmall. Consequently, finding an efficient simulation algorithm for the generalSchr¨odinger equation by using the Brownian motion seems a very challeng-ing problem, indeed.Let us then consider the special case of the Schr¨odinger equation, where q is of the form (1.2). This means that we have the Yukawa equation withdifferent parameters in different regions. In this case, we can sometimesavoid simulating the Brownian particles in a fine time-mesh by killing theparticle in the region with exponential clock.The idea of the simulation is comes from the following representationof the solutions, which follows from [8, Corollary 2.11], the strong Markovproperty of the Brownian motion and the memoryless property of the expo-nential distribution.3.2. Theorem.
Let q be given by (1.2) . Then the function u in (2.2) canbe give as (3.3) u ( x ) = E x [ f ( W τ D ); τ ∗ > τ D ] , where τ ∗ is an independent exponential clock with intensity λ m on sub-domain D m . Note that the formula (3.3) depends on the entire path for the Brownianmotion only through the killing time τ ∗ .Our estimator for u ( x ) is(3.4) ˆ u K ( x ) = 1 K (cid:88) k ∈ K ∗ ( λ ) f ( w kx ( τ k )) , where w kx , k = 1 , . . . , K are independent simulations of the trajectoriesBrownian particles starting from point x , and the set K ∗ ( λ ) ⊂ { , . . . , K } contains the particles that are not killed; τ k is the termination-step of thealgorithm.The optimistic idea of how to individual trajectories w x = w kx are gener-ated as follows: Set w x (0) = 0 and suppose w x (1) , . . . , w x ( j −
1) are gener-ated. Suppose w x ( j − ∈ D m , and not is near the boundary. Suppose the λ m -panharmonic Radon–Nikodym derivative and the harmonic measure forthe domain D m are known (this is the optimistic part). Generate w x ( j )on the boundary by using the harmonic measure H w x ( j − λ m ( D ; d y ). Now w j ( x ) = y ∈ ∂D m . Kill the particle with probability 1 − Z w x ( j − λ m ( D m ; y ).If the particle dies, the algorithm terminates and the particle adds 0 to thesum (3.4). If the particle survives and y ∈ ∂D , then the algorithm termi-nates and f ( w kx ( j )) is added to the sum (3.4). If y ∈ ∂D m (cid:48) \ ∂D , thenwe have to simulate the particle with fine time-mesh in order to allow theparticle to enter well inside the domain D m (cid:48) . Once the particle is well insidea domain D m (cid:48) , we can generate a new particle by using the correspondingharmonic measure and the λ m (cid:48) -panharmonic Radon–Nikodym derivative asbefore. IMULATION OF SCHR ¨ODINGER EQUATION 5
The harmonic measures and panharmonic Radon–Nikodym derivativesare not known for arbitrary domain. This limits the applicability of the ideaabove. However, for balls B ( x, r ) ⊂ R n the Radon–Nikodym derivative Z xλ ( B ( x, r ); y ), λ > y and by the self-similarity of the Brownian motion Z xλ ( B ( x, r ); y ) = ψ ( µr ) , where µ = √ λ. The function ψ is the Laplace transform of the first hittingtime of a Bessel process on the level 1 starting from 0 . This is well-known,see e.g. [9, Theorem 2] for n ≥ ψ ( µ ) = µ ν ν Γ( ν + 1) I ν ( µ ) , µ > , where I ν ( x ) = ∞ (cid:88) m =0 m !Γ( m + ν + 1) (cid:16) x (cid:17) m + ν is the modified Bessel function of the first kind of order ν = n/ −
1. Theharmonic measure for balls is, due to the rotational symmetry of Brownianmotion, simply the uniform distribution.This leads to the following simulation algorithm:3.6.
Algorithm (Killing walk on spheres (KWOS)) . (i) Set w x (0) = x and suppose w x (1) , . . . , w x ( j −
1) are generated.(ii) Generate w x ( j ) as follows:(a) If w x ( j − ∈ D m and if dist( w x ( j − , ∂D m ) < ε startgenerating the Brownian path with fine time-mesh δt (cid:28) ε :generate w x ( j ) = w x ( j −
1) + √ δtξ ( j ) , where ξ ( j ) an n -vector of independent standard normal ran-dom variables. Generate exponential killing: with probability p ( δt ) = 1 − e − λ δt . If killing occurs, the algorithm terminatesand adds 0 to the sum (3.4).(b) If w x ( j − ∈ D m and if dist( w x ( j − , ∂D m ) ≥ ε generate w x ( j ) ∈ ∂B ( w x ( j − , r ), r = dist( w x ( j − , ∂D m ) uniformly.Kill the particle with probability 1 − ψ ( µr ), where ψ ( µr ) isgiven by (3.5). If killing occurs, the algorithm terminates andadds 0 to the sum (3.4).(c) If dist( w x ( j − , ∂D ) < ε let w x ( j ) be the projection of w x ( j −
1) to the boundary ∂D . In this case the algorithmterminates and adds f ( w x ( j )) to the sum (3.4).4. Examples
In this section, we present examples to motivate our algorithm and to illus-trate its potential applications in one and two-dimensional settings. Theseexamples were computed by using a simple implementation algorithms on
XUXIN YANG, ANTTI RASILA, AND TOMMI SOTTINEN
Mathematica, and they mainly chosen from the point of view of visualisa-tion. It should be noted that our approach is obviously more much attractivein higher dimensions, where many other methods for solving such are notavailable, or lead into excessive computation times . One-dimensional example.
Next we consider the equation (1.1) in thecase where u : [ a, b ] → R and [ a, b ] is a real interval, i.e., ∆ u = u (cid:48)(cid:48) , and q ( x )is as in (1.2). Then we obtain the impulsive differential equation(4.1) 12 u (cid:48)(cid:48) = λ ( x ) u, λ ( x ) = λ j for x ∈ ( x j − , x j ) , j = 2 , . . . , M, where a = x < x < · · · < x M − < x M = b , M ≥ u ∈ C ([ a, b ]),and u (cid:48)(cid:48) is piecewise continuous on [ a, b ].One should note that the theory of impulsive differential equations ismuch richer than the corresponding theory of differential equations withoutimpulse effects. For example, initial value problems of such equations maynot, in general, possess any solutions at all even when the correspondingdifferential equation is smooth enough, fundamental properties such as con-tinuous dependence relative to initial data may be violated, and qualitativeproperties like stability may need a suitable new interpretation. Moreover,a simple impulsive differential equation may exhibit several new phenom-ena such as rhythmical beating, merging of solutions, and noncontinuabilityof solutions. See e.g. [6] for more information about impulsive differentialequations. However, these problems do not arise in the case of (4.1) as,the boundary value problem (2.1) has a continuous weak solution, and itcoincides piecewise with the classical solution.The following algorithm is for the special case (4.1) with a = x < x Algorithm (Gambler’s ruin with killing walk on spheres (GR-KWOS)) . (i) Set w (0) = x and suppose w (1) , . . . , w ( j − 1) are generated.(ii) Generate w ( j ) as follows:(a) If w ( j − ∈ ( x , x ) then it will exit the domain in x withprobability ( x − w ( j − / ( x − x ). The algorithm terminatesand gives u ( x ) to the sum (3.4). Otherwise w ( j ) = x andwe are entering the killing zone.(b) If w ( j − ∈ ( x , x ) then it will exit the domain in x withprobability 1 − ( x − w ( j − / ( x − x ). The algorithm termi-nates and gives u ( x ) to the sum (3.4). Otherwise w ( j ) = x an we are entering the killing zone. Time required for computing the examples in this section with a MacBook Air com-puter varied between seconds and minutes. Wolfram Mathematica 10.2 and very ba-sic implementations of the algorithms with no performance optimisations were used incomputations. IMULATION OF SCHR ¨ODINGER EQUATION 7 (c) If w ( j − ∈ [ x , x ] we are entering or in the killing zone. Inthis case use KWOS Algorithm 3.6.The above algorithm can further be refined by using the harmonic measureand the λ m -harmonic Radon–Nikodym derivative in all the domains D m , m = 1 , , ψ of (3.5) can beobtained from [1, (3.0.1)]: E x [ e − λT ] = cosh(( b + a − x ) (cid:112) λ/ b − a ) (cid:112) λ/ , where T = inf { t > W t (cid:54)∈ ( a, b ) } and W is the one-dimensional Brownianmotion by using parameters x = 0, a = − r , b = r , and r > Example. Let x = 0 , x = 1 , x = 2 , x = 3 . 5, and consider theboundary value problem for the equation (4.1) with u ( x ) = 1 , u ( x ) = 2, λ = λ = 0 and λ = 2. Then the solution of the problem is given by u ( x ) ≈ − . t for x ∈ (0 , , . e − t (274 . 757 + e t ) for x ∈ (1 , , − . . t for x ∈ (2 , . . The solution u , and approximate solution obtained through GR-KWOSAlgorithm 4.2, are illustrated in Figure 1. Figure 1. The solution u of the BVP of Example 4.3(dashed) and its approximation ˆ u obtained by using GR-KWOS algorithm. Two-dimensional examples. Next we consider examples where the twodimensional potential function of the equation (1.1) is computed by usingAlgorithm 3.6. We start by computing the pure Yukawa case. Potentialtheory of the Yukawa equation has been studied by Duffin [5]. Stochasticmethodology for studying Yukawa potentials was developed in [8].4.4. Example. (Yukawa equation) Let D = (cid:8) ( x, y ); 0 < x < min { , y + 1 } , < y < (cid:9) . XUXIN YANG, ANTTI RASILA, AND TOMMI SOTTINEN We solve the Yukawa equation, i.e., the equation (1.1) with constant q ≡ f ( x, y ) = x + y . The approximatesolution ˆ u of (3.4), where K = 1000, and a chain of balls (disks) used inWOS style simulation of the path of one Brownian particle, are illustratedin Figure 2. Figure 2. An approximation of the solution ˆ u of theBVP of Example 4.4 with KWOS algorithm (left), and chainof disks used by KWOS algorithm in simulation of an indi-vidual path (right).4.5. Example. (Mixed Laplace–Yukawa equation) Consider the problem(1.1) on the domain D = (cid:8) ( x, y ); 0 < x < . , < y < . (cid:9) , where q is as in (1.2), the boundary values are given by f ( x, y ) = x , D = (cid:8) ( x, y ); 0 < x < y, < y < (cid:9) ,D = (cid:8) ( x, y ); 0 < x < . , x < y < . (cid:9) ,λ = 0 and λ = 2 . The approximate solution ˆ u of (3.4), is illustratedin Figure 3. 5. Conclusion As a conclusion, the killing walk on spheres (KWOS) algorithm is a verysimple tool to compute a solution to the Schr¨odinger equation with piecewiseconstant positive potential. The algorithm is based on the classical walkon spheres. It avoids estimating the exit time from the ball by using theinterpretation of killing. If the potential has negative values, then the killinginterpretation is no longer available and one needs to estimate the exit time IMULATION OF SCHR ¨ODINGER EQUATION 9 Figure 3. The approximate solution ˆ u of the BVP ofExample 4.5 with KWOS algorithm.by using, e.g., the recent walk on moving spheres (WOMS) algorithm dueto Deaconu et al. [3, 4]. References [1] A. Borodin and P. Salminen : Handbook of Brownian Motion; Facts and Formulae ,2nd edition, Birkh¨auser, 2002.[2] K.L. Chung and Z. Zhao : From Brownian motion to Schr¨odinger’s equation , 2ndPrinting, Springer, 2001.[3] M. Deaconu and S. Herrmann : Hitting Time for Bessel Processes — Walk onMoving Spheres Algorithm (WOMS), Ann. Appl. Probab. , No. 6 (2013), 2259–2289.[4] M. Deaconu , S. Herrmann and S. Maire : The walk on moving spheres: A newtool for simulating Brownian motion’s exit time from a domain, Mathematics andComputers in Simulation (2015), article in press.[5] R.J. Duffin : Yukawan potential theory, J. Math. Anal. Appl. (1971), 105–130.[6] V. Lakshmikantham , D. Bainov , and P. Simeonov : Theory of impulsive differen-tial equations , World Scientific, Singapore, 1989.[7] M.E. Muller : Some continuous Monte Carlo methods for the Dirichlet problem. Ann. Math. Statist. (1956), 569–589.[8] A. Rasila and T. Sottinen : Yukawa Potential, Panharmonic Measure and Brown-ian Motion. Preprint. arXiv:1310.2167[9] J.G. Wendel : Hitting Spheres with Brownian Motion, Ann. Probab. (1), (1980),164–169.(1), (1980),164–169.