A frequency-dependent p -adaptive technique for spectral methods
aa r X i v : . [ m a t h . NA ] S e p A frequency-dependent p -adaptive technique for spectral methods Mingtao Xia a , Sihong Shao b, ∗ , Tom Chou a, ∗ a Department of Mathematics, UCLA, Los Angeles, CA 90095-1555, USA b LMAM and School of Mathematical Sciences, Peking University, Beijing 100871, CHINA
Abstract
When using spectral methods, a question arises as how to determine the expansion order,especially for time-dependent problems in which emerging oscillations may require adjust-ing the expansion order. In this paper, we propose a frequency-dependent p -adaptive tech-nique that adaptively adjusts the expansion order based on a frequency indicator. Usingthis p -adaptive technique, combined with recently proposed scaling and moving techniques,we are able to devise an adaptive spectral method in unbounded domains that can captureand handle diffusion, advection, and oscillations. As an application, we use this adaptivespectral method to numerically solve the Schr¨odinger equation in the whole domain andsuccessfully capture the solution’s oscillatory behavior at infinity. Keywords: unbounded domain, spectral method, adaptive method, Schr¨odingerequation, Jacobi polynomial, Hermite function, Laguerre function
1. Introduction
Unbounded domain problems arise in many scientific applications [1, 2] and adaptivenumerical methods are needed on many occasions, for instance, in solving the Schr¨odingerequation in unbounded domains when the solution’s behavior varies over time and we wishto capture the solution’s behavior in the whole domain. As an important class of numer-ical algorithms, adaptive methods have witnessed numerous advances in their efficiencyand accuracy [3, 4, 5, 6]. However, despite considerable progress that has been made forspectral methods in unbounded domains [7], there are few adaptive methods that apply inunbounded domains.In [8], adaptive scaling and moving techniques were proposed for spectral methodsin unbounded domains and it was noted that adjusting the expansion order is necessarywhen the function displays oscillatory behavior that varies over time. In this paper, we first ∗ Corresponding author.
Email addresses: [email protected] (Mingtao Xia), [email protected] (Sihong Shao), [email protected] (Tom Chou)
Preprint submitted to Elsevier October 6, 2020 evelop a frequency-dependent technique for spectral methods which adjusts the expansionorder N ( N + 1 basis functions are used to approximate the solution). This techniquetakes advantage of the frequency indicator defined in [8] and corresponds to p -adaptivity[9, 10, 11]. By adjusting the expansion order efficiently, our p -adaptive technique can beused to accurately solve problems with varying oscillatory.By combining this p -adaptive technique with scaling and moving methods, we developan adaptive spectral method that can capture diffusion, advection, and oscillations inunbounded domains. Since scaling and adjusting the expansion order both depend onthe frequency indicator, we also investigate the interdependence of these two techniques.We demonstrate that appropriately adjusting the expansion order can facilitate scalingto more efficiently distribute allocation points. In turn, proper scaling can help avoidunnecessary increases in the expansion order when it does not increase accuracy, therebyavoiding unnecessary computational burden.The significance of this adaptive spectral method is that it can capture the solution’sbehavior in the whole domain. We demonstrate the utility of our method by solvingSchr¨odinger’s equation in R . Here, the unboundedness and the oscillatory nature of thesolution pose two major numerical challenges [12]. Specifically, in the semiclassical regime,when the wavelength of the solution is small, the function becomes extremely oscillatory.Moreover, in certain situations, one has to work with a very large computational domainthat is difficult to automatically determine.Previous numerical methods which solve Schr¨odinger’s equation in unbounded domainsusually truncate the domain into a finite subdomain and impose artificial boundary con-ditions, which may be nonlocal and complicated [13, 12, 14, 15]. Our adaptive spectralmethod tackles the oscillatory problem directly in the original unbounded domain withoutthe need to truncate it or to devise an artificial boundary condition.This paper is organized as follows. In the next section, we first present a p -adaptivetechnique for spectral methods and use examples to illustrate its efficiency. In Section 3,we incorporate and study this technique within existing scaling and moving techniques anddevise an adaptive spectral method in unbounded domains. Application of our adaptivespectral methods to numerically solving Schr¨odinger’s equation is given in Section 4. Wesummarize our results in Section 5 and propose directions for future work.
2. Frequency-dependent p -adaptivity We present a frequency-dependent p -adaptive spectral method based on informationextracted from only the numerical solution of time-dependent problems. In [8], we showedthat a frequency indicator defined for spectral methods is particularly useful in measuringthe contribution of high frequency modes in the numerical solution. Because high frequencymodes decay more slowly, this indicator could be used to determine scaling in spectralmethods applied to unbounded domains. In this work, we will show that the frequency2ndicator can also be used to determine whether more or fewer basis functions are neededto refine or coarsen the numerical solution.Given a set of orthogonal basis functions { B i ( x ) } ∞ i =0 under a specific weight function ω ( x ) > I N u ( x ) = U N ( x ) = N X i =0 u i B i ( x ) (2.1)is defined as in [8] F ( U N ) := N P i = N − M +1 γ i u iN P i =0 γ i u i , (2.2)where γ i = R Λ B i ( x ) ω ( x )d x is the square of L ω -weighted norm of the basis function B i ( x ).This frequency indicator measures the contribution of the M highest-frequency componentsto the L ω -weighted norm of U N . Here M is often chosen to be [ N ] following the -rule[16, 17]. This indicator provides a lower bound for the error divided by the norm ofthe numerical solution k u −I N − M u k ω kI N u k ω which is illustrated in [8]. Thus, the quality of thenumerical interpolation U N can be measured by F ( U N ).For a time-dependent problem, the expansion order N may need adjusting dynamically,which can be reflected by the frequency indicator. If the frequency indicator increases, thelower bound for k u −I N − M u k ω kI N u k ω will also increase. On the other hand, as N increases, theerror k u − I N u k ω as well as F ( U N ) are expected to decrease. By sufficiently increasing theexpansion order N , the frequency indicator as well as the error can be kept small. If thefrequency indicator decreases, we can also consider decreasing N to relieve computationalcost without compromising accuracy, as was done in [9]. The pseudo-code of the proposed p -adaptive technique is given in Alg. 1.The p -adaptive spectral method in Alg. 1 for time-dependent problems consists of twoingredients: refinement (increasing N ) and coarsening (decreasing N ). It maintains ac-curacy when there are emerging oscillations by increasing the expansion order N . It alsodecreases N when the expansion order is larger than needed to avoid unnecessary com-putation. In Alg. 1, the frequency indicator subroutine is to calculate the frequencyindicator defined in Eq. (2.2) for the numerical solution U N while the evolve subroutineis to obtain the numerical solution U N ( t + ∆ t ) at the next timestep from U N ( t ).In Line 11 of Alg. 1, the refine subroutine uses U N to generate a new numerical solutionwith a larger expansion order U N +1 (refine), and in Line 20 the coarsen subroutine uses U N to generate a new numerical solution with a smaller expansion order U N − (coarsen).The refinement or coarsening is achieved by reconstructing the function values of U N +1 or U N − at the new set of collocation points { x i } :3 lgorithm 1 Pseudo-code of the p -adaptive technique which may increase (refine) or decrease(coarsen) the expansion order N . Initialize
N, N , γ ≥ , η = η >
1, ∆ t , T , α , β , U N (0), N max , N min t ← f ← frequency indicator ( U N ( t )) while t < T do U N ( t + ∆ t ) ← evolve ( U N ( t ) , ∆ t ) f ← frequency indicator ( U N ( t + ∆ t )) l ← if f > ηf then while f > ηf and l ≤ N max do l ← l + 1 U N +1 ← refine ( U N ( t + ∆ t )) N ← N + 1 f ← frequency indicator ( U N ) end while f ← f η ← γη η else if f < f /η then r ← False while f < f /η and N > N min and not r do ˜ U N − ( t + ∆ t ) ← coarsen ( U N ( t + ∆ t )) f ← frequency indicator ( ˜ U N − ( t + ∆ t )) if f < f then f ← f r ← True U N − ( t + ∆ t ) ← ˜ U N − ( t + ∆ t ) N ← N − end if end while if r then f ← f end if end if t ← t + ∆ t end while N ± ( x i , t ) = U N ( x i , t ) , i = 0 , ..., N ± , (2.3)where U N +1 uses N + 2 basis functions for refinement and U N − uses N basis functions forcoarsening.In Alg. 1, ηf is the refinement threshold such that if the current frequency indicator f > ηf , we increase the expansion order N . The while loop starting in Line 9 ensures weeither refine enough such that the frequency indicator, after increasing N , is smaller thanthe threshold ηf , or the maximal allowable expansion order increment within a single step N max is reached.After increasing N , f is renewed to be the current frequency indicator and η is multi-plied by a factor γ ≥
1, enabling us to dynamically adjust the refinement threshold for thenext refinement in order to prevent increasing N too fast without substantially increasingaccuracy. On the other hand, when an extremely large N is needed to match the increas-ingly oscillatory behavior of the numerical solution, we can set γ (cid:23) γ = 1, as wewill do in Examples 5 and 6. We have observed numerically, as expected, that the larger η , γ are, the more difficult it is to increase the expansion order.We also consider reducing N when a large expansion order is not really needed and f /η is the threshold for decreasing the expansion order. If the condition in Line 17 issatisfied and N > N min , the minimal allowable expansion order, and we have not increased N in the current step, on the contrary we consider decreasing the expansion order belowLine 17. As long as the frequency indicator of the new numerical solution with the de-creased expansion order F ( U N − ) is smaller than f , the frequency indicator recorded afterpreviously adjusting the expansion order, reducing the expansion order is accepted; elsereducing the expansion order is declined. Therefore, f after coarsening will not surpass f before coarsening. This procedure is described by the If condition in Line 22. If N isdecreased, f will also get renewed to be the latest frequency indicator.In addition, if the current frequency indicator f ∈ [ f η , ηf ], neither the refinement northe coarsening subroutine is activated.Alg. 1 can be generalized to higher dimensions in a dimension-by-dimension manner.The expansion order for each dimension can change simultaneously within each timestepby using the tensor product of one-dimensional basis functions, in much the same waymoving and scaling algorithms were generalized to higher dimensions [8]. For example, fora two-dimensional problem, given U ~N ( x, y ) := N x X i =0 N y X j =0 u i,j B i ( x ) B j ( y ) (2.4)where ~N = ( N x , N y ), the frequency indicator in the x -direction is defined as5 able 1: Typical choices of basis functions { B i } ∞ i =0 and computational domain Λ. Computational domain Bounded interval (0 , ∞ ) ( −∞ , ∞ )Basis functions Jacobi polynomials Laguerre polynomials/functions Hermite polynomials/functions F x ( U ~N ) := N x P i = N x − M x +1 N y P j =0 γ i γ j u i,jN x P i =0 N y P j =0 γ i γ j u i,j , (2.5)while the frequency indicator in y -direction is similarly defined. At each timestep, we keep N y fixed and use F x to judge whether or not to renew N x → ˜ N x ; simultaneously, we fix N x and use F y to renew N y → ˜ N y if adjusting the expansion order in y dimension is needed.Finally N x , N y are updated to ˜ N x , ˜ N y .In this work, the relative L ω -errorError = k U N − u k ω k u k ω , (2.6)is used to measure the quality of the spectral approximation U N ( x ) compared to the ref-erence solution u ( x ). Table 1 lists some typical choices of orthogonal basis functions fordifferent domains Λ that we use in this paper.We provide two examples of using this p -adaptive technique in Alg. 1 below, where thegeneralized Jacobi polynomials [18] are used. Theorem 3.41 in [18] gives an estimation forthe interpolation error of a function u in the Jacobi-weighted Sobolev space for α, β > − k ∂ lx ( I N,α,β u − u ) k ω α + l,β + l ≤ c r ( N − m + 1)! N ! N l − ( m +1) / k ∂ mx u k ω α + m,β + m , ≤ l ≤ m ≤ N + 1 , (2.7)where c is a positive constant independent of m, N and u . When m > l = 0, theleft hand side becomes the interpolation error k ( I N,α,β u − u ) k ω α,β which decreases with N .Therefore, by increasing the expansion order for the Jacobi polynomials it is generally truethat the interpolation will be more accurate. Theorem 7.16 and Theorem 7.17 in [18] givesimilar error estimates for Laguerre and Hermite interpolations, which reveals that undersome smoothness assumptions, the interpolation error decreases when the expansion order N increases.Since unbounded domain problems may involve diffusive and advective behavior, wediscuss and develop adaptive spectral methods in unbounded domains in the next section.6 -10 -5
10 2 4 610 -10 -5 Figure 1:
Numerically solving Eq. (2.8) with Chebyshev polynomials using Alg. 1. For solutionsthat become increasingly oscillatory, the p -adaptive technique can increase the expansion ordereffectively to capture the oscillations and maintain a small error by keeping the frequency indicatorlow. Using a fixed N fails to maintain the frequency indicator and results in a large error. xample 1. We numerically solve the PDE ∂ t u = (cid:18) x + 2 t + 1 (cid:19) ∂ x u, x ∈ [ − , , (2.8)with a Dirichlet boundary condition specified at x = 1 given as u (1 , t ) = cos 3( t + 1). ThisPDE admits an analytical solution u ( x, t ) = cos(( t + 1)( x + 2)) . (2.9)We solve it numerically by using Chebyshev polynomials with Chebyshev-Gauss-Robattoquadrature nodes and weights. The Chebyshev polynomials are orthogonal under theweight function ω ( x ) = (1 − x ) − , i.e. , they correspond to Jacobi polynomials with α = β = − . Since u ( x, t ) becomes increasingly oscillatory over time, an increasingexpansion order is required to capture these oscillations. We start with N = 10 at t = 0,the parameters η = 1 . , γ = 1 . , N max = 3 , N min = 0, and a timestep ∆ t = 0 . u ( x, t ) is plotted in Fig. 1(a). The increasing oscillations leadto a fast rise in the frequency indicator as the contribution from high frequency modesincreases. Keeping the same number of basis functions over time will fail as it will beeventually incapable of capturing the shorter wavelength oscillations.However, a much more accurate approximation can be obtained (see Fig. 1(b)) with our p -adaptive method which maintains the frequency indicator (see Fig. 1(c)) by increasing thenumber of basis functions (shown in Fig. 1(d)). Furthermore, the coarsening subroutine fordecreasing the expansion order described in the while loop in Line 17 will not be triggered(shown in Fig. 1(d)).When directly approximating the reference solution in Eq. (2.9), we can achieve 10 − accuracy with only 20 basis functions. However, when numerically solving Eq. (2.8), theerror will accumulate due to the increasing oscillatory behavior which will require evenmore basis functions to achieve the same accuracy as the direct approximation to Eq. (2.9).Thus, the oscillatory behavior of the solution poses additional difficulties and requires evenmore refinement when numerically solving a PDE.Next, we present an example of a two-dimensional problem in [ − , . Example 2.
We approximate the function u ( x, y, t ) = cos (cid:0) xy (5 − | t − | ) (cid:1) + y − | t − / | sin (cid:0) x (5 − | t − | ) (cid:1) , ( x, y ) ∈ [ − , (2.10)by Legendre polynomials (corresponding to Jacobi polynomials with α = β = 0) withLegendre-Gauss-Robatto quadrature nodes and weights in both dimensions. Within t ∈ [0 , ], the function becomes more oscillatory over time in both dimensions, requiring in-creasing expansion orders. For t ∈ [ , ], the error for approximation with fixed expansion8 -10 -5 -15 -10 -5 -15 -10 -5 Figure 2:
Using the p -adaptive technique to approximate the two-dimensional function in Eq. (2.10)with Legendre polynomials. Refinement is applied in each direction simultaneously to captureincreasing oscillations in both directions. Coarsening is applied when large expansion orders arenot needed. Anisotropic oscillatory behavior requires adjusting the expansion order in each directiondifferently. The frequency indicators in both dimensions are kept low, leading to a small error. x and y and the adjustment of expansion is anisotropic. We show that Alg. 1 can appropri-ately increase N x , N y when t < and reduce N x , N y when t ≥ . We take N x = N y = 36at t = 0 with a timestep ∆ t = 0 .
01, and γ x = γ y = 1 . , η x = η y = 1 . , N max ,x = N max ,y =3 , N x, min = N y, min = 0.It is clear from Fig. 2(a) that fixing the number of basis functions in each dimen-sion leads to an approximation which deteriorates while the proposed p -adaptive spectralmethod is able to keep the error small. Furthermore, when t ∈ [ ,
5] we see that with fixedexpansion order the approximation error decreases, indicating that coarsening may be per-formed to relieve computational burden while maintaining accuracy. Alg. 1 first tracks theincreasing oscillation by increasing expansion orders in both x and y dimensions. When t ≥ , Alg. 1 senses a decrease in frequency indicator and decreases both N x and N y adaptively (shown in Figs. 2(b)) without compromising accuracy, as shown in Fig. 2(a).Since sin(4 x (5 − | t − | )) is the most oscillatory term in u ( x, y, t ), the function becomesmore oscillatory in x than in y when t ∈ [0 , ]. Because the function displays different os-cillatory behavior in x - and y -directions, the expansion order should be adjusted anisotrop-ically and N x needs increasing more than N y in order to maintain F x small. Both F x and F y are maintained well for the p -adaptive approximation over time (shown in Fig. 2(c) and(d)), leading to satisfactory error control. Overall, Alg. 1 preserves accuracy for all timeswhile still avoids using excessive values of N x and N y when they are not needed.
3. Adaptive spectral methods in unbounded domains
Unbounded domain problems are often more difficult to numerically solve than boundeddomain problems. Diffusion and advection in unbounded domains necessitates knowledgeof the solution’s behavior at infinity. To distinguish and handle diffusive and advectivebehavior in unbounded domains, techniques for scaling and moving basis functions areproposed in [8]. When combining scaling, moving, refinement and coarsening, we candevise a comprehensive adaptive spectral approach for unbounded domains. A flow chartof our overall approach is given in Fig. 3. The scaling, refinement and coarsening techniquesall rely on a common frequency indicator.As is stated in [8], advection may cause a false increase in the frequency indicator.Thus, we must first compensate for advection by the moving technique before we considereither scaling or adjusting the expansion order. Next, as the cost of changing the scalingfactor is lower than increasing the expansion order, we implement scaling before adjustingthe expansion order. Only if scaling cannot maintain the frequency indicator below therefinement threshold do we consider increasing the expansion order. Coarsening is also con-sidered after scaling if the frequency indicator decreases below the threshold for coarseningwhile no refinement is performed in the current timestep.10 nitialize N , ∆ t , T , β , U ( β ) N (0), x L , x R = x ( β )[ N +23 ] f , f ← frequency indicator ( U ( β ) N,x L ( x, t )) e ← exterior error indicator ( U ( β ) N,x L ( x, , x R ) t < T End MOVE?Renew e , x L , x R SCALE?Renew β, f , x R REFINE or COARSEN? t = t + ∆ t REFINE?Renew η Renew e , f , f , x R , N YesNo NoYesYes NoNoYes NoYes
Figure 3:
Flow chart of an adaptive spectral method in unbounded domains which consists of moving,scaling, refinement and coarsening techniques.
11s we have done in [8], we also decrease the scaling factor β by multiplying it by acommon ratio q < f > νf . The scaling we perform here contains an additional step: when the currentfrequency indicator decreases and is below f , we consider increasing the scaling factor β by dividing it by the common ratio q as long as the frequency indicator decreases afterincreasing β . When f ∈ [ f , νf ], β is neither increased nor decreased. Thus, at each step,the scaling factor β may be either increased or decreased as long as the frequency indicatordecreases after adjusting the scaling factor. A decrease in the scaling factor indicates thatthe allocation points are more efficiently distributed. These changes avoid unnecessarycomputational burden that may arise if N is excessively increased. We briefly describe ourmodified scaling subroutine for one timestep in Alg. 2.For simplicity, we assume that the function is moving rightward so we need to movethe basis functions rightward. Therefore, ( x R , ∞ ) is the “exterior domain” of the spectralapproximation on which we wish to control the error as illustrated in [8]. For Laguerrepolynomials/functions the parameter x L in the algorithm in Fig. 3 denotes the startingpoint for the approximation, while for Hermite polynomials/functions x L represents thetranslation of Hermite polynomials/functions, i.e. , we use {H i ( x − x L ) } or { ˆ H i ( x − x L ) } .After the expansion order N has changed, we need to renew both the threshold for scalingand the threshold for adjusting the expansion order. We also need to renew the thresholdfor moving, denoted by e , after N has changed because different N leads to different x R . U ( β ) N,x L is the spectral approximation with the scaling factor β . Here the exterior-errorindicator for the semi-unbounded domain is defined in [8] and we can generalize it to R when using Hermite polynomials/functions E ( U ( β ) N,x L , x R ) = k ∂ x U ( β ) N,x L · I ( x R , + ∞ ) k ω β k ∂ x U ( β ) N,x L · I ( −∞ , + ∞ ) k ω β , (3.1)where ω β is the weight function and x R is taken to be x ( β )[ N +23 ] for Hermite functions/polynomialsand x ( α,β )[ N +23 ] for Laguerre functions/polynomials [8] in view of the often-used -rule. Thedifference between the choices of x R for Hermite and Laguerre basis functions arises be-cause the allocation points for Hermite functions are symmetrically distributed aroundtheir center while those for Laguerre functions are one-sided, to the right of the startingpoint x L in the axis.For the scaling subroutine we need the following parameters: the common ratio q < ν ; a predetermined lower bound forthe scaling factor β and an upper bound β . For the moving subroutine, required parametersinclude the minimal displacement for the moving technique δ , the maximal displacementwithin a single timestep d max , and the parameter of the threshold for activating the moving12 lgorithm 2 Pseudo-code of the frequency-dependent scaling technique which may increase ordecrease the scaling factor β . f ← frequency indicator ( U ( α,β ) N ( t + ∆ t )) if f > νf then β ˜ β ← qβ U ( α, ˜ β ) N ← scale ( U ( α,β ) N ( t + ∆ t ) , ˜ β ) ˜ f ← frequency indicator ( U ( α, ˜ β ) N ) while ˜ f ≤ f and ˜ β ≥ β do β ← ˜ β U ( α,β ) N ( t + ∆ t ) ← U ( α, ˜ β ) N f ← ˜ f f ← ˜ f ˜ β ← qβ U ( α, ˜ β ) N ← scale ( U ( α,β ) N ( t + ∆ t ) , ˜ β ) ˜ f ← frequency indicator ( U ( α, ˜ β ) N ) end while else if f < f then β ˜ β ← β/q U ( α, ˜ β ) N ← scale ( U ( α,β ) N ( t + ∆ t ) , ˜ β ) ˜ f ← frequency indicator ( U ( α, ˜ β ) N ) while ˜ f ≤ f and ˜ β ≤ β do β ← ˜ β U ( α,β ) N ( t + ∆ t ) ← U ( α, ˜ β ) N f ← ˜ f f ← ˜ f ˜ β ← β/q U ( α, ˜ β ) N ← scale ( U ( α,β ) N ( t + ∆ t ) , ˜ β ) ˜ f ← frequency indicator ( U ( α, ˜ β ) N ) end while end if able 2: Error, β , and N at t = 5 for different η and γ with both p -adaptive and scaling techniques. γ η technique µ . The scale and move in Fig. 3 are the scaling and moving subroutinesand exterior error indicator calculates the exterior-error indicator for the movingsubroutine. Detailed discussions about scaling and moving techniques are given in [8].After first applying the moving technique, adjusting the expansion order and scal-ing both depend on the frequency indicator and aim to keep the frequency indicator lowto control the error. The relationship and interdependence between them is key to under-standing and justifying the first-scaling-then-adjusting expansion-order procedure in Fig. 3.Thus, we need to investigate how the proposed scaling technique will affect our p -adaptivetechnique and how these two techniques interact with each other. We use two examplescontaining both diffusive and oscillatory behavior to investigate how the two techniqueswill be activated and influence each other. In Example 3, both refinement and reducing β are needed for matching increasing oscillatory and diffusive behavior of the solution; inExample 4, a less oscillatory and diffusive solution over time implies that coarsening andincreasing β may be considered. Example 3.
We approximate the function u ( x, t ) = exp (cid:20) − x ( bt + a ) (cid:21) cos x, t ∈ R + (3.2)with the generalized Laguerre function basis { ˆ L ( α,β ) i ( x ) } ∞ i =0 discussed in [8] with the pa-rameter α = 0. The magnitude of oscillations for this function, exp( − x ( bt + a ) ), increasesover time, requiring proper scaling. Under a variable transformation y = xbt + a , u ( x, t ) canbe rewritten as u ( y, t ) = cos (( bt + a ) y ) exp( − y ), indicating that the solution is increas-ingly oscillatory in y as time increases. Thus, if we reduce the scaling factor β to matchthe diffusive behavior of the solution, proper refinement is also required. In other words,diffusive and oscillatory behavior is coupled in this example. We carry out numerical ex-periments using the algorithm described in Fig. 3 with different ( η, γ ) to investigate howscaling and refinement influence each other. We deactivate the moving technique by set-ting d max = 0 since the solution exhibits no intrinsic advection. Even if we had allowed14 able 3: Error and N at t = 5 for different η and γ with the p -adaptive technique but without thescaling technique, β = 4. γ η moving, it was hardly activated. We set ∆ t = 10 − , N = 50 at t = 0 and a = 2 , b = 0 . q = v − = 0 . , β = 0 . , β = 10 , N min = 0 , N max = 3 and choose the initial scaling factor β = 4.In Table 2 and Table 3 the error in R + is recorded in the lower-left part of each entrywhile the scaling factor β and expansion order N at t = 5 is recorded in the upper-right.By comparing entries in each column/row for smaller η, γ , both tables show the expansionorder N is likely to be increased more when the threshold for refinement ηf is lower.It can be observed from Table 2 that with more refinement β tends to be smaller.This interaction between p -adaptivity and scaling arises because more refinement leads toa larger expansion order N and a smaller scaling threshold νf . Since scaling will only beperformed if the frequency indicator after scaling decreases, proper refinement is not likelyto lead to over-scaling.Moreover, by comparing N at t = 5 between Tables 2 and 3, we see that N tends tobe smaller with the scaling technique for the same γ, η . This implies that without scaling,the refinement procedure is more often activated, leading to a larger N to compensate forthe incapability of scaling alone to maintain a low frequency indicator. This results ina larger computational burden without an improvement in accuracy. This behavior hasbeen expected from the design of Alg. 3 since we put scaling before refinement so thatredistribution of collocation points is tried first to avoid unnecessary refinement when theincrease in frequency indicator results from diffusion instead of oscillation. Example 4.
We approximate the function u ( x, t ) = exp [ − ( bt + a ) x ] cos x, x, t ∈ R + (3.3)with the generalized Laguerre function basis with the parameter α = 0. The magnitudeof oscillations for this function, exp( − ( bt + a ) x ), decreases over time and increasing thescaling factor β to more densely redistribute the allocation points is needed. Furthermore,under the variable transformation y = ( bt + a ) x , u ( x, t ) can be rewritten as u ( y, t ) =cos( ybt + a ) exp( − y ). Since the oscillations decrease with y , one can reduce the expansion15 able 4: Error, β and N at t = 5 for different η and γ with/without scaling for the p -adaptivetechnique. η order. We consider coarsening with or without scaling to investigate whether increasing β can facilitate coarsening (and save computational effort) or result in higher accuracy.We carry out numerical experiments using the algorithm described in Fig. 3 and different( η, γ ) and also deactivate the moving technique by setting d max = 0 since the solutionexhibits no intrinsic advection. We set ∆ t = 10 − , N = 50 at the beginning and set theparameters a = , b = 0 . q = v − = 0 . , β = 0 . , β = 10 , N min = 0 , N max = 3 andinitial scaling factor β = 4. We use a different threshold η for coarsening and we havechecked numerically that the parameter γ in the refinement subroutine will not affect thecoarsening subroutine in this example.In Table 4 the error in R + is recorded in the lower-left part of each entry while the scalingfactor β and expansion order N at t = 5 is recorded in the upper-right. By comparingentries in each row we see that a smaller η will lead to easier coarsening and a smaller N at t = 5. Since the approximation with larger N is always better, whether we can achieve thesame level of accuracy with a smaller expansion order N if proper scaling is implemented isof interest. The initial approximation error is 1 . × − and the approximation will notworsen after coarsening regardless of η because in the p -adaptive subroutine coarseningis allowed only when the post-coarsening frequency indicator remains below the previousthreshold f . Moreover, by comparing the two rows in Table 4 we see that if the solutionconcentrates and becomes less diffusive, increasing β and more efficiently redistributing theallocation points allows the scaling technique to achieve high accuracy with fewer expansionorders than without the scaling technique.The errors and expansion orders over time are plotted in Fig. 4 where the p -adaptivemethod is compared with the non- p -adaptive method when scaling is applied. FromFigs. 4(a) and (c) we can observe that both scaled and unscaled methods maintain theerror below the initial approximation error. Yet, upon comparing Fig. 4(b) to Fig. 4(d)it is readily seen that the scaled method leads to appropriate coarsening while succeedingin maintaining low error, but the unscaled method will increase N when increasing theexpansion order is not actually needed, resulting in additional unnecessary computationalburden. In Fig. 4(e) the scaled and p -adaptive spectral method with η = 4 is comparedwith the scaling-only spectral method. We see that the errors for both methods are almostthe same but the p -adaptive method can reduce unnecessary computation by decreasing N adaptively while still maintaining a low error, and the approximation error for the p -adaptive method fluctuates due to a decreasing N .16 -15 -10 -5 -15 -10 -5 -15 -10 -5 Figure 4:
Approximation to Eq. (3.3) with scaling and p -adaptive spectral methods. Increasing β byscaling can save computational burden while maintaining accuracy by more efficiently redistributingallocation points. The approximation error is controlled below the initial approximation error forboth scaled and unscaled p -adaptive methods, but the expansion order of the scaled method issmaller. On the other hand, adjusting the frequency indicator without decreasing N will notachieve higher accuracy even with a much larger expansion order. β is increased more in the p -adaptive method,implying that the reason why the p -adaptive method can achieve the same accuracy as non- p -adaptive method with a smaller expansion order is that it can redistribute the allocationpoints more efficiently.Finally, we conclude that all three methods: scaling, p -adaptive+scaling, and p -adaptivemethods can maintain the error well below the initial approximation error, but the com-bined p -adaptive+scaling method can achieve this with the smallest expansion order andis therefore the most efficient method among them.
4. Applications in solving the Schr¨odinger equation
In this section, we apply our adaptive spectral methods described in Fig. 3 to solve theSchr¨odinger equation in unbounded domains i∂ t ψ ( x, t ) = − ∂ x ψ ( x, t ) + V ( x ) ψ ( x, t ) + V ex ( x, t ) ψ ( x, t ) , x ∈ R , (4.1)which is equivalent to the PDE discussed in [14] i∂ t u ( x, t ) = (cid:2) − ( ∂ x + iA ( x, t )) + V ( x, t ) (cid:3) u ( x, t ) (4.2)under the transformation u ( x, t ) = e i R t V ex ( x,s ) d s ψ ( x, t ). Here, we shall use spectral meth-ods with the Hermite function basis. The solution is complex, so in the spectral decom-position the coefficients of the basis functions are complex. The major difference here isthat in [14] the Schr¨odinger equation is solved in a bounded domain ( x − , x + ) with absorb-ing boundary conditions. Using spectral methods, we are able to solve the Schr¨odingerequation without truncating the domain.We solve the weak form of Eq. (4.1)( ∂ t ψ, v ) = − i ( ∂ x ψ, ∂ x v ) + (( V ( x ) + V ex ( x, t )) ψ, v ) , v ∈ L ( −∞ , ∞ ) , (4.3)which is to find Ψ βN,x L ( t, x ) := P Ni =0 ψ βi,x L ( t ) ˆ H βi ( x − x L ) in V βN,x L = span { ˆ H βi ( x − x L ) } Ni =0 satisfying the initial condition and( ∂ t Ψ βN,x L , v ) + i ( ∂ x Ψ βN,x L , ∂ x v ) = − i (( V ( x ) + V ex ( x, t ))Ψ N,x L , v ) , ∀ v ∈ V βN,x L . (4.4)We denote ψ βN,x L ( t ) := ( ψ β ,x L ( t ) , ..., ψ βN,x L ( t )), which can be analytically solved to advancetime ψ βN,x L ( t n +1 ) = exp (cid:20) i Z t n +1 t n ( D βN + V βN,x L ( t ))d t (cid:21) ψ βN,x L ( t n ) (4.5)18here D βN ∈ R ( N +1) × ( N +1) is a symmetric matrix with entries( D βN ) ℓj = β p ℓ ( ℓ + 1) j = ℓ + 2 ,β p ( ℓ − ℓ − j = ℓ − ,β ℓ j = ℓ, , (4.6)and the matrix V βN,x L ( t ) ∈ R ( N +1) × ( N +1) has entries( V βN,x L ( t )) ℓj = Z ∞−∞ ( V ( x ) + V ex ( x, t )) ˆ H βℓ − ( x − x L ) ˆ H βj − ( x − x L )d x. (4.7)The evaluation of exp( i R t n +1 t n ( D βN + V βN,x L ( t ))d t ) ψ βN,x L ( t n ) is performed as follows. First,we denote ˜ V βN,x L ≈ R t n +1 t n V βN,x L ( t )d t where the integration is evaluated by Gauss-Legendreformula. Therefore, when calculating the matrix-vector product ˜ V βN,x L X N for a vector X N := ( X , ..., X N ) ∈ R N +1 , its ℓ th component is( ˜ V N,x L X N ) ℓ = N X j =0 N X s =0 ˆ H βℓ − ( x βs ) ˆ H βj ( x βs ) (cid:20) V ( x βs + x L ) + V ex ( x βs + x L , t n + (1 − q )d t )+ V ex ( x βs + x L , t n + d t ) + V ex ( x βs + x L , t n + (1 + q )d t ) (cid:21) X j ∆ t (4.8)where ∆ t = t n +1 − t n . We can first calculate N X j =0 ˆ H βj ( x βs ) (cid:20) V ( x βs + x L ) + V ex ( x βs + x L , t n + (1 − q )d t )+ V ex ( x βs + x L , t n + d2 ) + V ex ( x βs + x L , t n + (1 + q )d t ) (cid:21) X j ∆ t (4.9)for each subindex s ; then, evaluating ( ˜ V βN,x L X N ) ℓ for each subindex ℓ will only require an O ( N ) operation. In this way, given any arbitrary potentials V ( x ) , V ex ( x, t ) we can calculate˜ V βN,x L X N in O ( N ) operations without explicitly calculating entries in ˜ V βN,x L . We approxi-mate the matrix-vector product exp h i ( D βN ∆ t + ˜ V N,x L ) i ψ βN,x L ( t n ) in the following way: werewrite exp h i ( D N ∆ t + ˜ V N,x L ) i ψ βN,x L ( t n ) = exp h m i ( D N ∆ t + ˜ V N,x L ) i m ψ βN,x L ( t n ), whichis introduced as the “scaling and squaring” method in [19], and approximate the matrix-19ector product exp h m i ( D N ∆ t + ˜ V N,x L ) i X N by truncating the infinite Taylor expansionseries P ∞ j =0 1 m j j ! h i ( D N ∆ t + ˜ V N,x L ) i j X N . Here, we take m = 6.As mentioned in Section 1, two main numerical difficulties when solving the Schr¨odingerequation are the unboundedness and oscillatory behavior of the solutions. In fact, thesolution may be increasingly oscillatory behavior at infinity over time, making it very hardto solve in the unbounded domain. However, with our adaptive spectral methods, we canefficiently solve the Schr¨odinger equation in unbounded domains accurately and capturethese oscillations.We shall revisit the two numerical examples discussed in [14]. In the following examples,curves labeled “adaptive” indicate that scaling, moving, and p -adaptive techniques are allapplied as in the algorithm described in Fig. 3, while curves labeled “ p -adaptive” mean thatAlg. 1 is used without scaling or moving; similarly, curves labeled “scaling & p -adaptive”are obtained when the p -adaptive and scaling techniques are applied. Curves labeled“scaling & moving” or “scaling” correspond to applying the scaling and moving techniquesor the scaling technique without the p -adaptive subroutine. The “non-adaptive” curves areobtained when we do not apply any of the scaling, moving, or p -adaptive techniques. Example 5.
We numerically solve the Schr¨odinger equation which is solved in Example1 of [14] and take V = V ex = 0 in Eq. (4.1), admitting the analytic solutionΨ( x, t ) = 1 √ ζ + it exp (cid:20) ik ( x − kt ) − ( x − kt ) ζ + it ) (cid:21) , (4.10)where k is related to the propagation speed of the beam and ζ determines the widthof the beam. The absolute values of the real part of Ψ( x, t = 0 , . ,
1) are plotted inFig. 5(a), illustrating the increasingly oscillatory and diffusive behavior in the rightwardpropagating solution. Treatment of this solution will thus require scaling, moving, and p -adaptive techniques. The imaginary parts of the reference solution (not plotted) overtime are also increasingly oscillatory. We shall apply the algorithm described in Fig. 3.We set ζ = 0 . , k = 1, and initialize N = 50 at t = 0. Other parameters are set to q = ν − = 0 . , µ = 1 . , d = 0 . , β = 0 . , β = 2 , d max = 0 . , N max = 6 , N min =0 , η = 1 . , γ = 1 .
05, and ∆ t = 0 . ψ βN ( t n +1 ) = exp( iD βN d t ) ψ βN ( t n ) . (4.11)When all four techniques are applied, the error is the smallest (shown in Fig. 5(b)) sincewe can keep the exterior error indicator in ( x R , ∞ ) small (shown in Fig. 5(c)) by matchingthe solution’s intrinsic advection. We can simultaneously prevent the frequency indicatorfrom growing too fast (shown in Fig. 5(d)), thus ensuring a small error bound.From the reference solution it can be observed that increasing the expansion order overtime is an intrinsic requirement and failure to do so prevents the capture of the increasing20
15 -9 -3 3 9 1510 -40 -30 -20 -10 -10 -5
10 0.2 0.4 0.6 0.8 110 -5 -10 -5 Figure 5:
Numerically solving the Schr¨odinger equation with vanishing potentials. Applying scaling,moving, and p -adaptive techniques can successfully capture diffusive, advective, and oscillatorybehavior of the solution and yields an accurate numerical solution that prevents the frequencyindicator from growing too fast. The exterior-error indicator is also kept small by moving the basisfunctions rightward to avoid a deteriorating approximation at ∞ . Failure to incorporate any of themoving, scaling, or p -adaptive techniques results in a much larger error. x → ∞ , moving the basis rightward requires correspondingly more refinement (shown inFig. 5(e)). However, the p -adaptive method alone cannot compensate for the inability tocapture diffusion and advection, resulting in an inaccurate approximation. We have alsochecked that apart from what is shown in Fig. 5, applying any single scaling, moving or p -adaptive technique, or combining any two of them will all result in a much larger errorthan employing all three techniques indicated in Fig. 3.Finally, we numerically solve the Schr¨odinger equation with non-vanishing potentials. Example 6.
We numerically solve the following standard Schr¨odinger equation Eq. (4.1)equivalent to Example 2 in [14] with potentials V ex ( x, t ) = 50 √ π sin(10 t ) Z x −∞ exp( − z )d z, V ( x ) = − h e − x − + e − x +1) i . (4.12)Given an even function as the initial condition for Example 2 in [14], the solution is alsoan even function and the solution of Eq. (4.1) obeys | ψ ( − x, t ) | = | ψ ( x, t ) | . No bias towards −∞ or + ∞ is preferred. Therefore, we use the Hermite function basis and apply thealgorithm described in Fig. 3 but deactivate the moving technique by setting d max = 0.We set the initial condition to be the same as that of Example 5 and set η = 1 . , γ =1 , q = 0 . , ν = q − , N min = 0 , N = 200 , β = 0 . , β = 2 and β = 1 . t = 0 with themaximal expansion order increment for each step N max = 20.The reason why we set γ = 1 is that the expansion order N needs to be increased quicklyto catch up with the highly increasingly oscillatory behavior of the numerical solution. Weset a uniform timestep ∆ t = 0 .
01. We use the numerical solution solved with a fixed N = 2500 and only the scaling technique activated as the reference solution. For the p -adaptive method, we have added an additional restriction that the expansion order cannotsurpass the expansion order of the reference solution N = 2500.We can easily see that the spectral method with both scaling and p -adaptive techniquesoutperforms the non-adaptive spectral method or with only one of these two techniquesemployed (shown in Fig. 6(a)). The frequency indicator of using both scaling and p -adaptive techniques is also the smallest (Fig. 6(b)), and the similarity between the frequencyindicator and error is again confirmed as stated in [8]. Moreover, the unscaled method willresult in a larger expansion order (shown in Fig. 6(a)), leading to excessive refinementwith no improvement in accuracy (shown in Fig. 6(a)). In this example, the coarseningprocedure will not lead to a large increase of frequency indicator and does not significantlycompromise accuracy (shown in Figs. 6(b) and (c)). Finally, the scaling factors of the p -adaptive spectral method and the non- p -adaptive spectral method trend similarly overtime; they both decrease after experiencing an initial, transient increase (Fig. 6(d)).22 -10 -5 -10 -5
10 1 2 3 4 520060010001400180022002600 0 1 2 3 4 50.40.81.21.62
Figure 6:
Numerically solving the Schr¨odinger equation with non-vanishing potentials. Rapidlyincreasing oscillations of the solution over time requires much refinement and proper scaling tomaintain accuracy. It is again verified that proper scaling can avoid unnecessary refinement andavoid unnecessary computational burden by adaptively adjusting the scaling factor. Without scal-ing, the expansion order soon reaches the upper bound for N (the expansion order of the referencesolution) and the approximation soon deteriorates due to an inability to further increase N or ad-just β and maintain a low frequency indicator. Failure to accommodate the p -adaptive techniquewill also result in a larger error because of an inability to capture the oscillatory behavior.
5. Summary and Conclusion
In this paper, we proposed a frequency-dependent p -adaptive technique that adjusts theexpansion order for spectral methods. We demonstrated its applicability to time-dependentproblems with varying oscillatory behavior. In order to develop efficient numerical methodsfor problems requiring solutions in unbounded domains, we also combined the p -adaptivetechnique with scaling ( r -adaptivity) and moving ( h -adaptivity) methods to devise a com-plete adaptive spectral method that can successfully deal with diffusion, advection, andoscillation. 23he relationship between scaling and p -adaptive techniques for spectral methods inunbounded domains, both of which depend on the same frequency indicator, is also in-vestigated. We successfully applied our adaptive spectral method to numerically solveSchr¨odinger’s equation. The associated solutions are highly oscillatory in the whole domain,posing numerical difficulties for existing numerical methods that truncate the domain.For future research, the relationship among the adaptive techniques for spectral meth-ods, scaling, moving, refinement and coarsening, can be further studied and rigorous nu-merical analysis for these techniques should be investigated. Furthermore, fast algorithmswith mapped Chebyshev polynomials for solving PDEs in unbounded domains have beendeveloped using the fast Fourier transform [20], but there lacks fast and efficient algorithmsexploiting Laguerre and Hermite basis functions, particularly for higher-dimensional prob-lems. Thus, generalizing these adaptive methods for mapped Jacobi polynomials may bea compelling future research direction. Acknowledgments
MX and TC acknowledge support from the National Science Foundation through grantDMS-1814364 and the Army Research Office through grant W911NF-18-1-0345. SS ac-knowledges financial support from the National Natural Science Foundation of China(Nos. 11822102, 11421101) and Beijing Academy of Artificial Intelligence (BAAI). Compu-tational resources were provided by the High-performance Computing Platform at PekingUniversity.