Efficient quantum algorithm for dissipative nonlinear differential equations

Edited by Anthony Leggett, University of Illinois at Urbana–Champaign, Urbana, IL, and approved July 19, 2021 (received for review March 6, 2021)
August 26, 2021
118 (35) e2026805118

Significance

Nonlinear differential equations appear in many domains and are notoriously difficult to solve. Whereas previous quantum algorithms for general nonlinear differential equations have complexity exponential in the evolution time, we give the first quantum algorithm for dissipative nonlinear differential equations that is efficient provided the dissipation is sufficiently strong relative to nonlinear and forcing terms and the solution does not decay too rapidly. We also establish a lower bound showing that differential equations with sufficiently weak dissipation have worst-case complexity exponential in time, giving an almost tight classification of the quantum complexity of simulating nonlinear dynamics. Furthermore, numerical results for the Burgers equation suggest that our algorithm may potentially address complex nonlinear phenomena even in regimes with weaker dissipation.

Abstract

Nonlinear differential equations model diverse phenomena but are notoriously difficult to solve. While there has been extensive previous work on efficient quantum algorithms for linear differential equations, the linearity of quantum mechanics has limited analogous progress for the nonlinear case. Despite this obstacle, we develop a quantum algorithm for dissipative quadratic n-dimensional ordinary differential equations. Assuming R<1, where R is a parameter characterizing the ratio of the nonlinearity and forcing to the linear dissipation, this algorithm has complexity T2qpoly(logT,logn,log1/ϵ)/ϵ, where T is the evolution time, ϵ is the allowed error, and q measures decay of the solution. This is an exponential improvement over the best previous quantum algorithms, whose complexity is exponential in T. While exponential decay precludes efficiency, driven equations can avoid this issue despite the presence of dissipation. Our algorithm uses the method of Carleman linearization, for which we give a convergence theorem. This method maps a system of nonlinear differential equations to an infinite-dimensional system of linear differential equations, which we discretize, truncate, and solve using the forward Euler method and the quantum linear system algorithm. We also provide a lower bound on the worst-case complexity of quantum algorithms for general quadratic differential equations, showing that the problem is intractable for R2. Finally, we discuss potential applications, showing that the R<1 condition can be satisfied in realistic epidemiological models and giving numerical evidence that the method may describe a model of fluid dynamics even for larger values of R.
Models governed by both ordinary differential equations (ODEs) and partial differential equations (PDEs) arise extensively in the natural and social sciences, medicine, and engineering. Such equations characterize physical and biological systems that exhibit a wide variety of complex phenomena, including turbulence and chaos. We focus here on differential equations with nonlinearities that can be expressed with quadratic polynomials, which include many archetypal models in biology, fluid dynamics, and plasma physics.
Quantum algorithms offer the prospect of rapidly characterizing the solutions of high-dimensional systems of linear ODEs (13) and PDEs (410). Such algorithms can produce a quantum state proportional to the solution of a sparse (or block-encoded) n-dimensional system of linear differential equations in time poly(logn) using the quantum linear system algorithm (QLSA) (11).
Early work on quantum algorithms for differential equations already considered the nonlinear case (12). This work gave a quantum algorithm for ODEs that simulates polynomial nonlinearities by storing multiple copies of the solution. The complexity of this approach is polynomial in the logarithm of the dimension but exponential in the evolution time, scaling as O(1/ϵT). While some heuristic quantum algorithms for nonlinear ODEs have been studied more recently (13, 14), these works do not make precise statements about concrete implementations or running times of quantum algorithms. The recent preprint (15) also describes a quantum algorithm to solve a nonlinear ODE by linearizing it using a different approach from the one taken here. However, a proof of correctness of their algorithm involving a bound on the condition number and probability of success is not given. The authors also do not describe how barriers such as those of ref. 16 could be avoided in their approach.
Constructing efficient quantum algorithms for general classes of nonlinear dynamics has been considered largely out of reach since the linearity of quantum mechanics makes it challenging to efficiently represent nonlinear dynamics. Furthermore, nonlinear modifications of quantum mechanics generically enable quickly solving hard computational problems such as performing unstructured search among n items in time poly(logn). Indeed, the search lower bound (17) can be used to show that nonlinear dynamics is exponentially difficult to simulate in general (16, 18, 19).
In this article, we design and analyze a quantum algorithm that overcomes this limitation using Carleman linearization (2022). This approach embeds polynomial nonlinearities into an infinite-dimensional system of linear ODEs and then truncates it to obtain a finite-dimensional linear approximation. We discretize the finite ODE system in time using the forward Euler method and solve the resulting linear equations with the QLSA (11, 23). Subject to the condition R<1, where the quantity R (defined in Problem 1 below) characterizes the relative strength of the nonlinear and dissipative linear terms, we show that the total complexity of this quantum Carleman linearization (QCL) algorithm is sT2qpoly(logT,logn,log1/ϵ)/ϵ, where s is the sparsity, T is the evolution time, q is the ratio of the norm of the initial vector to that of the solution vector at time T, n is the dimension, and ϵ is the allowed error (Theorem 1). In the regime R<1, this is an exponential improvement over ref. 12, which has complexity exponential in T.
Note that the solution cannot decay exponentially in T for the algorithm to be efficient, as captured by the dependence of the complexity on q—a known limitation of quantum ODE algorithms (2). For homogeneous ODEs with R<1, the solution necessarily decays exponentially in time (SI Appendix, Eq. S30), so the algorithm is not asymptotically efficient. However, even for solutions with exponential decay, we still find an improvement over the best previous result O(1/ϵT) (12) for sufficiently small ϵ. Thus, our algorithm might provide an advantage over classical computation for studying evolution for short times. More significantly, our algorithm can handle inhomogeneous quadratic ODEs, for which it can remain efficient in the long-time limit since the solution can remain asymptotically nonzero (for an explicit example, see the discussion just before the proof of Lemma 1 in SI Appendix, Section 3.A.1) or can decay slowly [i.e., q can be poly(T)]. Inhomogeneous equations arise in many applications, including, for example, the discretization of PDEs with nontrivial boundary conditions.
We also show that the requirement R<1 cannot be significantly improved in general (Theorem 2). Specifically, we provide a quantum lower bound showing worst-case hardness of simulating strongly nonlinear dynamics for R2, establishing a limitation on the ability of quantum computers to simulate general nonlinear dynamics.
Our quantum algorithm could potentially be applied to study models governed by quadratic ODEs arising in biology and epidemiology as well as in fluid and plasma dynamics. In particular, the celebrated Navier–Stokes equation with linear damping, which describes many physical phenomena, can be treated by our approach if the Reynolds number is sufficiently small. While the formal validity of our arguments assumes R<1, we find in numerical experiments that our proposed approach remains valid for larger R in certain cases.
We emphasize that as in related quantum algorithms for linear algebra and differential equations, instantiating our approach requires an implicit description of the problem that allows for efficient preparation of the initial state and implementation of the dynamics. Furthermore, since the output is encoded in a quantum state, readout is restricted to features that can be revealed by efficient quantum measurements. More work remains to understand how these methods might be applied, as we discuss further in Discussion.

Problem Statement

We focus on an initial value problem described by the n-dimensional quadratic ODE
dudt=F2u2+F1u+F0(t),  u(0)=uin.
[1]
Here u=[u1,,un]TRn, u2=[u12,u1u2,,u1un,u2u1,,unun1,un2]TRn2, each uj=uj(t) is a function of t on the interval [0,T] for j[n]{1,,n}, F2Rn×n2, F1Rn×n are time-independent matrices, and the inhomogeneity F0(t)Rn is a C1 continuous function of t.
The main computational problem we consider is as follows.
Problem 1. In the quantum quadratic ODE problem, we consider an n-dimensional quadratic ODE as in [1]. We assume F2, F1, and F0 are s-sparse (i.e., have at most s nonzero entries in each row and column), F1 is diagonalizable, and the eigenvalues λj of F1 satisfy Re(λn)Re(λ1)<0. We parametrize the problem in terms of the quantity
R1|Re(λ1)|uinF2+F0uin.
[2]
For some given T>0, we assume the values Re(λ1), F2, F1, F0(t) for each t[0,T], and F0maxt[0,T]F0(t), F0maxt[0,T]F0(t) are known, and that we are given oracles OF2, OF1, and OF0 that provide the locations and values of the nonzero entries of F2, F1, and F0(t) for any specified t, respectively, for any desired row or column. We are also given the value uin and an oracle Ox that maps 000Cn to a quantum state proportional to uin. Our goal is to produce a quantum state proportional to u(T) for some given T>0 within some prescribed 2 error tolerance ϵ>0.
When F0(t)=0 (i.e., the ODE is homogeneous), the quantity R=uinF2|Re(λ1)| characterizes the ratio of the nonlinearity to the linear dissipation. It is qualitatively similar to the Reynolds number in fluid dynamics.
Our algorithm uses a rescaling of Problem 1 that we introduce in SI Appendix, Section 1.

Main Results

For Problem 1 with R<1, we present a quantum algorithm that improves the complexity from an exponential dependence on T to a nearly quadratic dependence. This algorithm is based on the concept of Carleman linearization (2022). Carleman linearization is a method for converting a finite-dimensional system of nonlinear differential equations into an infinite-dimensional linear one. This is achieved by introducing powers of the variables into the system, allowing it to be written as an infinite sequence of coupled linear differential equations. We then truncate the system to N equations, where the truncation level N depends on the allowed approximation error, giving a finite linear ODE system.
Given a system of quadratic ODEs [1], we apply the Carleman procedure to obtain the system of linear ODEs
dŷdt=A(t)ŷ+b(t),  ŷ(0)=ŷin
[3]
with the tridiagonal block structure
ddtŷ1ŷ2ŷN=A11A21A12A22ANN1AN1NANNŷ1ŷ2ŷN+F0(t)00,
[4]
where ŷj=ujRnj, ŷin=[uin;uin2;;uinN], and Aj+1jRnj×nj+1, AjjRnj×nj, Aj1jRnj×nj1 for j[N] satisfying
Aj+1j=F2Ij1+IF2Ij2++Ij1F2,
[5]
Ajj=F1Ij1+IF1Ij2++Ij1F1,
[6]
Aj+1j=F0(t)Ij1+IF0(t)Ij2++Ij1F0(t).
[7]
To construct a system of linear equations, we next divide [0,T] into m=T/h time steps and apply the forward Euler method on [3], letting
yk+1=[I+A(kh)h]yk+b(kh)
[8]
where ykRn+n2++nN approximates ŷ(kh) for each k[m+1]0{0,1,,m}, with y0=yinŷ(0)=ŷin, and letting all yk be equal for k[m+p+1]0\[m+1]0, for some sufficiently large integer p. (It is unclear whether another discretization could improve performance, as discussed further in Discussion.)
In the above system, the first n components of yk for k[m+p+1]0 (i.e., y1k) is ϵ-close to the exact solution u(T), up to normalization. We apply the high-precision QLSA (23) to [8] and postselect on k to produce y1k/y1k for some k[m+p+1]0\[m]0. SI Appendix, Section 2, gives more details about the system of linear equations and the measurement procedure.
We control the additive approximation error ϵ by combining a convergence theorem for Carleman linearization that improves previous analysis (22) with a bound for the global error of the Euler method. To analyze the complexity of this approach, we upper bound the condition number of the linear system and lower bound the success probability of the final measurement. We state our main algorithmic result as follows.
Theorem 1. (Informal version.) Consider an instance of the quantum quadratic ODE problem as defined in Problem 1, with its Carleman linearization as defined in [3]. Assume R<1. Let quin/u(T). There exists a quantum algorithm producing a state that approximates u(T)/u(T) with error at most ϵ1, succeeding with probability Ω(1), with query and gate complexity
sT2qϵpoly(logT,logn,log1/ϵ).
[9]
SI Appendix, Section 2, states the formal version of Theorem 1 with more details.
We also provide a quantum lower bound for the worst-case complexity of simulating strongly nonlinear dynamics, showing that the problem is intractable for R2. Following the approach of refs. 16, 18, we construct a protocol for distinguishing two states of a qubit driven by a certain quadratic ODE. Provided R2, this procedure distinguishes states with overlap 1ϵ in time poly(log(1/ϵ)). Since nonorthogonal quantum states are hard to distinguish, this implies a lower bound on the complexity of the quantum ODE problem. This gives the following limitation on quantum algorithms for nonlinear ODEs.
Theorem 2. Assume R2. Then there is an instance of the quantum quadratic ODE problem defined in Problem 1 such that any quantum algorithm for producing a quantum state approximating u(T)/u(T) with bounded error must have worst-case time complexity exponential in T.

Algorithm Analysis

In this section we describe the main ingredients used to prove Theorem 1.
First, we provide an upper bound for the error from Carleman linearization for arbitrary evolution time T. To the best of our knowledge, the first and only explicit bound on the error of Carleman linearization appears in ref. 22. However, the authors only consider homogeneous quadratic ODEs, and to bound the error for arbitrary T, the logarithmic norm of F1 is assumed to be negative (theorems 4.2 and 4.3 of ref. 22), which is too strong for our case. Instead, we give an analysis under milder conditions, providing a convergence guarantee for general inhomogeneous quadratic ODEs.
Lemma 1. Consider an instance of the quadratic ODE [1], with its corresponding Carleman linearization as defined in [3]. As in Problem 1, assume that the eigenvalues λj of F1 satisfy Re(λn)Re(λ1)<0. Then for any j[N], the error ηj(t)uj(t)ŷj(t) satisfies
ηj(t)tNF2uinN+1
[10]
where we require R<1 as defined in [2].
We present the proof of Lemma 1 (and an improved version for homogeneous equations) and explain the necessity of our assumptions in SI Appendix, Section 3.A.1.
Next, we provide an upper bound for the error incurred by approximating [3] with the forward Euler method.
Lemma 2. Consider an instance of the quantum quadratic ODE problem as defined in Problem 1, with R<1 as defined in [2]. Choose a sufficiently small time step h. Suppose the error from Carleman linearization satisfies η(t)u(T)/4. Then the global error from the forward Euler method [8] on the interval [0,T] for [3] satisfies
ŷ1(T)y1(T)3N2.5Th[(F2+F1+F0)2+F0].
[11]
SI Appendix, Section 3.A.2, gives the proof of Lemma 2.
Next we upper bound the condition number of the linear system.
Lemma 3. Consider an instance of the quantum quadratic ODE problem as defined in Problem 1. Apply the forward Euler method [8] to the Carleman linearization [3]. Then the condition number of the linear system satisfies
κ3(m+p+1).
[12]
SI Appendix, Section 3.B, gives the proof of Lemma 3.
We next describe a procedure for preparing the right-hand side of the linear system.
Lemma 4. Assume we are given the value uin, and let Ox be an oracle that maps 000Cn to a normalized quantum state uin proportional to uin. Assume we are also given the values F0(t) for each t[0,T], and let OF0 be an oracle that provides the locations and values of the nonzero entries of F0(t) for any specified t. Then a quantum state encoding the right-hand side of the linear system can be prepared using O(N) queries to Ox and O(m) queries to OF0, with gate complexity larger by a poly(logN,logn) factor.
The proof of Lemma 4 appears in SI Appendix, Section 3.C, along with a discussion of why the encoding of the initial state suffices.
After applying the QLSA to [8], we perform a measurement to extract a final state of the desired form. We now consider the probability that this measurement succeeds.
Lemma 5. Consider an instance of the quantum quadratic ODE problem defined in Problem 1, with the QLSA applied to the linear system using the forward Euler method [8]. Suppose the error from Carleman linearization satisfies η(t)u(T)/4, and the error from Euler method satisfies ŷ(T)ymu(T)/4. Then the probability of measuring a state y1k for k=[m+p+1]0\[m+1]0 satisfies
Pmeasurep+19(m+p+1)Nq2.
[13]
We prove Lemma 5 in SI Appendix, Section 3.D.
The proof of the main algorithmic result (Theorem 1) follows by applying the above lemmas with appropriate parameters and carefully analyzing the complexity. The detailed analysis appears in SI Appendix, Section 3.E.

Lower Bound

In this section, we establish a limitation on the ability of quantum computers to solve the quadratic ODE problem when the nonlinearity is sufficiently strong. We quantify the strength of the nonlinearity in terms of the quantity R defined in [2]. Whereas there is an efficient quantum algorithm for R<1 (as shown in Theorem 1), we show here that the problem is intractable for R2. We establish this result by showing how the nonlinear dynamics can be used to distinguish nonorthogonal quantum states, a task that requires many copies of the given state.
Previous work on the computational power of nonlinear quantum mechanics shows that the ability to distinguish nonorthogonal states can be applied to solve unstructured search (and other hard computational problems) (16, 18, 19). Here we show a similar limitation using an information-theoretic argument.
Lemma 6. Let |ψ,|ϕ be states of a qubit with |ψ|ϕ|=1ϵ. Suppose we are given either a black box that prepares ψ or a black box that prepares ϕ. Then any bounded-error protocol for determining whether the state is ψ or ϕ must use Ω(1/ϵ) queries.
The proof of Lemma 6 is presented in SI Appendix, Section 4.A.
Lemma 6 can be used to establish limitations on the ability of quantum computers to simulate general nonlinear dynamics since nonlinear dynamics can be used to distinguish nonorthogonal states. Whereas previous work considers models of nonlinear quantum dynamics [such as the Weinberg model (18, 19) and the Gross-Pitaevskii equation (16)], here we aim to show the difficulty of efficiently simulating more general nonlinear ODEs—in particular, quadratic ODEs with dissipation—using quantum algorithms.
Lemma 7. There exists an instance of the quantum quadratic ODE problem as defined in Problem 1 with R2, and two states of a qubit with overlap 1ϵ (for 0<ϵ<13/10) as possible initial conditions, such that the two final states after evolution time T=O(log(1/ϵ)) have an overlap no larger than 3/10.
SI Appendix, Section 4.B, presents the proof of Lemma 7.
The proof of our main lower bound result (Theorem 2) follows by combining Lemma 6 and Lemma 7, and is presented in SI Appendix, Section 4.C.
Note that exponential time is achievable since our QCL algorithm can solve the problem by taking N to be exponential in T, where N is the truncation level of Carleman linearization. [The algorithm of Leyton and Osborne also solves quadratic differential equations with complexity exponential in T but requires the additional assumptions that the quadratic polynomial is measure-preserving and Lipschitz continuous (12).]

Applications

Our quantum algorithm could be applied to study models governed by inhomogeneous quadratic ODEs with linear dissipation. Such equations arise in biology and epidemiology as well as in fluid and plasma dynamics.
One concrete example is the so-called SEIR model of an epidemic (24). The model
dPSdt=ΛPSPrvacPS+ΛrtraPSPIP
[14]
dPEdt=ΛPEPPETlat+rtraPSPIP
[15]
dPIdt=ΛPIP+PETlatPITinf
[16]
dPRdt=ΛPRP+rvacPS+PITinf
[17]
describes the population of P individuals composed of four components: susceptible (PS), exposed (PE), infected (PI), and recovered (PR). Here rtra is the rate of transmission from an infected to a susceptible person, rvac is the vaccination rate, Tlat is the latent time until an exposed person becomes infectious, Tinf is the infectious time that an infectious person can infect others, and Λ is a flux of individuals constantly refreshing the population. This type of model has been used to describe the early spread of the COVID-19 virus (24). Realistic parameters from that study give a value R<1 (as discussed in more detail in SI Appendix, Section 5), showing that the assumptions of our algorithm can correspond to some real-world problems that are only moderately nonlinear. This example can be generalized to a high-dimensional system of ODEs that models the early spread over a large number of cities with interaction (25, 26).
We also consider the celebrated Navier–Stokes equation
tu+(u)u+βu=ν2u+f
[18]
with linear damping and a forcing term. Equations of the form [18] are ubiquitous in fluid mechanics (27), and related models such as those studied in refs. 2830 are used to describe the formation of large-scale structure in the universe. Similar equations also appear in magnetohydrodynamics (31) and in models that describe the motion of free particles that stick to each other upon collision (32). In the inviscid case, ν=0, the resulting Euler-type equations with linear damping are also of interest, both for modeling micromechanical devices (33) and for their intimate connection with viscous models (34).
As a specific example, consider the one-dimensional forced viscous Burgers equation
tu+uxu=νx2u+f,
[19]
which is the one-dimensional case of Eq. 18 with β=0. This is often used as a simple model of convective flow (35). For concreteness, let the initial condition be u(x,0)=U0sin(2πx/L0) on the domain x[L0/2,L0/2], and use Dirichlet boundary conditions u(L0/2,0)=u(L0/2,0)=0. We force this equation using a localized off-center Gaussian with a sinusoidal time dependence, given by f(x,t)=U0exp(xL0/4)22(L0/32)2cos(2πt). To solve this equation using the Carleman method, we discretize the spatial domain into nx points and use central differences for the derivatives to get a system of quadratic ODEs. This equation is of the form [1] and can thus generate the Carleman system [4]. The resulting linear ODE can then be integrated using the forward Euler method, as shown in Fig. 1. Even though the parameters used in this example result in R44, which violates the requirement R<1 of the QCL algorithm, we see from the absolute error plot in Fig. 1 that the maximum absolute error over time decreases exponentially as the truncation level N is incremented (in this example, the maximum Carleman truncation level considered is N=4). Surprisingly, this suggests that in this example, the error of the classical Carleman method converges exponentially with N, even though R>1.
Fig. 1.
Integration of the forced viscous Burgers equation using Carleman linearization on a classical computer (source code available in ref. 38). The viscosity is set so that the Reynolds number Re=U0L0/ν=20. The parameters nx=16 and nt=4,000 are the number of spatial and temporal discretization intervals, respectively. The corresponding Carleman convergence parameter is R=43.59. (Left) Initial condition, source, and solutions plotted at a third of the nonlinear time 13Tnl=L03U0. The direct solution is obtained from the ode45 solver in MATLAB. (Center) l2 norm of the absolute error between the Carleman solutions at various truncation levels N and the direct solution. (Right) Convergence of the time-maximum error.

Discussion

In this paper we have presented a QCL algorithm for a class of quadratic nonlinear differential equations. Compared to the previous approach of ref. 12, our algorithm improves the complexity from an exponential dependence on T to a nearly quadratic dependence, under the condition R<1 as defined in [2]. Qualitatively, this means that the system must be dissipative and that the nonlinear and inhomogeneous effects must be small relative to the linear effects. We have also provided numerical results suggesting the classical Carleman method may work on certain PDEs that do not strictly satisfy the assumption R<1. Furthermore, we established a lower bound showing that for general quadratic differential equations with R2, quantum algorithms must have worst-case complexity exponential in T. We also discussed several potential applications arising in biology and fluid and plasma dynamics.
It is natural to ask whether the result of Theorem 1 can be achieved with a classical algorithm, i.e., whether the assumption R<1 makes differential equations classically tractable. Clearly, a naive integration of the truncated Carleman system [4] is not efficient on a classical computer since the system size is Θ(nN). However, furthermore, it is unlikely that any classical algorithm for this problem can run in time polylogarithmic in n. Indeed, if we consider Problem 1 with λ1T1 but let the nonlinearity and forcing be even smaller such that R<1, then in the asymptotic limit we have a linear differential equation with no dissipation. Hence, any classical algorithm that could solve Problem 1 could also solve nondissipative linear differential equations, which is a BQP-hard problem even when the dynamics are unitary (36). In other words, an efficient classical algorithm for this problem would imply efficient classical algorithms for any problem that can be solved efficiently by a quantum computer, which is considered unlikely.
Our upper and lower bounds leave a gap in the range 1R<2, for which we do not know the complexity of the quantum quadratic ODE problem. We hope that future work will close this gap and determine for which R the problem can be solved efficiently by quantum computers in the worst case.
Furthermore, the complexity of our algorithm has nearly quadratic dependence on T, namely, T2poly(logT). It is unknown whether the complexity for quadratic ODEs must be at least linear or quadratic in T. Note that sublinear complexity is impossible in general because of the no-fast-forwarding theorem (37). However, it should be possible to fast-forward the dynamics in special cases, and it would be interesting to understand the extent to which dissipation enables this.
The complexity of our algorithm depends on the parameter q defined in Theorem 1, which characterizes the decay of the final solution relative to the initial condition. This restricts the utility of our result since we must have a suitable initial condition and terminal time such that the final state is not exponentially smaller than the initial state. However, it is unlikely that such a dependence can be significantly improved since renormalization of the state can be used to implement postselection, which would imply the unlikely consequence BQP=PP (see section 8 of ref. 2 for further discussion). As discussed in the introduction, the solution of a homogeneous dissipative equation necessarily decays exponentially in time, so our method is not asymptotically efficient. However, for inhomogeneous equations the final state need not be exponentially smaller than the initial state even in a long-time simulation, suggesting that our algorithm could be especially suitable for models with forcing terms.
The quantum part of the algorithm might also be improved. In this paper we limit ourselves to the first-order Euler method to discretize the linearized ODEs in time. This is crucial for the analysis in Lemma 2, which states that the global error increases at most linearly with T (see the formal version of Lemma 2 in SI Appendix, Section 3.A.2). However, it is unclear how to give a similar bound for higher-order numerical schemes. If this obstacle could be overcome, the error dependence of the complexity might be improved.
It is also natural to ask whether our approach can be improved by taking features of particular systems into account. Since the Carleman method has only received limited attention and has generally been used for purposes other than numerical integration, it seems likely that such improvements are possible. In fact, the numerical results discussed in Applications (see in particular Fig. 1) suggest that the condition R<1 is not a strict requirement for the viscous Burgers equation, since we observe convergence even though R44. This raises the possibility that our approach could be applicable in a broader context than our analytical results might suggest. We leave a detailed investigation of this for future work.
A related question is whether our algorithm can efficiently simulate systems exhibiting dynamical chaos. The condition R<1 might preclude chaos, but we do not have a proof of this. More generally, the presence or absence of chaos might provide a more fine-grained picture of the algorithm’s efficiency.
When contemplating applications, it should be emphasized that our approach produces a state vector that encodes the solution without specifying how information is to be extracted from that state. Simply producing a state vector is not enough for an end-to-end application since the full quantum state cannot be read out efficiently. In some cases it may be possible to extract useful information by sampling a simple observable, whereas in other cases, more sophisticated postprocessing may be required to infer a desired property of the solution. Our method does not address this issue but can be considered as a subroutine whose output will be parsed by subsequent quantum algorithms. We hope that future work will address this issue and develop end-to-end applications of these methods.
Finally, the algorithm presented in this paper might be extended to solve related mathematical problems on quantum computers. Obvious candidates include initial value problems with time-dependent coefficients and boundary value problems. Carleman methods for such problems are explored in ref. 21, but it is not obvious how to implement those methods in a quantum algorithm. It is also possible that suitable formulations of problems in nonlinear optimization or control could be solvable using related techniques.

Data Availability

Source code data have been deposited in GitHub (38).

Acknowledgments

We thank Paola Cappellaro for valuable discussions and comments. J.-P.L. thanks Aaron Ostrander for inspiring discussions. H.Ø.K. thanks Bernhard Paus Graesdal for introducing him to the Carleman method. H.K.K. and N.F.L. thank Seth Lloyd for preliminary discussions on nonlinear equations and the quantum linear systems algorithm. We also thank Dong An and anonymous referees for pointing out the exponential decay of solutions to dissipative homogeneous equations. A.M.C. and J.-P.L. did part of this work while visiting the Simons Institute for the Theory of Computing in Berkeley and gratefully acknowledge its hospitality. A.M.C. and J.-P.L. acknowledge support from the Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Quantum Algorithms Teams, and Accelerated Research in Quantum Computing programs and from the NSF (CCF-1813814). H.Ø.K., H.K.K., and N.F.L. were partially funded by award DE-SC0020264 from the Department of Energy.

Supporting Information

Appendix (PDF)

References

1
D. W. Berry, High-order quantum algorithm for solving linear differential equations. J. Phys. A Math. Theor. 47, 105301 (2014).
2
D. W. Berry, A. M. Childs, A. Ostrander, G. Wang, Quantum algorithm for linear differential equations with exponentially improved dependence on precision. Commun. Math. Phys. 356, 1057–1081 (2017).
3
A. M. Childs, J. P. Liu, Quantum spectral methods for differential equations. Commun. Math. Phys. 375, 1427–1457 (2020).
4
B. D. Clader, B. C. Jacobs, C. R. Sprouse, Preconditioned quantum linear system algorithm. Phys. Rev. Lett. 110, 250504 (2013).
5
Y. Cao, A. Papageorgiou, I. Petras, J. Traub, S. Kais, Quantum algorithm and circuit design solving the Poisson equation. New J. Phys. 15, 013021 (2013).
6
A. Montanaro, S. Pallister, Quantum algorithms and the finite element method. Phys. Rev. A (Coll. Park) 93, 032324 (2016).
7
P. C. S. Costa, S. Jordan, A. Ostrander, Quantum algorithm for simulating the wave equation. Phys. Rev. A (Coll. Park) 99, 012323 (2019).
8
A. M. Childs, J. P. Liu, A. Ostrander, High-precision quantum algorithms for partial differential equations. arXiv [Preprint] (2020) https://arxiv.org/abs/2002.07868 (Accessed 18 February 2020).
9
A. Engel, G. Smith, S. E. Parker, Quantum algorithm for the Vlasov equation. Phys. Rev. A (Coll. Park) 100, 062315 (2019).
10
N. Linden, A. Montanaro, C. Shao, Quantum vs. classical algorithms for solving the heat equation. arXiv [Preprint] (2020). https://arxiv.org/abs/2004.06516 (Accessed 14 April 2020).
11
A. W. Harrow, A. Hassidim, S. Lloyd, Quantum algorithm for linear systems of equations. Phys. Rev. Lett. 103, 150502 (2009).
12
S. K. Leyton, T. J. Osborne, A quantum algorithm to solve nonlinear differential equations. arXiv [Preprint] (2008). https://arxiv.org/abs/0812.4423 (Accessed 5 March 2020).
13
I. Joseph, Koopman-von Neumann approach to quantum simulation of nonlinear classical dynamics. arXiv [Preprint] (2020). https://arxiv.org/abs/2003.09980 (Accessed 22 March 2020).
14
I. Y. Dodin, E. A. Startsev, On applications of quantum computing to plasma simulations. arXiv [Preprint] (2020). https://arxiv.org/abs/2005.14369 (Accessed 29 May 2020).
15
S. Lloyd et al., Quantum algorithm for nonlinear differential equations. arXiv [Preprint] (2020). https://arxiv.org/abs/2011.06571 (Accessed 12 November 2020).
16
A. M. Childs, J. Young, Optimal state discrimination and unstructured search in nonlinear quantum mechanics. Phys. Rev. A (Coll. Park) 93, 022314 (2016).
17
C. H. Bennett, E. Bernstein, G. Brassard, U. Vazirani, Strengths and weaknesses of quantum computing. SIAM J. Comput. 26, 1510–1523 (1997).
18
D. S. Abrams, S. Lloyd, Nonlinear quantum mechanics implies polynomial-time solution for NP-complete and #P problems. Phys. Rev. Lett. 81, 3992 (1998).
19
S. Aaronson, NP-complete problems and physical reality. ACM SIGACT News 36, 30–52 (2005).
20
T. Carleman, Application de la théorie des équations intégrales linéaires aux systèmes d’équations différentielles non linéaires. Acta Math. 59, 63–87 (1932).
21
K. Kowalski, W. H. Steeb, Nonlinear Dynamical Systems and Carleman Linearization (World Scientific, 1991).
22
M. Forets, A. Pouly, Explicit error bounds for Carleman linearization. arXiv [Preprint] (2017). https://arxiv.org/abs/1711.02552 (Accessed 5 March 2020).
23
A. M. Childs, R. Kothari, R. D. Somma, Quantum algorithm for systems of linear equations with exponentially improved dependence on precision. SIAM J. Comput. 46, 1920–1950 (2017).
24
C. Wang et al., Evolving epidemiology and impact of non-pharmaceutical interventions on the outbreak of coronavirus disease 2019 in Wuhan, China. JAMA 323, 1915–1923 (2020).
25
D. Bichara, Y. Kang, C. Castillo-Chávez, R. Horan, C. Perrings, SIS and SIR epidemic models under virtual dispersal. Bull. Math. Biol. 77, 2004–2034 (2015).
26
R. Qurratu Aini, D. Aldila, K. Sugeng, Basic reproduction number of a multi-patch SVI model represented as a star graph topology. AIP Conf. Proc. 2023, 020237 (2018).
27
P. G. Lemarié-Rieusset, The Navier-Stokes Problem in the 21st Century (CRC Press, 2018).
28
L. Kofman, D. Pogosian, S. Shandarin, Structure of the universe in the two-dimensional model of adhesion. Mon. Not. R. Astron. Soc. 242, 200–208 (1990).
29
S. F. Shandarin, Y. B. Zeldovich, The large-scale structure of the universe: Turbulence, intermittency, structures in a self-gravitating medium. Rev. Mod. Phys. 61, 185–220 (1989).
30
M. Vergassola, B. Dubrulle, U. Frisch, A. Noullez, Burgers’ equation, devil’s staircases and the mass distribution for large-scale structures. Astron. Astrophys. 289, 325–356 (1994).
31
P. A. Davidson, An Introduction to Magnetohydrodynamics (Cambridge Texts in Applied Mathematics) (Cambridge University Press, 2001).
32
Y. Brenier, E. Grenier, Sticky particles and scalar conservation laws. SIAM J. Numer. Anal. 35, 2317–2328 (1998).
33
M. H. Bao, Micro Mechanical Transducers: Pressure Sensors, Accelerometers and Gyroscopes (Elsevier, 2000).
34
C. M. Dafermos, Hyperbolic Conservation Laws in Continuum Physics (Springer, 2005), vol. 3.
35
J.M. Burgers, A mathematical model illustrating the theory of turbulence. Adv. Appl. Mech. 1, 171–199 (1948).
36
R. P. Feynman, Quantum mechanical computers. Opt. News 11, 11–20 (1985).
37
D. W. Berry, G. Ahokas, R. Cleve, B. C. Sanders, Efficient quantum algorithms for simulating sparse Hamiltonians. Commun. Math. Phys. 270, 359–371 (2007).
38
H. Ø. Kolden, Carleman solution of the viscous Burgers equation. GitHub. https://github.com/hermankolden/CarlemanBurgers. Deposited 1 March 2021.

Information & Authors

Information

Published in

Go to Proceedings of the National Academy of Sciences
Go to Proceedings of the National Academy of Sciences
Proceedings of the National Academy of Sciences
Vol. 118 | No. 35
August 31, 2021
PubMed: 34446548

Classifications

Data Availability

Source code data have been deposited in GitHub (38).

Submission history

Published online: August 26, 2021
Published in issue: August 31, 2021

Keywords

  1. quantum algorithm
  2. nonlinear differential equations
  3. Carleman linearization
  4. Navier–Stokes equation
  5. plasma dynamics

Acknowledgments

We thank Paola Cappellaro for valuable discussions and comments. J.-P.L. thanks Aaron Ostrander for inspiring discussions. H.Ø.K. thanks Bernhard Paus Graesdal for introducing him to the Carleman method. H.K.K. and N.F.L. thank Seth Lloyd for preliminary discussions on nonlinear equations and the quantum linear systems algorithm. We also thank Dong An and anonymous referees for pointing out the exponential decay of solutions to dissipative homogeneous equations. A.M.C. and J.-P.L. did part of this work while visiting the Simons Institute for the Theory of Computing in Berkeley and gratefully acknowledge its hospitality. A.M.C. and J.-P.L. acknowledge support from the Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Quantum Algorithms Teams, and Accelerated Research in Quantum Computing programs and from the NSF (CCF-1813814). H.Ø.K., H.K.K., and N.F.L. were partially funded by award DE-SC0020264 from the Department of Energy.

Notes

This article is a PNAS Direct Submission.

Authors

Affiliations

Joint Center for Quantum Information and Computer Science, University of Maryland, College Park, MD 20742;
Institute for Advanced Computer Studies, University of Maryland, College Park, MD 20742;
Department of Mathematics, University of Maryland, College Park, MD 20742;
Department of Physics, Norwegian University of Science and Technology, NO-7491 Trondheim, Norway;
Plasma Science and Fusion Center, Massachusetts Institute of Technology, Cambridge, MA 02139;
Hari K. Krovi
Quantum Engineering and Computing, Raytheon BBN Technologies, Cambridge, MA 02138;
Nuno F. Loureiro
Plasma Science and Fusion Center, Massachusetts Institute of Technology, Cambridge, MA 02139;
Konstantina Trivisa
Department of Mathematics, University of Maryland, College Park, MD 20742;
Institute for Physical Science and Technology, University of Maryland, College Park, MD 20742;
Joint Center for Quantum Information and Computer Science, University of Maryland, College Park, MD 20742;
Institute for Advanced Computer Studies, University of Maryland, College Park, MD 20742;
Department of Computer Science, University of Maryland, College Park, MD 20742

Notes

1
To whom correspondence may be addressed. Email: [email protected].
Author contributions: J.-P.L., H.Ø.K., H.K.K., N.F.L., K.T., and A.M.C. designed research; J.-P.L., H.Ø.K., H.K.K., N.F.L., K.T., and A.M.C. performed research; J.P.L. led the theoretical analysis; H.Ø.K. led the numerical experiments; and J.-P.L., H.Ø.K., H.K.K., N.F.L., K.T., and A.M.C. wrote the paper.

Competing Interests

The authors declare no competing interest.

Metrics & Citations

Metrics

Note: The article usage is presented with a three- to four-day delay and will update daily once available. Due to ths delay, usage data will not appear immediately following publication. Citation information is sourced from Crossref Cited-by service.


Citation statements




Altmetrics

Citations

If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click Download.

Cited by

    Loading...

    View Options

    View options

    PDF format

    Download this article as a PDF file

    DOWNLOAD PDF

    Get Access

    Login options

    Check if you have access through your login credentials or your institution to get full access on this article.

    Personal login Institutional Login

    Recommend to a librarian

    Recommend PNAS to a Librarian

    Purchase options

    Purchase this article to get full access to it.

    Single Article Purchase

    Efficient quantum algorithm for dissipative nonlinear differential equations
    Proceedings of the National Academy of Sciences
    • Vol. 118
    • No. 35

    Media

    Figures

    Tables

    Other

    Share

    Share

    Share article link

    Share on social media