WKB-based scheme with adaptive step size control for the Schrödinger equation in the highly oscillatory regime
aa r X i v : . [ m a t h . NA ] F e b WKB-based scheme with adaptive step size control forthe Schr¨odinger equation in the highly oscillatory regime
Jannis K¨orner, Anton Arnold ∗ , Kirian D¨opfner Institute of Analysis and Scientific Computing, Technische Universit¨at Wien, WiednerHauptstr. 8-10, 1040 Wien, Austria
Abstract
This paper is concerned with an efficient numerical method for solving the1D stationary Schr¨odinger equation in the highly oscillatory regime. Be-ing a hybrid, analytical-numerical approach it does not have to resolve eachoscillation, in contrast to standard schemes for ODEs. We build upon theWKB-based marching method from [1] and extend it in two ways: By com-paring the O ( h ) and O ( h ) methods from [1] we design an adaptive stepsize controller for the WKB method. While this WKB method is very ef-ficient in the highly oscillatory regime, it cannot be used close to turningpoints. Hence, we introduce for such regions an automated method cou-pling, choosing between the WKB method for the oscillatory region and astandard Runge-Kutta-Fehlberg 4(5) method in smooth regions.A similar approach was proposed recently in [2, 3], however, only for a O ( h )-method. Hence, we compare our new strategy to their method on twoexamples (Airy function on the spatial interval [0 , ] with one turning pointat x = 0 and a parabolic cylinder function having two turning points), andillustrate the advantages of the new approach w.r.t. accuracy and efficiency. Keywords:
Schr¨odinger equation, highly oscillatory wave functions, higherorder WKB-approximation, adaptive step size control, Airy function,parabolic cylinder function ∗ Corresponding author.
Email addresses: [email protected] (Jannis K¨orner), [email protected] (Anton Arnold), [email protected] (KirianD¨opfner)
Preprint submitted to Computer Physics Communications February 8, 2021 . Introduction
This paper deals with the numerical solution of the highly oscillatory 1DSchr¨odinger equation ε ϕ ′′ ( x ) + a ( x ) ϕ ( x ) = 0 , x ∈ R . (1.1)Here, 0 < ε ≪ ε := ~ √ m ) and ϕ the (pos-sibly complex valued) Schr¨odinger wave function. The real valued coefficientfunction a ( x ) := E − V ( x ) is related to the potential V ( x ), and E denotesthe given energy of a considered particle. We shall assume here that a isbounded away from zero, i.e. a ( x ) ≥ τ for some τ >
0. The (local) wavelength of a solution ϕ to (1.1) is given by λ ( x ) = πε √ a ( x ) . Hence, for smallvalues of ε the solution is highly oscillatory, especially in the semi-classicallimit ε → E from the right boundary into an electronic device(e.g., diode), modeled on the interval [0 , − ε ψ ′′ E ( x ) + V ( x ) ϕ ( x ) = Eψ E ( x ) , x ∈ (0 , ,ψ ′ E (0) + i k (0) ψ E (0) = 0 ,ψ ′ E (1) − i k (1) ψ E (1) = − k (1) , (1.2)where k := ε − p E − V ( x ) is the wave vector and V denotes the electrostaticpotential. Note that our assumption a ( x ) ≥ τ > E > V ( x ), so thesolution ψ E becomes oscillatory. Then, one is often interested in macroscopicquantities like the electron density n and the current density j , which aregiven by n ( x ) = Z ∞ f ( k ) | ψ E ( k ) ( x ) | d k , j ( x ) = ε Z ∞ f ( k ) ℑ (cid:16) ψ E ( k ) ( x ) ψ ′ E ( k ) ( x ) (cid:17) d k , (1.3)where f represents the injection statistics of the electrons. In order to com-pute these quantities the Schr¨odinger equation (1.2) has to be solved many2imes, as a fine grid in E is needed. Hence, efficient methods for the solutionof (1.2) are of great interest in such applications.In [1, 4, 5], efficient and accurate WKB-based (named after the physicistsWentzel, Kramers, Brillouin; cf. [6]) numerical schemes have been developedfor (1.1) in the oscillatory regime. They allow the computation of a solutionusing a coarse spatial grid with step size h > λ . In fact, the grid limitationcan there be reduced to at least h = O ( √ ε ). Now, this work includes anadaptive step size control into the algorithm from [1] as well as a couplingmechanism, which allows the algorithm to switch to a standard ODE method(e.g., Runge-Kutta) during the computation in order to avoid technical orefficiency problems in regions where the coefficient function a ( x ) is very smallor indeed equal to zero. We recall that the WKB-approximation is not validclose to turning points, i.e. where a ( x ) = 0. A coupling mechanism wasalso used in [2], where the authors presented another WKB-based numericalscheme for the initial value problem corresponding to (1.1). Therefore, onegoal of this work is to compare numerical results from our method with theresults given by the method from [2], by considering two examples where theanalytical solution is known.This paper is organized as follows: In Section 2 we give a short review ofthe second order (w.r.t. ε ) WKB-marching method from [1]. Section 3 thendescribes the adaptive step size control algorithm as well as the couplingmechanism. In Section 4 we recap the Runge-Kutta-WKB method from [2]and point out the difference between their step size control and the one usedin this paper. In Section 5 we include numerical investigations on the errorestimators of our algorithm as well as a comparison of the numerical resultsof our method and the method from [2]. We then conclude in Section 6.
2. The WKB-marching method
We aim at solving the Schr¨odinger equation (1.1), augmented with theinitial conditions ϕ ( x ) = ϕ , ϕ ′ ( x ) = ϕ ′ (2.1)with some x ∈ R . First we shall review the basics of the second order(w.r.t. ε ) WKB-marching method from [1] with focus on the algorithm. Themotivation for this method was the construction of a numerical scheme thatis uniformly correct in ε and sometimes even asymptotically correct, i.e.3he numerical error goes to zero with ε → h remainsconstant. For further details see [1]. The method consists of two parts:1. Analytic pre-processing of (1.1) by transforming the equation into asmoother (i.e. less oscillatory) problem that can be solved accuratelyand efficiently on a coarse grid.2. Obtaining a numerical scheme by discretization of the smoother prob-lem. Analytic pre-processing.
The well-known WKB-approximation (cf. [6])for the oscillatory regime where a ( x ) ≥ τ for some τ >
0, is based on insertingthe ansatz ϕ ( x ) ∼ exp ε ∞ X p =0 ε p φ p ( x ) ! (2.2)into equation (1.1). After a comparison of ε -powers one obtains the firstthree functions φ p ( x ) as φ ( x ) = ± i Z xx p a ( y ) d y , (2.3) φ ( x ) = ln( a ( x ) − ) , (2.4) φ ( x ) = ∓ i Z xx b ( y ) d y , b ( x ) := − a ( x ) (cid:16) a ( x ) − (cid:17) ′′ . (2.5)Truncation of the sum (2.2) after p = 2 leads to the second order (w.r.t. ε )asymptotic WKB-approximation of ϕ ( x ): ϕ ( x ) = c exp (cid:0) i ε φ ε ( x ) (cid:1) a ( x ) + c exp (cid:0) − i ε φ ε ( x ) (cid:1) a ( x ) , (2.6)with the phase φ ε ( x ) := Z xx (cid:16)p a ( y ) − ε b ( y ) (cid:17) d y . (2.7)In the WKB-marching method of [1] this second order WKB-approximationis used to transform (1.1) into a smoother problem. To clarify our termi-nology, we point out that this method (as well as the Runge-Kutta-WKB4ethod of Section 4) has both a WKB-order (w.r.t. ε ; referring to the usedcut-off in the asymptotic expansion (2.2)) and a numerical order (w.r.t. thestep size h ; referring to the convergence order). Firstly, using the notation U ( x ) = (cid:18) u ( x ) u ( x ) (cid:19) := a ( x ) ϕ ( x ) ε ( a ( x ) / ϕ ( x ) ) ′ √ a ( x ) , (2.8)the second order differential equation (1.1) can be reformulated as a systemof first order differential equations: ( U ′ ( x ) = (cid:2) ε A ( x ) + ε A ( x ) (cid:3) U ( x ) , x > x ,U ( x ) = U I . (2.9)Here, the two matrices A and A are given by A ( x ) := p a ( x ) (cid:18) − (cid:19) ; A ( x ) := (cid:18) b ( x ) 0 (cid:19) . Then, the first order system (2.9) for U ( x ) is transformed by the change ofvariables Z ( x ) = (cid:18) z ( x ) z ( x ) (cid:19) := exp (cid:18) − i ε Φ ε ( x ) (cid:19) P U ( x ) , (2.10)with the two matrices P := 1 √ (cid:18) i 11 i (cid:19) ; Φ ε ( x ) := (cid:18) φ ε ( x ) 00 − φ ε ( x ) (cid:19) , where φ ε is the phase function defined in (2.7). This leads to the system ( Z ′ ( x ) = ε N ε ( x ) Z ( x ) , x > x ,Z ( x ) = Z I = P U I , (2.11)where N ε ( x ) is a matrix with only off-diagonal non-zero entries: N ε , ( x ) = b ( x ) e − ε φ ε ( x ) , N ε , ( x ) = b ( x ) e ε φ ε ( x ) . { x n , n ∈ N } . Thenthe original solution can be recovered by the inverse transformation U ( x ) = P − exp (cid:18) i ε Φ ε ( x ) (cid:19) Z ( x ) . It should be noted that, throughout the whole transformation from U ( x )to Z ( x ), the phase integral φ ε ( x ) is assumed to be known exactly. For ageneralization using a spectral method to numerically compute the phase see[7]. Numerical scheme.
The derivation of the second order (in h ) scheme for(2.11) is obtained via the second order Picard approximation e Z n +1 := e Z n + ε Z x n +1 x n N ε ( x ) d x e Z n + ε Z x n +1 x n N ε ( x ) Z xx n N ε ( y ) d y d x e Z n . (2.12)Since the entries of N ε ( x ) are highly oscillatory, (2.12) involves (iterated)oscillatory integrals. With φ ε assumed to be known exactly, they are thenapproximated using similar techniques as the asymptotic method in [8]. Thefirst order (in h ) scheme for (2.11) is derived by only taking into accountthe first two terms from (2.12). For both of these schemes we introduce thefollowing notations: b ( y ) := b ( y )2 (cid:16)p a ( y ) − ε b ( y ) (cid:17) ; b k +1 ( y ) := 12 ( φ ε ( y )) ′ d b k d y ( y ) , k = 0 , , h ( y ) := e i y − h ( y ) := e i y − − i y . Further, let { x , x , ..., x N } be a grid we want to compute the solution on,and h := max ≤ n ≤ N | x n − x n − | be the step size. Then both schemes read asfollows:First order scheme: Let Z := Z I be the initial condition and let n = 0 , ..., N −
1. Then the algorithm updates as Z n +1 = (cid:0) I + A n (cid:1) Z n , (2.13)6ith the matrix A n := ε b ( x n +1 ) − ε φ ε ( x n ) h (cid:0) − ε s n (cid:1) e ε φ ε ( x n ) h (cid:0) ε s n (cid:1) ! − i ε b ( x n ) e − ε φ ε ( x n ) − b ( x n +1 ) e − ε φ ε ( x n +1 ) b ( x n +1 ) e ε φ ε ( x n +1 ) − b ( x n ) e ε φ ε ( x n ) ! , and the phase increments s n := φ ε ( x n +1 ) − φ ε ( x n ) . Second order scheme: Let Z := Z I be the initial condition and let n =0 , ..., N −
1. Then the algorithm updates as Z n +1 = (cid:0) I + A mod,n + A n (cid:1) Z n , (2.14)with the matrices A mod,n := − i ε b ( x n ) e − ε φ ε ( x n ) − b ( x n +1 ) e − ε φ ε ( x n +1 ) b ( x n +1 ) e ε φ ε ( x n +1 ) − b ( x n ) e ε φ ε ( x n ) ! + ε b ( x n +1 ) e − ε φ ε ( x n +1 ) − b ( x n ) e − ε φ ε ( x n ) b ( x n +1 ) e ε φ ε ( x n +1 ) − b ( x n ) e ε φ ε ( x n ) ! + i ε b ( x n +1 ) − e − ε φ ε ( x n ) h (cid:0) − ε s n (cid:1) e ε φ ε ( x n ) h (cid:0) ε s n (cid:1) ! − ε b ( x n +1 ) − ε φ ε ( x n ) h (cid:0) − ε s n (cid:1) e ε φ ε ( x n ) h (cid:0) ε s n (cid:1) ! , and A n := − i ε ( x n +1 − x n ) b ( x n +1 ) b ( x n +1 ) + b ( x n ) b ( x n )2 (cid:18) − (cid:19) − ε b ( x n ) b ( x n +1 ) (cid:18) h (cid:0) − ε s n (cid:1) h (cid:0) ε s n (cid:1)(cid:19) + ε b ( x n +1 ) [ b ( x n ) − b ( x n +1 )] (cid:18) h (cid:0) − ε s n (cid:1) − h (cid:0) ε s n (cid:1)(cid:19) . . Error estimation and step size control The WKB scheme is efficient in the highly oscillatory regime, but notapplicable close to turning points, i.e. zeros of a ( x ), see [9]. This is evi-dent already from the transformation (2.8), which does not make sense when a ( x ) ≤
0. For mixed problems, e.g., the Airy equation on R +0 (see Section5.1), which has a turning point at x = 0, it is therefore convenient to coupletwo different methods: a method for highly oscillatory ODEs (e.g., WKB-based) away from the turning point, and a standard ODE method (e.g.,Runge-Kutta) close to the turning point, where the solution is smooth any-how. Here, we choose the well-known Runge-Kutta-Fehlberg 4(5) (RKF45)scheme (cf. [10]) as the standard ODE method. The exact coupling mech-anism as well as the introduction of an adaptive step size control to thealgorithm will be described in the two following subsections. In order to compute the solution efficiently, an adaptive step size con-troller, based on an estimator for the local truncation error, will be addedto the numerical methods. This control allows the step size to increase ordecrease while aiming to keep the error estimator as close as possible withina given error tolerance. To illustrate the functionality of this step size con-troller, we shall consider a numerical scheme of order k . We are then able toapply this step size control individually to the different methods mentionedabove.Let Y ( k ) n = ( ϕ ( k ) n , ϕ ′ ( k ) n ) and Y ( k − n = ( ϕ ( k − n , ϕ ′ ( k − n ) be two numerical so-lutions of order k and k − Y ( x n ) =( ϕ ( x n ) , ϕ ′ ( x n )) of the initial value problem (1.1), (2.1). E.g., one could choosethe WKB schemes (2.13) and (2.14) of h -order 1 and 2, respectively. Nextwe want to decide whether to accept the numerical solution Y ( k ) n or rather toretry the computation with a modified step size. To this end we define theestimator for the local (relative) truncation error asest n := k Y ( k ) n − Y ( k − n k ∞ k Y ( k ) n k ∞ . (3.1)Let h n,trial := x n − x n − be the (trial) step size which was used to computethe solutions at the current step n . We then use the common approach ofvarying the step size via the multiplicative control h new := θ n · h n,trial . θ n to be based on the so-called elementary con-troller (e.g., see [11]). Additionally, we introduce limitations in such way thatthe step size controller responds smoother to abrupt changes in the solutionbehaviour, that is, the ratio between two consecutive step sizes should notbe exorbitantly large or small. That said, we choose the factor similar to [12,p. 310] as θ n := max . , min , (cid:18) RTolest n (cid:19) k +1 !! , (3.2)where RTol is a given relative error tolerance and the values 0 . θ n ≥ θ n ≥
1, we accept the n -th step with the step size defined as h n := h n,trial and define the trial step size for the next step as h n +1 ,trial := h new . However, if θ n <
1, the n -th step gets rejected and a reattempt is done withthe smaller (trial) step size h n,trial by updating its value as h n,trial → h new . As already mentioned above, the algorithm shall automatically switchbetween two numerical methods, the WKB method in the oscillatory regimeand another method, which is valid close to turning points. To realize thisdynamical switching mechanism we follow a similar strategy as [3]. To il-lustrate this procedure we now consider two numerical schemes of order k (1) and k (2) , where the superscripts (1) and (2) correspond to the two schemes.In each step, the adaptive step size algorithm from the previous section isapplied to both schemes individually up to the definition of the quantities θ (1) n and θ (2) n , i.e., we just evaluate (3.1) and (3.2). At this point the couplingmechanism intervenes and it selects the numerical method that yields thelarger value of θ ( i ) n , for i = 1 ,
2, hence yielding the larger proposed new stepsize. We thus favor the method with the smaller error estimator, discountedby its respective order k ( i ) . Hence we defineΘ n := max (cid:0) θ (1) n , θ (2) n (cid:1) , n -th step. Throughthis procedure the algorithm does not only use the error estimator to findthe next step size, but also to decide between the two methods.The adaptive algorithm then continues with checking the acceptance con-dition Θ n ≥
1. If the step gets accepted, the algorithm sets h n := h n,trial ,h n +1 ,trial := Θ n · h n,trial . For Θ n <
1, the step is rejected, and a reattempt is done with the smaller(trial) step size h n,trial by updating its value as h n,trial → Θ n · h n,trial .
4. The Runge-Kutta-WKB method
In this section we give a short review of the Runge-Kutta-WKB (RK-WKB) method presented in [2] for the initial value problem (1.1), (2.1). Themethod is based on a dynamic switching mechanism between a standardRunge-Kutta scheme and a new stepping procedure that uses the WKB-ansatz (2.2) as an approximation of the true solution. This stepping proce-dure reads as follows: ϕ n +1 := γ + f + ( x n + h n ) + γ − f − ( x n + h n ) , (4.1) ϕ ′ n +1 := δ + f ′ + ( x n + h n ) + δ − f ′− ( x n + h n ) , (4.2) x n +1 := x n + h n , (4.3)where γ ± := ϕ ′ n f ∓ ( x n ) − ϕ n f ′∓ ( x n ) f ′± ( x n ) f ∓ ( x n ) − f ′∓ ( x n ) f ± ( x n ) , (4.4) δ ± := ϕ ′′ n f ′∓ ( x n ) − ϕ ′ n f ′′∓ ( x n ) f ′′± ( x n ) f ′∓ ( x n ) − f ′′∓ ( x n ) f ′± ( x n ) , (4.5) ϕ ′′ n = − ε a ( x n ) ϕ n . (4.6)Here, f ± are chosen as WKB-approximations (2.2) of some finite order. Forinstance, one gets the second (WKB-)order method by setting f ± ( x ) := exp (cid:0) ± i ε φ ε ( x ) (cid:1) a ( x ) . f ± withhigher (WKB-)orders. However, the stepping procedure (4.1)-(4.6) is alwaysa first order (in h ) numerical method. This holds because the coefficients γ ± and δ ± are chosen in such way that one finds ϕ n +1 = ϕ n + ϕ ′ n h + O ( h ) ,ϕ ′ n +1 = ϕ ′ n + ϕ ′′ n h + O ( h ) , from (4.1)-(4.3), as stated in [2].It is also worth noting that in [2] the authors use a slightly different dy-namic switching mechanism and a different step size control compared to thealgorithm presented in Section 3. Firstly, their error estimator (3.1) withinthe WKB-stepping procedure uses the difference between two numerical so-lutions of different WKB-orders instead of different h -orders, simply sincethey do not have schemes of two different h -orders at their disposal. Thealgorithm then decides between a WKB step and a RK step by choosingthe method with smaller error estimator. In their more recent paper [3] theauthors use a similar step size control and switching mechanism as presentedin Section 3, but they do not limit the ratio of two consecutive step sizes asdone in (3.2).The goal of the following section is to compare numerical results fromthe WKB-marching method to results one gets using the RKWKB method.To both methods we will apply (exactly) the step size control and switchingmechanism from Section 3, for the sake of comparability. Since the WKB-stepping procedure of the RKWKB method is always of first order w.r.t. h ,the definition of the error estimator (3.1) does not make sense. Thereforewe shall use two different WKB orders (instead of h -orders) to be able tocompute the error estimator in this case, as also done in [2]. Further, sinceour algorithm consists of a different step size update formula and switchingcriterium, this may cause slightly different numerical results for the RKWKBmethod compared to [3, 2].
5. Numerical results: WKB-marching method vs. the Runge-Kutta-WKB method
In this section we will compare numerical results of the WKB-marchingmethod and of the RKWKB method by applying both algorithms to twoexamples, where exact analytical solutions are available. The first example11orresponds to a linear coefficient function a ( x ) and is taken from [2, 3],whereas the second example involves a quadratic function a ( x ) and appearsin [9]. In both examples the phase integral (2.7) in the WKB basis func-tions (2.3)-(2.5) can easily be computed exactly, hence we do not use anynumerical integration routine here. By contrast, in [3] they evaluate theWKB basis functions with a numerical integration routine and therefore getanother source of error in their method. Moreover, we will always use thesecond order WKB-marching method, since no scheme of higher WKB-orderhas been developed yet, see [1]. All simulations in this paper are done with Matlab version 9.7.0.1319299 (R2019b).
The first example investigated in [2] is the Airy equation ϕ ′′ ( x ) + xϕ ( x ) = 0 , x > ,ϕ (0) = − / +i 3 − / Γ( ) ,ϕ ′ (0) = − / − i 3 / Γ( ) , (5.1)which results if one chooses the coefficient function a ( x ) = x as well as ε = 1.Here, the exact solution is given by ϕ exact ( x ) = Ai( − x ) + i Bi( − x ) , where Ai and Bi denote the Airy functions of first and second kind, re-spectively. This example demonstrates very well the advantages of a WKBmethod, since the solution becomes more and more oscillatory for large valuesof x . While standard adaptive ODE methods, e.g., Runge-Kutta methods,would need to decrease the step size more and more, a WKB-based methoddoes not have to resolve the individual oscillations. Actually, it even allowsto increase the step size for large x and is therefore highly efficient. This canbe seen in Figures 4-6.As a first step we shall now test the reliability of our choice of errorestimator (3.1) – but only for the WKB steps of the two methods from Section2 and Section 4. Since we know the exact solution to this problem, we cancompute the local truncation error in each step and are able to compare itto the error estimator. The results can be seen in Figures 1 and 2.12 -10 -5 Figure 1: The error estimator (3.1) in comparison to the actual (relative) local truncation error computedusing the WKB-marching method of first order in h , (2.13), and second order in h , (2.14). -10 -5 Figure 2: The error estimator (3.1) in comparison to the actual (relative) local truncation error computedusing the RKWKB method. Since the method is of first order (in h ), the estimator as well as the localerrors were computed by using different WKB-orders instead (here order 2 and 3). According to the results of Figures 1 and 2, the error estimator is inexcellent agreement with the local truncation error (for this example). Hence,it seems to be an adequate choice. We also find a very good agreement byplotting the respective relative errors for one single step as functions of thestep size h for a fixed starting point x , see Figure 3. For the RKWKBmethod the error estimator and the true local error are indistinguishable onthe shown scale. 13 -8 -7 -6 -5 -4 Figure 3: The error estimators (3.1) in comparison to the actual (relative) local truncation errors for bothmethods in one single step as functions of the step size h for a fixed starting point x = 10. For thecomputation of the estimator with the RKWKB method a third order WKB-ansatz was used. Remark 5.1.
The computation of the local truncation error involves theevaluation of the Airy function, which seems to have an accuracy problem forlarge values of x when using the standard routine airy() in Matlab . Hencewe used a modified implementation of the Airy functions, which consists ofthe function airy() from
Matlab for small values of x and asymptoticexpressions for the Airy function for large values of x . More precisely, theevaluations were performed as given in Table 1. The asymptotic expressionsare based on well-known expansions, which can be found in the Appendix.The order for truncating the expansions was set to K = 3 . airy() asymptoticsAi x ∈ [0 , x ∈ (500 , ∞ )Ai ′ x ∈ [0 , x ∈ (400 , ∞ )Bi x ∈ [0 , x ∈ (500 , ∞ )Bi ′ x ∈ [0 , x ∈ (400 , ∞ ) Table 1: Intervals for evaluating the Airy function of first and second kind as well as their derivatives.For small values of x the original function from Matlab was used, and for large x the evaluation wasperformed using the asymptotic expressions (A.1)-(A.4) with truncation after K = 3. Now, we will give numerical results for solving the initial value problem(5.1) for the Airy equation. We recall that the algorithm automaticallychooses between RKF45 and the respective WKB steps. To illustrate this14ifference in the Figures 4-10, we will mark RKF45 steps with red dots, WKBsteps using the WKB-marching method with blue squares, and WKB stepsusing the RKWKB method with green triangles. -6 -5 -4 -3 Figure 4: Real part of the numerical solution obtained by using the WKB-marching method includingthe Runge-Kutta switching mechanism compared to the (exact) reference solution (solid line). The initialstep size was set to h ,trial = 0 . − . At x ≈ . -6 -5 -4 -3 Figure 5: Real part of the numerical solution obtained by using the RKWKB method compared to the(exact) reference solution (solid line). The initial step size was set to h ,trial = 0 . − . A third order WKB-ansatz was used. At x ≈ . According to Figures 4 and 5 the WKB-marching method seems to per-form slightly better than the RKWKB method in this example, in matters15f global error. But at the same time the WKB-marching method only needshalf as many steps as the RKWKB method. Also, within the algorithm us-ing the WKB-marching method, the switch from RKF45 steps to WKB stepshappens earlier as can be seen in Figure 4. In Figure 6, the relative errors ofboth algorithms are compared on [0 , ]. We find that the algorithm usingthe WKB-marching method made fewer steps (57 vs. 86) while producinga lower global error at the same time. The almost identical grid spacing ofboth methods from x ≈
300 onwards is due to limiting the quotient of twoconsecutive step sizes, imposed in both methods. -6 -5 -4 -3 Figure 6: Relative error comparison of the WKB-marching method and RKWKB on [0 , ] using aninitial step size h ,trial = 0 .
5. The relative error tolerance was set to RTol = 10 − . A third order WKB-ansatz was used for the RKWKB method. Overall the algorithm using the WKB-marching method made57 steps, whereas RKWKB made 86 steps. Furthermore, Figures 7 and 8 show the ratio between the WKB- and RK-error estimators as well as the ratio of the proposed step sizes (by the WKBand RKF schemes) for each step of the above simulations.16 -5 Figure 7: Numerical computation using the adaptive WKB-marching method: The magenta line indicatesthe ratio of the two error estimators (due to the WKB and RKF scheme). The green line gives theanalogous ratio of the step sizes (locally) proposed by these two methods. -5 Figure 8: Numerical computation using the RKWKB method: The magenta line indicates the ratio of thetwo error estimators (due to the WKB and RKF scheme). The green line gives the analogous ratio of thestep sizes (locally) proposed by these two methods.
For both WKB methods (from Section 2 and Section 4), these plotsdemonstrate well the superiority of the RK scheme in the less oscillatoryregime, i.e. for x small, as the ratio of the error estimators is very largethere. The switching points to the WKB schemes are well defined (again forboth WKB methods) – due to the monotonous behaviour of both the ratioof error estimators and the ratio of proposed step sizes. Note also that the17tep size ratios are bounded from above and below because of the limitationsin (3.2). Remark 5.2.
The evaluation of the Airy function is numerically ill posed forlarge values of x , even with the modified implementation using the asymptoticexpansions. This is because the function φ ( x ) gets very large and thereforethe evaluation of (2.6) leads to high inaccuracies within double precision. Forinstance, computing exp(i x ) for x = 10 within double precision resultsin a relative error of approximately . Therefore, we recommend to solvethe Airy equation (5.1) within double precision only until x = 10 . The sameproblem arises for trigonometric functions such as sin or cos .5.2. Second Example: Parabolic cylinder function The second example can be found in [9] and includes a quadratic coeffi-cient function a ( x ) as well as the parameter 0 < ε ≪ ε ϕ ′′ ( x ) + (cid:0) − x + x (cid:1) ϕ ( x ) = 0 , x ∈ [0 , ,ϕ (0) = κU ( ν, z (0)) ,ϕ ′ (0) = − κ − ε − U ′ ( ν, z (0)) , (5.2)with ν := − √ ε , z ( x ) := 2 √ ε (1 − x ) , and κ := 2 U ( ν, − i √ ε U ′ ( ν, . The exact solution reads ϕ exact ( x ) = κU ( ν, z ( x )) , where U ( ν, z ) denotes the parabolic cylinder function (PCF) (cf. [13, § x = 0 and x = 2. Therefore, weexpect the two WKB methods (from Section 2 and Section 4) to use RKF45steps near the turning points and WKB steps between them. Numericalresults for the specific choice ε = 2 − are presented in Figures 9-11: As beforewe compare the two WKB methods (WKB-marching method, RKWKB).18oreover, we also include the result when only using (standard) Runge-Kutta steps. -6 -4 -2 Figure 9: ε = 2 − : Real part of the numerical solution obtained by using the WKB-marching methodcompared to the reference solution (solid line). The initial step size was set to h ,trial = 0 .
05 and therelative error tolerance was set to RTol = 10 − . The algorithm made 73 steps. -6 -4 -2 Figure 10: ε = 2 − : Real part of the numerical solution obtained by using the RKWKB method comparedto the reference solution (solid line). The initial step size was set to h ,trial = 0 .
05 and the relative errortolerance was set to RTol = 10 − . A third order WKB-ansatz was used. The algorithm made 129 steps. -6 -4 -2 Figure 11: ε = 2 − : Real part of the numerical solution obtained by using just the RKF45 methodcompared to the reference solution (solid line). The initial step size was set to h ,trial = 0 .
05 and therelative error tolerance was set to RTol = 10 − . The algorithm made 228 steps. These simulations illustrate the superiority of the WKB-marching methodin several respects: For this example it exhibits the lowest global relative errorat the end of the interval (among the three methods presented), but it alsoneeds significantly fewer steps. Similar to the previous example there werealmost only half of the steps needed in comparison to the RKWKB method,while now producing even a significantly smaller global error (by a factor ∼ ). Also note the large oscillations of the global error within the moreoscillatory part of the solution, when using either the pure Runge-Kutta-Fehlberg method or the RKWKB method. Remark 5.3.
The computation of the reference solution involves the evalu-ation of the PCF, which is not readily available in
Matlab . But the PCFcan be related to the Kummer confluent hypergeometric function F (see [13, § Matlab as hypergeom() .
6. Conclusion
In this paper we have introduced an extension to the WKB-marchingmethod from [1] by including into the algorithm an adaptive step size con-troller as well as a coupling mechanism to a standard ODE-solver. Thiscoupling mechanism ensures well defined switching points between WKBsteps and Runge-Kutta-Fehlberg 4(5) steps for oscillatory and, respectively,smoother regions of the solution to the ODE. In numerical simulations based20n two examples this method yielded smaller and smoother global errorsin comparison to the Runge-Kutta-WKB method from [2, 3], an alterna-tive WKB-based scheme. Especially for the Airy equation on the spatialinterval [0 , ] the efficiency of the method is demonstrated very well, asthe scheme skips millions of oscillations within one step, while staying accu-rate at the same time. A corresponding Matlab program is available in aGitHub repository , which also offers the possibility to compute the phase(2.7) numerically, as described in [7]. Acknowledgements
The authors A. Arnold and K. D¨opfner were partially supported by thebinational FWF-project I3538-N32. Moreover, the authors J. K¨orner andK. D¨opfner have been (partially) supported from the Austrian Science Fund(FWF) through grant number W1245.
Appendix A. Asymptotic formulas for Airy functions
For real-valued x and x → ∞ , asymptotic expansions for the Airy func-tions and their first derivatives are given in [13, § − x ) ∼ √ πx cos (cid:16) ζ − π (cid:17) ∞ X k =0 ( − k u k ζ k + sin (cid:16) ζ − π (cid:17) ∞ X k =0 ( − k u k +1 ζ k +1 ! , (A.1)Ai ′ ( − x ) ∼ x √ π sin (cid:16) ζ − π (cid:17) ∞ X k =0 ( − k v k ζ k − cos (cid:16) ζ − π (cid:17) ∞ X k =0 ( − k v k +1 ζ k +1 ! , (A.2)Bi( − x ) ∼ √ πx − sin (cid:16) ζ − π (cid:17) ∞ X k =0 ( − k u k ζ k + cos (cid:16) ζ − π (cid:17) ∞ X k =0 ( − k u k +1 ζ k +1 ! , (A.3)Bi ′ ( − x ) ∼ x √ π cos (cid:16) ζ − π (cid:17) ∞ X k =0 ( − k v k ζ k + sin (cid:16) ζ − π (cid:17) ∞ X k =0 ( − k v k +1 ζ k +1 ! . (A.4) https://github.com/JannisKoerner/adaptive-WKB-marching-method ζ := x and the coefficients u k and v k are given by (see [13, § u = v = 1 ,u k = (2 k + 1) · (2 k + 3) · (2 k + 5) · ... · (6 k − k k != (6 k − · (6 k − · (6 k − k − k u k − , k = 1 , , ... ,v k = 6 k + 11 − k u k , k = 1 , , ... . References [1] A. Arnold, N. Ben Abdallah, C. Negulescu,
WKB-based schemes for theoscillatory 1D Schr¨odinger equation in the semi-classical limit , SIAM J.Numer. Anal. 49 (2011), no. 4, pp. 1436-1460.[2] W. J. Handley, A. N. Lasenby, and M. P. Hobson,
The Runge-Kutta-Wentzel-Kramers-Brillouin method , arXiv:1612.02288[3] F. J. Agocs, W. J. Handley, A. N. Lasenby, and M. P. Hobson,
An effi-cient method for solving highly oscillatory ordinary differential equationswith applications to physical systems , Phys. Rev. Res. 2 , 013030 (2020),arXiv:1906.01421.[4] T. Jahnke, C. Lubich,
Numerical integrators for quantum dynamics closeto the adiabatic limit , Numerische Mathematik, 94 (2003), pp. 289-314.[5] K. Lorenz, T. Jahnke, C. Lubich,
Adiabatic integrators for highly oscil-latory second-order linear differential equations with time-varying eigen-decomposition , BIT, 45 (2005), no. 1, pp. 91-115.[6] L. D. Landau, E.M. Lifschitz,
Quantenmechanik , Akademie-Verlag,Berlin, 1985[7] A. Arnold, C. Klein, B. Ujvari,
WKB-method for the 1D Schr¨odingerequation in the semi-classical limit: enhanced phase treatment , submit-ted, (2019). arxiv.org/abs/1808.01887[8] A. Iserles, S.P. Norsett, S. Olver,
Highly oscillatory quadrature: Thestory so far , in A. Bermudez de Castro, ed.,
Proceeding of ENuMath,Santiago de Compostella (2006) , Springer Verlag, (2006), pp. 97-118.229] A. Arnold, K. D¨opfner,
Stationary Schr¨odinger equation in the semi-classical limit: WKB-based scheme coupled to a turning point , Calcolo57, No. 1 (2020) 44 pp.[10] E. Hairer, S. Nørsett, G. Wanner,
Solving Ordinary Differential Equa-tions I , Springer, Berlin, 2000[11] G. S¨oderlind,
Automatic Control and Adaptive Time-Stepping , Numeri-cal Algorithms. 31., pp. 281-310, 2002[12] J. C. Butcher,
Numerical Methods for Ordinary Differential Equations ,Wiley-Interscience, Second Edition, New York, 2008[13] F. W. Olver, D. W. Lozier, R. F. Boisvert, C. W. Clark,