Time-dependent spectral renormalization method

TL;DR

This paper introduces a novel time-dependent spectral renormalization method for simulating linear and nonlinear evolution equations by converting them into integral equations and solving via fixed-point iteration, incorporating physical properties like conservation laws.

Contribution

It extends the spectral renormalization method to the time domain, enabling efficient simulation of evolution equations with physical constraints.

Findings

Successfully applied to NLS, PT-symmetric nonlocal NLS, and viscous Burgers' equations.

Demonstrates convergence and stability of the method.

Incorporates physics through conservation laws and dissipation rates.

Abstract

The spectral renormalization method was introduced by Ablowitz and Musslimani in 2005, [Opt. Lett. 30, pp. 2140-2142] as an effective way to numerically compute (time-independent) bound states for certain nonlinear boundary value problems. % of the nonlinear Schr\"odinger (NLS), Gross-Pitaevskii and water wave type equations to mention a few. In this paper, we extend those ideas to the time domain and introduce a time-dependent spectral renormalization method as a numerical means to simulate linear and nonlinear evolution equations. The essence of the method is to convert the underlying evolution equation from its partial or ordinary differential form (using Duhamel's principle) into an integral equation. The solution sought is then viewed as a fixed point in both space and time. The resulting integral equation is then numerically solved using a simple renormalized fixed-point iteration method. Convergence is achieved by introducing a time-dependent renormalization factor which is numerically computed from the physical properties of the governing evolution equation. The proposed method has the ability to incorporate physics into the simulations in the form of conservation laws or dissipation rates. This novel scheme is implemented on benchmark evolution equations: the classical nonlinear Schr\"odinger (NLS), integrable $PT$ symmetric nonlocal NLS and the viscous Burgers' equations, each of which being a prototypical example of a conservative and dissipative dynamical system. Numerical implementation and algorithm performance are also discussed.

Tables

Table 1. Table 1: List of renormalization factors and their corresponding conserved quantity.

	$R^2\left(t\right)$	Conserved Quantity
1.	$\frac{C_1}{{\Vert \varphi \Vert }_2^2}$	$J_1$
2.	$\frac{C_2}{\text{Im}\int 𝑑x\varphi {\varphi }_x^{\ast }}$	$J_2$
3.	$\frac{{\Vert {\varphi }_x\Vert }_2^2\pm \sqrt{{\Vert {\varphi }_x\Vert }_2^4-4C_3{\Vert \varphi \Vert }_4^4}}{2{\Vert \varphi \Vert }_4^4}$	$J_3$
4.	$\frac{C_1+C_2}{{\Vert \varphi \Vert }_2^2+\text{Im}\int 𝑑x\varphi {\varphi }_x^{\ast }}$	$J_1+J_2$
	$\frac{\kappa \pm \sqrt{{\kappa }^2-4{\alpha }_3{\Vert \varphi \Vert }_4^4𝒞}}{2{\alpha }_3{\Vert \varphi \Vert }_4^4},𝒞={\alpha }_1{C}_1+{\alpha }_2{C}_2+{\alpha }_3{C}_3,$	${\alpha }_1{J}_1+{\alpha }_2{J}_2+{\alpha }_3{J}_3$
5.	$\kappa ={\alpha }_1{\Vert \varphi \Vert }_2^2+{\alpha }_2\text{Im}\int 𝑑x\varphi {\varphi }_x^{\ast }+{\alpha }_3{\Vert {\varphi }_x\Vert }_2^2$	${\alpha }_3\ne 0$
6.	$\sqrt{\frac{C_1{C}_2}{{\Vert \varphi \Vert }_2^2\text{Im}\int 𝑑x\varphi {\varphi }_x^{\ast }}}$	$J_1{J}_2$
7.	$\frac{C_1{\Vert {\varphi }_x\Vert }_2^2-C_3{\Vert \varphi \Vert }_2^2}{C_1{\Vert \varphi \Vert }_4^4}$	$J_3/J_1$
8.	$\frac{C_2{\Vert {\varphi }_x\Vert }_2^2-C_3\text{Im}\int 𝑑x\varphi {\varphi }_x^{\ast }}{C_2{\Vert \varphi \Vert }_4^4}$	$J_3/J_2$

Equations

$$\frac{\partial \psi }{\partial t}=\mathcal{L}\psi +\mathcal{N}[\psi ],$$

$${\int }_{\Omega }{Q}_j[\psi (x,t)]dx=C_j,$$

$$\frac{d}{dt}{\int }_{\Omega }{X}_j[\psi (x,t)]dx={\int }_{\Omega }{Y}_j[\psi (x,t)]dx,$$

$$\mathcal{L}=\sum _{n=1}^N{a}_n\frac{{\partial }^n}{\partial x^n},$$

$$\psi (x,t=0)=f(x),$$

$$\hat{\psi }(k,t)=F[\psi (x,t)]\equiv \frac{1}{\sqrt{2\pi }}{\int }_{-\infty }^{+\infty }\psi (x,t)e^{-ikx}dx,$$

$$\psi (x,t)=F^{-1}[\hat{\psi }(k,t)]\equiv \frac{1}{\sqrt{2\pi }}{\int }_{-\infty }^{+\infty }\hat{\psi }(k,t)e^{ikx}dk.$$

$$\psi (x,t)=S(t)f(x)+{\int }_0^{t}d\tau S(t-\tau )\mathcal{N}\left[\psi (x,\tau )\right],$$

$$S(t)\equiv e^{t\mathcal{L}},t\ge 0.$$

$$S(t)\eta (x)=F^{-1}\left[e^{t\hat{\mathcal{L}}(k)}\hat{\eta }(k)\right],$$

$$\hat{\mathcal{L}}(k)\equiv \sum _{n=1}^N{a}_n(ik{)}^n.$$

$$\psi (x,t)=R(t)\varphi (x,t).$$

$$\varphi (x,t)$$

$${\varphi }_{n+1}(x,t)$$

$${\int }_{\Omega }{Q}_j[R(t)\varphi (x,t)]dx=C_j.$$

$$\frac{d}{dt}{\int }_{\Omega }{X}_j[R(t)\varphi (x,t)]dx={\int }_{\Omega }{Y}_j[R(t)\varphi (x,t)]dx.$$

$$G(x,\tau )\equiv \mathcal{N}[R(\tau )\varphi (x,\tau )].$$

$$I(x,t)\equiv {\int }_0^{t}d\tau S(t-\tau )G(x,\tau ),$$

$$\hat{I}(k,t)={\int }_0^{t}d\tau e^{(t-\tau )\hat{\mathcal{L}}(k)}\hat{G}(k,\tau ).$$

$$\hat{I}(k,t_{m+1})$$

$$\hat{G}(k,\tau )$$

$$\hat{I}(k,t_{m+1})=$$

$$\displaystyle+A\hat{G}(k,t_{m})+B\hat{G}(k,t_{m+1})\bigg{]},$$

$$A=\frac{e^{-\Delta t\hat{\mathcal{L}}(k)}+\Delta t\hat{\mathcal{L}}(k)-1}{\Delta t{\hat{\mathcal{L}}}^2(k)},$$

$$B=\frac{1-e^{-\Delta t\hat{\mathcal{L}}(k)}\left(\Delta t\hat{\mathcal{L}}(k)+1\right)}{\Delta t{\hat{\mathcal{L}}}^2(k)}.$$

$${\varphi }_{n+1}(x,t)=\frac{1}{R_n(t)}S(t)f(x)+\frac{1}{R_n(t)}{I}_n(x,t),$$

$${\hat{\varphi }}_{n+1}(k,t)=\frac{1}{R_n(t)}{e}^{t\hat{\mathcal{L}}(k)}\hat{f}(k)+\frac{1}{R_n(t)}{\hat{I}}_n(k,t),$$

$${\psi }_{n+1}(x,t)=S(t)f(x)+{\int }_0^{t}d\tau S(t-\tau )\mathcal{N}\left[{\psi }_n(x,\tau )\right],$$

$${\psi }_t=i{\psi }_{xx}+2i|\psi {|}^2\psi ,$$

$${\psi }_{\mathrm{ex}}(x,t)=\mathrm{sech}(x+2\xi t)e^{-i[\xi x+({\xi }^2-1)t]},$$

$$\text{power:}{J}_1\left(\psi \right)\equiv {\int }_{-\infty }^{\infty }dx|\psi {|}^2=C_1,$$

Full text

††journal: Physica D

1 Introduction

Computational methods play an indispensable role in various branches of the physical [1], chemical and biological [2, 3] sciences. In many cases the underlying phenomenon being studied is modeled by either a single or a set of ordinary or partial differential equations. Well-known examples include the quantum Schrödinger equation governing the time evolution of a quantum wave function [4], reaction-diffusion type systems [5] that describe chemical reaction or population dynamics (such as swimming microorganism), and the Navier-Stokes equation modeling the motion of an incompressible Newtonian fluid [6] to mention a few.

Over the years, numerical tools have played an ever increasing role in the advancement of scientific discoveries. They provide a unique opportunity to tackle many challenging scientific problems that are otherwise difficult to solve. This is the case, for instance, in complex turbulent flows, weather prediction, stochastic neural dynamics and many body physics with large degrees of freedom. Surprisingly enough, computational methods and numerical simulations can actually ignite new fundamental ideas and ultimately lead to the development of new theories. One example that stands out is the concept of a soliton which emerged as a result of numerical experiments. Solitons, shape-invariant nonlinear waves that exhibit particle-like behavior upon collision, were discovered in 1965 by Zabusky and Kruskal [7] while performing numerical simulations on the Korteweg-de Vries equation [8]. Their findings sparked intense research interest in many areas of physics and mathematics which subsequently led to the establishment of the inverse scattering transform and integrable nonlinear evolution equations [9, 10]. Soon thereafter, the notion of soliton or solitary wave, spread to many diverse areas in the physical, chemical and biological sciences. Examples include, optical spatial and temporal solitons [11, 12], atomic Bose-Einstein condensates [13, 14], atomic chains [15], molecular and biophysical systems [16] and electrical lattices [17].

A unifying theme among such diverse fields is the development of computational methods capable of accurately and efficiently capturing the physics under study. As such, numerous numerical tools have been developed in the last few decades to simulate evolution equations of the nonlinear Schrödinger and Gross-Pitaevskii type [18, 19, 20, 21, 22, 23, 24]. Perhaps the simplest schemes are finite difference methods such as the Crank-Nicolson algorithm [20, 21, 22]. Among the easiest to implement are the well-known Runge-Kutta integrators [23]. For evolution equations that can be written as a coupled system of linear and nonlinear equations, both of which can be solved exactly, time-splitting methods are fast and typically quite accurate [18, 24]. For stiff problems the so-called exponential differencing methods [25, 26] are effective. Here, the underlying evolution equation is written in an integral form which is then solved using exponential time-stepping schemes [27, 28, 29]. An issue of great importance is how to devise numerical schemes capable of enforcing physics into the simulations. For example, when it comes to simulating a conservative dynamical system one faces the challenge of imposing conservation laws. In the time-independent setting multiple methods are known that yield stationary soliton solutions that conserve total power [30, 31]. If the system under study happens to be Hamiltonian, a commonly used method is the geometric or symplectic integrators [32, 33, 34]. They are capable of exactly preserving the symplectic area in phase space and the Hamiltonian.

In this paper we introduce the time-dependent spectral renormalization method as an effective and simple computational tool to numerically simulate linear and nonlinear evolution equations. The proposed algorithm has the capability to incorporate the underlying laws of physics in the form of conservation laws or dissipation/rate equations. The idea is to convert the given dynamical system from its evolution equation into an integral form. This approach allows us to think of the solution sought as being a fixed point in space and time of the resulting integral equation. To compute the fixed point a time-dependent renormalization factor is introduced which is found using the physical properties of the original governing evolution equation such as conservation laws or dissipation rates. The solution is then obtained from a renormalized fixed-point iteration scheme. This novel time-dependent spectral renormalization approach allows the flexibility to “integrate” physics into the simulations. The proposed method is applied on benchmark problems: the classical nonlinear Schrödinger, integrable $PT$ symmetric nonlocal nonlinear Schrödinger and viscous Burgers’ equations, each being a prototypical representative of conservative and dissipative evolution equations. Our scheme extendes the time-independent spectral renormalization concept introduced by Ablowitz and Musslimani in 2005 [35] to the time domain. We point out that the steady state spectral renormalization has been widely used in many applications related to nonlinear optics [36, 37, 38, 39], Bose-Einstein condensation [40] and water waves [41].

The paper is organized as follows. In Sec. 2 we introduce the time-dependent spectral renormalization scheme and detail its implementation to general evolution equations. The numerical implementation of the algorithm to conservative systems is given in Sec. 3 where the classical and integrable $PT$ symmetric NLS equations are used to study important properties of the method. Application of the scheme to dissipative evolution equations is presented in Sec. 4 where the Burgers’ equation is used as a benchmark problem to highlight the generalities of the proposed scheme. We conclude in Sec. 5.

2 Time-dependent spectral renormalization

To begin, we consider a complex or real valued wave function $\psi(x,t)$ that depends on space variable $x$ and time $t\geq 0.$ The time evolution of $\psi$ is given by

[TABLE]

where $\mathcal{L}$ is a linear differential operator and $\mathcal{N}$ is a nonlinear function of $\psi.$ The spatial domain $\Omega$ over which the partial differential equation (1) is defined can be either bounded or unbounded. Like most physically relevant evolution equations, Eq. (1) represents either a conservative or dissipative dynamical system. In the former case, this implies the existence of conserved quantities given by

[TABLE]

where all $C_{j},j=1,2,3,\cdots$ are constant in time. In the latter case, a set of dissipation or rate equations induced from the dynamical system (1) are derived in the form

[TABLE]

where $X_{j}$ and $Y_{j},j=1,2,3,\cdots$ are functions of $\psi(x,t)$ and referred to as the density and flux, respectively. To highlight the effectiveness of our method and emphasize its usefulness, we shall limit the discussion to linear differential operators $\mathcal{L}$ with constant coefficients given by

[TABLE]

where all $a_{n}$ are constants in space and time. The treatment of a more general case for which all or some of the coefficients $a_{n}$ are allowed to depend on $x$ and/or $t$ is also possible. In this situation, any such term(s) would be included in the nonlinearity $\mathcal{N}[\psi].$ Equation (1) is supplemented with the initial condition

[TABLE]

and either periodic, $\psi(x+L,t)=\psi(x,t),$ with $L$ being the domain size, or rapidly decaying to zero boundary conditions. Throughout the rest of the paper, we shall use the symmetric forward and inverse Fourier transforms defined by

[TABLE]

[TABLE]

For problems formulated on a bounded spatial interval $\Omega=[0,L]$ we instead represent the field $\psi(x,t)$ in terms of its Fourier series representation. The implementation of the time-dependent spectral renormalization algorithm to a general evolution equation of the type given in (1) follows several simple steps which we next outline in detail.

Rewrite Eq. (1) in Duhamel’s integral form:

[TABLE]

with $S(t)$ being the so-called time propagator (or semi-group) defined by

[TABLE]

Equation (8) is often used as a departing point for the derivation of many time-stepping schemes such as Runge-Kutta [42, 43, 44, 45] and exponential differencing methods [28], among others. Importantly, formula (8) implies that the solution $\psi(x,t)$ can be viewed as a fixed point in space and time. The propagator $S(t)$ can be computed with the aid of the Fourier transform or Fourier series. Thus, we have the following representation

[TABLE]

where $\eta(x)$ is an arbitrary square-integrable or periodic function and $\hat{\mathcal{L}}(k)$ is the Fourier symbol corresponding to the linear operator $\mathcal{L}$ defined by (here $i=\sqrt{-1}$)

[TABLE] 2. 2.

Introduce a time-dependent renormalization factor $R(t)$ by rewriting the function as

[TABLE]

Here, $\phi(x,t)$ is either a real or complex valued function depending on the nature of the partial differential equation (1). At this stage, $R(t)$ is an unknown function of time and is assumed to be real and nonzero at all times. A remark about complex renormalization is discussed in Sec. 3 in connection to the classical and $PT$ symmetric NLS equations. Substituting (12) into Eq. (8) gives an integral representation for the new field

[TABLE]

Equation (2) constitutes the basis for the spectral renormalization method. The solution $\phi(x,t)$ is numerically found from the fixed-point iteration

[TABLE]

for $n=1,2,3,\cdots$. This iteration is seeded with an initial guess $\phi_{1}(x,t).$ 3. 3.

Compute the renormalization factor $R(t)$ from the associated conservation law(s) or dissipation (rate) equation(s) induced from dynamical system (1). Specifically speaking, we have the following:

• (a)

If Eq. (1) is conservative, then the renormalization factor is computed from (2), i.e.,

[TABLE]

• (b)

If Eq. (1) is dissipative, then the renormalization factor is derived from the rate equation (3) which gives an evolution equation for $R(t)$ in the form

[TABLE]

The calculation of $R(t)$ is exemplified in Secs. 3 and 4 when applied to the NLS and viscous Burgers’ equations. This in turn allows the simulation “to keep in touch” with the original evolution equation and incorporate relevant physics into the integrator. 4. 4.

Evaluate the time integral that appears in Eq. (2). For ease of presentation we define the quantity

[TABLE]

Thus, the integral we are interested in computing is given by

[TABLE]

which, in Fourier space, takes the form

[TABLE]

To numerically evaluate the integral (19) we divide the time interval $[0,T]$ with $T$ being the end evolution time into $M$ (even) equally spaced intervals each of size $\Delta t\equiv T/M.$ Then, denote the time levels by $t_{j}\equiv j\Delta t,j=0,1,2,3,\cdots,M.$ Clearly, $t_{0}=0$ (initial time) and $t_{M}=T$ (end evolution time). After some algebra, we find that $\hat{I}(k,t_{m})$ satisfies the recursion relation

[TABLE]

valid for $m=0,1,2,\cdots,M-1.$ Our goal next is to derive an approximate formula for the integral that appears in Eq. (4). This can be accomplished by replacing the integrand $\hat{G}(k,\tau)$ by either a constant, linear, quadratic or higher order polynomial of the time variable $\tau$. This is referred to as Filon integration [46]. See also [25, 28] for applications of the Filon method in exponential time differencing. In this paper, we present results only for linear interpolants. On each time interval $t_{m}\leq\tau\leq t_{m+1},m=0,1,2,\cdots,M-1,$ we interpolate the function $\hat{G}(k,\tau)$ by the linear function

[TABLE]

Substituting Eq. (4) into (4) and performing integration by parts, we arrive at the result

[TABLE]

for $m=0,1,2,\cdots,M-1$ with $\hat{I}(k,t_{0})=\hat{I}(k,0)\equiv 0.$ Here, we define

[TABLE]

[TABLE]

Since the coefficients $A$ and $B$ depend only on $\Delta t$ and the Fourier mode $k$ (but not on the solution itself or the iteration index $m$), they are pre-computed only once. For certain problems, such as $\hat{\mathcal{L}}(k)=-k^{2}$ at $k=0$, the operator $\hat{\mathcal{L}}(k)$ vanishes at some value(s) of $k$. In this situation formulas (23) and (24) are undefined. To remedy this, their values are replaced by their corresponding limits $A\rightarrow\Delta t/2$ and $B\rightarrow\Delta t/2$ as $\hat{\mathcal{L}}(k)\rightarrow 0$. This idea has already been implemented in the context of exponential time differencing [28]. Another way to avoid dividing by zero is to use the Cauchy integral formula [29]. 5. 5.

Implement the time-dependent spectral renormalization fixed-point iteration based on Eqs. (2) and (18). Below we write it both in physical and Fourier spaces respectively:

[TABLE]

[TABLE]

where $n=1,2,3,\cdots.$ At every iteration step $n$ the function $I(x,t)$ or alternatively $\hat{I}(k,t)$ is obtained from Eq.  (22).

The iteration schemes (25) or (26) are seeded with an “arbitrary” initial guess $\phi_{1}(x,t)$ that satisfies the boundary conditions. Typical examples include: (i) the initial condition (5), (ii) $f(x)h(t)$ with $h(t)$ an “arbitrary” function of time and (iii) a random function in space and time. Convergence is achieved when the relative error between successive iterations is less than a prescribed level of error tolerance $\epsilon$, i.e., $\max_{x,t}|\phi_{n+1}-\phi_{n}|\leq\epsilon$ as $n\rightarrow\infty$. We remark that the time-dependent spectral renormalization approach allows one to choose the size of the total time interval $T$ to be large enough such that the standard (unrenormalized, or $R(t)=1$) fixed point iteration

[TABLE]

fails to converge. So far, we have presented a general framework where a solution to an evolution equation of the type given in Eq. (1) can be numerically obtained by a renormalized time-dependent fixed-point iteration. In the next two sections we implement this scheme on two important and physically relevant examples that represent conservative and dissipative dynamical systems.

3 Conservative case

3.1 Classical NLS equation

As a prototypical example representing a conservative dynamical system, we consider the classical nonlinear Schrödinger equation [47]

[TABLE]

where $\psi=\psi(x,t)$ is a complex-valued function of the real variables $x$ and time $t\geq 0.$ Subscripts indicate partial derivatives with respect to $x$ and $t.$ Equation (28) is posed on the whole real line and is supplemented with the following boundary and initial conditions respectively: $\psi(x,t)$ goes to zero sufficiently fast as $|x|\rightarrow\infty$ and $\psi(x,t=0)=f(x)$ which is assumed to be smooth and square-integrable on the whole real line. There are two important reasons to consider the NLS equation as a test bed for the performance of the time-dependent spectral renormalization scheme. The first is due to its wide applications in optics, condensed matter physics (such as Bose-Einstein condensation, superfluidity and superconductivity) and fluid mechanics (deep water waves). The second reason is tied to the fact that Eq.  (28) is an integrable evolution equation and admits an infinite number of conserved quantities [9]. Furthermore, it admits the exact moving soliton solution

[TABLE]

where $\xi$ is a real constant used to generate momentum. In the case of more general NLS-type equations several spectrally accurate numerical methods exist to find soliton modes [48, 49]. For the NLS equation, we have the following three physically relevant conserved quantities:

[TABLE]

where all $C_{j}$ are constants of motion. Next, we list the necessary steps to implement the time-dependent spectral renormalization method. To this end, we have $\mathcal{L}=i\frac{\partial^{2}}{\partial x^{2}},\hat{\mathcal{L}}(k)=-ik^{2},\mathcal{N}=2i|\psi|^{2}\psi$, the propagator (or semigroup) defined in formula (9) given by $S(t)=\exp\left(it\frac{\partial^{2}}{\partial x^{2}}\right),$ and the renormalized nonlinear term $G(x,t)=2iR^{3}(t)|\phi(x,t)|^{2}\phi(x,t).$ Solutions to Eq. (28), corresponding to initial condition (5) are then computed from the iterative scheme (26) with the renormalization factor $R(t)$ taken from Table 1 where

[TABLE]

Each choice of $R$ would correspond to a certain conservation law.

It is obvious that the “space of options” to single out one “favorite” formula for the renormalization factor is rather large. The decision as to which expression to pick is based on various factors such as convenience, the physics of the problem and the numerical validity of that specific formula.

We have implemented the time-dependent spectral renormalization scheme on the NLS equation corresponding to three different types of initial conditions: (i) $f(x)={\rm sech}(x)$ which gives rise to a stationary soliton, (ii) $f(x)={\rm sech}(x)e^{-2ix}$ leading to a moving pure soliton carrying nonzero momentum and (iii) a Gaussian wavepacket $f(x)=e^{-x^{2}/2}.$ The constants $C_{j},j=1,2,3$ are calculated directly from the initial condition.

Several iterations of the solution intensity $|\psi_{n}(x,t)|^{2}$ are shown in Fig. 1 which is obtained from iterating Eqn. (26) using $R(t)$ that exactly conserves power (row 1 in Table 1). To examine the effectiveness of the method, the iteration is seeded with four different initial guesses. In the first experiment, we chose an initial guess identical to the NLS initial condition, i.e., $\phi_{1}(x,t)={\rm sech}(x).$ This choice is somewhat natural to begin with. After a few rounds the wave function has corrected itself and locked on to the soliton solution given in Eq. (29) with $\xi=0$ – see Fig. 1(a). To test the scheme’s sensitivity to initial guesses, we kept the same initial condition as before but now seeded the iteration with a random field in both space and time chosen to be uniformly distributed on the interval $[0,1]\times[0,1].$ Interestingly enough, the outcome of the run was independent of the initial input as it again converged to the numerical soliton solution. This robust behavior is depicted in Fig. 1(b). We next conduct some numerical tests using different initial conditions. Several simulations are shown in Figs. 1(c) and (d) for a traveling soliton and Gaussian solutions, respectively. While the results shown in this paper are carried out using conservation of power only (as is the case in Fig. 1) we point out that other conservation quantities have been used to converge to the same solution. All numerical tests have been performed on a large spatial domain $(L=100)$ to insure the problem is indeed posed on the whole real line and guarantee decaying (zero) boundary conditions at all times. The total time interval used to produce Fig. 1 is $T=10$ which, without the use of any renormalization technique, would cause the iterative method in Eq. (27) to diverge.

One interesting behavior we observe is that the renormalization factor approaches zero when dealing with large time intervals. This is the case, for example, when computing stationary solitons $(\xi=0)$ while using real $R(t).$ Figure 2 (a), (b) and (c) depicts such behavior. However, this is not the case when solving the NLS equation with Gaussian input. Here, the (real) renormalization factor remained on the same order as it started for all times – see Figure 2 (d). We point out that the spectral renormalization scheme converges even when one chooses a rather large time grid spacing $\Delta t.$ Figure 3 shows maximum error between successive iterations defined by $\max_{x,t}|\psi_{n+1}(x,t)-\psi_{n}(x,t)|$ as a function of the iteration number $n.$ For small $\Delta t$ the scheme converged after some $25$ iterates giving rise to an accurate numerical solution. Interestingly, the iteration converged even for a relatively large $\Delta t$ (on the order of 0.5), however, the resulting solution shows poor accuracy.

To quantify the performance of the scheme, such as convergence, we have measured at each simulation run the maximum spatial difference between the (convergent) numerically obtained solution $\psi_{\rm num}(x,t)$ and the exact soliton given in (29). Thus, we define

[TABLE]

Furthermore, at each simulation run, we have monitored the difference between each numerically computed conserved quantity and its initial value as a function of time. Precisely, we checked the time evolution of the following quantities:

[TABLE]

[TABLE]

[TABLE]

For a fixed grid size $\Delta t$ we ran three different types of simulations all using the same NLS initial condition $f(x)={\rm sech}(x)e^{-2ix}$ ($\xi=2$), but different $R(t)$ that are computed from either conservation of power, momentum, or Hamiltonian. Once convergence is achieved, the output $\psi_{{\rm num}}(x,t)$ is recorded and used to compute the maximum value of the quantities $\mathsf{E}(t),\mathsf{P}(t),\mathsf{M}(t)$ and $\mathsf{H}(t)$ given in Eqns.(34)–(37) over the entire time domain. We then doubled the number of grid points (halved the grid size $\Delta t$) and repeated the same experiment again. This procedure was performed eight times (corresponding to eight different $\Delta t$). The findings are summarized in Fig. 4. As one can see from Fig. 4(a), the maximum error in the solution is converging to the exact soliton and is doing so at a second-order rate. Furthermore, Fig. 4(b), (c) and (d) show the errors in the power, momentum and Hamiltonian to be near machine precision when computed using their respective renormalization factor alone for each $\Delta t$. Otherwise, the error does not seem to be exceptional.

Finally, we have monitored the time evolution of the quantities $\mathsf{E}(t),\mathsf{P}(t),\mathsf{M}(t)$ and $\mathsf{H}(t)$ given in Eqns. (34)–(37). They are displayed in Fig. 5 for a traveling wave soliton. The error in the computed solution compared to the exact one, $E(t)$, is found to grow with time. Importantly, the errors in power, momentum and Hamiltonian are machine precision accurate when computed using their respective renormalization factor (see Fig. 5).

3.2 Comparison to other methods

It is of considerable interest to compare the time-dependent spectral renormalization method to other well known time-integrators. For that purpose, we choose to solve the NLS equation (28) with initial condition $f(x)={\rm sech}(x)e^{-2ix}$ using two different second-order accurate schemes: split-step Fourier and exponential time-differencing Runge-Kutta (ETDRK2). In our spectral renormalization simulations we enforced conservation of power only and computed the renormalization factor $R(t)$ from row 1 in Table 1. To quantify the performance of each scheme, we repeat the calculations that were performed in Fig. 4 using the three schemes mentioned above. In Fig. 6(a) the solution errors obtained from all three integrators are observed to be on the same order of magnitude (with the split-step error slightly lower than the others) and converging at the same rate. The story is different when it comes to conservation of power. Here, both the spectral renormalization and split-step methods are found to be numerically exact (see in Fig. 6(b)). Note that the split-step method and certain finite-difference schemes are known to conserve power exactly [50, 51]. The error obtained using these two methods are favorable in comparison to exponential time-differencing scheme. For the parameters considered here, the spectral renormalization method yields power values that are around seven orders of magnitude more accurate. On the other hand, the momentum and Hamiltonian are most accurate when computed using the split-step method. We point out that for the ETDRK2 scheme considered here, the physical quantities that any solution must satisfy appear to be constrained by the accuracy of the numerical solution. This is clearly not the case for the spectral renormalization scheme which can enforce certain relevant physical properties to excellent accuracy regardless of how accurate the solution is.

In comparison to standard method of lines type integrators, the spectral renormalization algorithm requires considerable storage for large $t.$ Specifically, one must store $(M+1)N$ data points, where $N,M$ denote the number of spatial and temporal grid points. In two spatial dimensions this storage requirement becomes $(M+1)N^{2}$. To alleviate this storage burden, as well as to increase the algorithm speed, one possibility is to divide the full time interval $[0,T]$ into $\mathcal{M}$ sub-intervals $[0,T]=[0,T_{1}]\cup[T_{1},T_{2}]\cdots\cup[T_{M-1},T_{M}]$. This cuts the required memory inside each sub-interval to $(M+1)N/\mathcal{M}$. Each time interval length is chosen sufficiently large such that the spectral renormalization algorithm is effective, efficient and fast. On the time span $[0,T_{1}]$ the proposed algorithm is implemented using the initial condition $\psi(x,0)=f(x)$ and continues until convergence to a fixed point in space and time is reached. The scheme is then applied to the time interval $[T_{1},T_{2}]$ with initial condition $f(x)=\psi(x,t=T_{1}).$ The code is repeated until the whole time span is covered. Typically, for problems integrated on time interval $[0,100]$ ten partitions would be sufficient for good performance. Finally, using the fast Fourier transform to approximate derivatives results in CPU run times that are respectable. For example, to reach the converged solutions shown in Fig. 1 on a standard laptop using Matlab took about one second. Different initial guesses resulted in similar run times.

3.3 Complex Renormalization Factor

So far, we have implemented and examined the performance of the time-dependent spectral renormalization scheme under the assumption that $R(t)$ is real valued. In fact this restriction is not needed for the formulation of the scheme. It is rather rooted in the fact that all NLS conservation laws give the magnitude but not the phase of the renormalization factor. To remedy this issue, we propose here one possibility of lifting this constraint by deriving an expression for the renormalization factor directly from the NLS equation. With this approach, one is still able to incorporate conservation laws. To this end, we substitute $\psi(x,t)=R(t)\phi(x,t)$ into Eq. (28), integrate over the entire real line while utilizing the fast decay of the wave function (localized boundary conditions) and obtain

[TABLE]

where

[TABLE]

[TABLE]

The constant $\alpha(0)$ is computed from (39) whereas $R(0)$ is found from Eq. (12) and (5). The result is

[TABLE]

Formula (38) provides an alternative expression for the renormalization factor $R$ without any assumption of being real. However, it is valid so long the quantity given in (39) is not zero. Furthermore, this approach restricts the choices for $\phi(x,0)$ as it could cause $\alpha(0)$ to become zero. Last, but not least, since $\phi$ is in general complex function, the term inside the exponent in Eq. (38) can lead to exponential growth or decay which will eventually put serious strain on the time numerical integration. To incorporate physics into (38) we replace $|R(\tau)|^{2}$ that appears inside the exponent by any of the expressions given in Table 1 where $R^{2}(\tau)$ is now replaced by $|R(\tau)|^{2}.$ We have implemented the time-dependent spectral renormalization scheme for the complex $R$ case and found that it converged to the solution from a large selection of initial guesses. The main difference between the real and complex renormalization approach is that now $R$ seems not to go to zero; it rather oscillates in time (see Fig. 2).

3.4 $PT$ symmetric integrable nonlocal NLS equation

In this section we apply the time-dependent spectral renormalization method to the $PT$ symmetric nonlocal nonlinear Schrödinger equation

[TABLE]

where $\psi(x,t)$ is a complex valued function of the real variables $x\in\mathbb{R}$ and time $t\geq 0.$ Equation (42) was first introduced in [69] and shown to be an integrable infinite-dimensional Hamiltonian dynamical system [70, 71, 72]. As such, it admits an infinite number of conserved quantities. Of particular interest is the conservation of the so-called “quasi-power”

[TABLE]

Furthermore, in [69] a one soliton solution was obtained in the form of a breathing pure one soliton

[TABLE]

where $\eta,\bar{\eta}$ are positive constants. Interest in wave propagation in $PT$ symmetric media has been at the forefront of research in physics and mathematics [52, 53, 54, 55, 56, 60, 58, 57, 61, 62, 63, 64, 65, 67, 68, 59, 66]. Here, we show how to use and implement the renormalization scheme to obtain a time-periodic soliton solution. Solutions to Eq. (42) corresponding to initial condition (5) are then computed from the iterative scheme given in Eq. (26). The renormalization factor $R(t)$ is computed from the conservation law (43) and is given by (assuming it is real)

[TABLE]

Note that the quantity in the denominator is real valued but not necessarily of a definite sign.

An alternative formula for the renormalization that assumes it be complex can be derived in a manner similar to what we presented for the classical NLS equation. To do so, we substitute $\psi(x,t)=R(t)\phi(x,t)$ into the nonlocal NLS equation and, after some algebra, find

[TABLE]

where $\alpha(\tau)$ is defined in (39) and

[TABLE]

[TABLE]

The value of $R(0)$ can be obtained by multiplying Eq. (5) by $f^{*}(-x)$ and integrating over the entire spatial domain

[TABLE]

We have implemented the time-dependent spectral renormalization scheme to solve the initial boundary value associated with the nonlocal $PT$ symmetric NLS Eq. (42) corresponding to initial condition given by

[TABLE]

The iteration scheme is seeded with an initial guess $\phi_{1}(x,t)=f(x).$ After some number of iteration, the method converged and locked on the breathing soliton given in (44). In order to quantify these numerical results we compare our converged solution to the exact one using the global max error given in (34). Additionally, since we derive the renormalization factor directly from conservation of quasi-power (45) we also measure the error in the quasi-power via the quantity

[TABLE]

A summary of our findings using the renormalization factor defined in Eq. (45) is shown in Fig. 7. For this set of parameters, namely $\eta\not=\bar{\eta}$, the soliton mode approaches a singularity at time $t_{S}=\pi/\left[4(\eta^{2}-\bar{\eta}^{2})\right]\approx 3.74$. The solution intensity is shown in Fig. 7(a) where the soliton peak is rapidly growing in time. We next compare the numerical solutions obtained from our scheme with other second-order accurate time-integration techniques. Figure 7(b) reveals that all the methods exhibit second-order convergence. Particularly, in Fig. 7(c) we display the error in the quasi-power (here given by $\Gamma_{1}=4.2$) obtained from three different methods. The split-step and spectral renormalization methods are found to preserve the quasi-power exactly, whereas the error computed using the ETDRK2 method is many orders of magnitude larger. Lastly, we compare the real (45) and complex (46) versions of the renormalization factor. The two different forms of $R(t)$ are shown in Fig. 8 for identical parameters. As with the classical NLS above, when $R$ is real the magnitude rapidly decays to zero, while when $R$ is complex it is observed to oscillate in time.

4 Dissipative Case: Burger’s equation

So far we have explained and implemented the time-dependent spectral renormalization algorithm in conservative systems and demonstrated its flexibility in terms of incorporating conserved physical properties into the simulation. In this section, we turn our attention to dissipative evolution equations and apply the scheme on the viscous Burgers’ equation

[TABLE]

Here, $\psi=\psi(x,t)$ is a real valued function and $\nu>0$ is the viscosity coefficient. The problem will be considered on the bounded domain $[0,2\pi]$ with periodic boundary conditions: $\psi(x+2\pi,t)=\psi(x,t)$ and $\psi_{x}(x+2\pi,t)=\psi_{x}(x,t).$ The integration time interval is $[0,T].$ The limit as $\nu\rightarrow 0$ is referred to as the inviscid Burgers’ equation. It was first proposed by Burgers [73] as a model for turbulent flows and later shown to play a fundamental role in many branches of the nonlinear sciences. This equation exhibits rich mathematical structures such as finite time blow-up (singularity), shock formation and wave breaking [47]. Physically speaking, it finds wide applications in gas dynamics, acoustics and traffic flows to name a few [74]. The viscous Burgers’ ($\nu\neq 0$) is often used in fluid mechanics as a simplified version of the Navier-Stokes equation in one space dimension under the assumption of incompressibility and no pressure gradient. Importantly, Eq. (52) can be solved exactly with the help of the so-called Cole-Hopf transformation [75, 76]. Indeed, if $\psi(x,t)$ solves Eq. (52) then the function $u(x,t)$ defined by

[TABLE]

satisfies the linear heat equation

[TABLE]

A spatially periodic exact solution to Eq. (52) corresponding to the initial condition

[TABLE]

is given by

[TABLE]

The viscous Burgers’ equation is dissipative in nature. Indeed, multiplying Eq. (52) by $\psi$ and integrating over the entire spatial domain we find

[TABLE]

This relation describes the rate at which the total energy of the system is dissipated.

With all this at hand, the initial boundary value problem (52) corresponding to initial condition (55) is solved using the iterative scheme given in Eq. (26). For the viscous Burgers’ equation we have the linear operator $\hat{\mathcal{L}}(k)=-\nu k^{2}$ and nonlinearity $\mathcal{N}=-\left(\psi^{2}\right)_{x}/2$. To derive a formula for the renormalization factor, we substitute $\psi(x,t)=R(t)\phi(x,t)$ into (57); integrate over the entire spatial domain and obtain

[TABLE]

where

[TABLE]

[TABLE]

Solving Eq. (58) for nontrivial solutions we obtain

[TABLE]

Unlike the NLS case, the quantities $\theta_{j}(t),j=1,2$ are non-negative and real. To evaluate the cumulative time integral in Eq. (61) a second-order accurate trapezoidal method is used.

Since we are directly enforcing dissipation of energy in Eq. (57) we are interested to learn how the energy in the numerical solution compares to that of the exact one. For that purpose, we define the quantity

[TABLE]

where $\psi_{\rm num}(x,t)$ is the numerical solution obtained from iterating Eq. (26) using the initial condition (55). The error convergence rates are shown in Fig. 9. Again, the global solution error in Fig. 9(a) is observed to converge at a second-order rate. This is also the case for the ETDRK2 scheme. Next, the error in the (time-dependent) dissipation rate is compared with the exact. Doing so reveals that, unlike the time-independent conserved quantities above, the dissipation rate is not numerically exact, but instead converges at second-order rate [see Fig. 9(b)]. The reason for this is that the time integral in Eq. (61) is approximated by a second-order accurate quadrature method. A similar convergence rate is found for the Runge-Kutta scheme.

5 Conclusions

In this paper, we have presented a new numerical method to simulate evolution equations which we term here the time-dependent spectral renormalization scheme. The proposed scheme is rooted in the (time-independent) spectral renormalization method introduced in 2005 by Ablowitz and Musslimani [35] to numerically compute stationary nonlinear bound states for certain nonlinear boundary value problems. In this regard, the proposed method can be thought of as an extension to the time domain. The idea here is to rewrite the given evolution equation (formulated in an ordinary or partial differential equation form) into an integral equation. The solution is then viewed as a fixed point in both space and time, rather than being a solution to an evolution equation. The resulting integral equation is then numerically solved using a simple renormalized fixed-point iteration. Convergence of the method is achieved by introducing a time-dependent renormalization factor which is numerically computed from the physical properties of the governing evolution equations. This novel time-dependent spectral renormalization scheme has the ability to incorporate physics into the numerical simulations and allows the numerical time integration to “keep in touch” with the original evolution equation. The proposed method is applied to two benchmark problems: the nonlinear Schrödinger and the Burgers’ equations each of which being a prototypical example of a conservative and dissipative system respectively.

6 Acknowledgment

The work of J.T.C. and Z.H.M was supported in part by NSF grant number DMS-0908599. Z.H.M thanks the Lady-Davis Trust for the financial support during his visit to the Technion-Israel Institute of Technology.