# 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)

## 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 $\text{R}<1$, where $\text{R}$ is a parameter characterizing the ratio of the nonlinearity and forcing to the linear dissipation, this algorithm has complexity ${T}^{2}q\hspace{0.17em}\text{poly}\left(\mathrm{log}T,\mathrm{log}n,\mathrm{log}1/\u03f5\right)/\u03f5$, where $T$ is the evolution time, $\u03f5$ 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 $\text{R}\ge \sqrt{2}$. Finally, we discuss potential applications, showing that the $\text{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 $\text{R}$.

### Sign up for PNAS alerts.

Get alerts for new articles, or get an alert when an article is cited.

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 (1–3) and PDEs (4–10). 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 $\text{poly}\left(\mathrm{log}n\right)$ 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\left(1/{\u03f5}^{T}\right)$. 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 $\text{poly}\left(\mathrm{log}n\right)$. 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 (20–22). 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 $\text{R}<1$, where the quantity $\text{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 $s{T}^{2}q\hspace{0.17em}\text{poly}\left(\mathrm{log}T,\mathrm{log}n,\mathrm{log}1/\u03f5\right)/\u03f5$, 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 $\u03f5$ is the allowed error (*Theorem 1*). In the regime $\text{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\left(1/{\u03f5}^{T}\right)$ (12) for sufficiently small $\u03f5$. 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 $\text{poly}\left(T\right)$]. Inhomogeneous equations arise in many applications, including, for example, the discretization of PDEs with nontrivial boundary conditions.We also show that the requirement $\text{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 $\text{R}\ge \sqrt{2}$, 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 $\mathrm{R}<1$, we find in numerical experiments that our proposed approach remains valid for larger $\mathrm{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 ODEHere $u={\left[{u}_{1},\dots ,{u}_{n}\right]}^{T}\in {\mathbb{R}}^{n}$, ${u}^{\otimes 2}={\left[{u}_{1}^{2},{u}_{1}{u}_{2},\dots ,{u}_{1}{u}_{n},{u}_{2}{u}_{1},\dots ,{u}_{n}{u}_{n-1},{u}_{n}^{2}\right]}^{T}\in {\mathbb{R}}^{{n}^{2}}$, each ${u}_{j}={u}_{j}\left(t\right)$ is a function of $t$ on the interval $\left[0,T\right]$ for $j\in \left[n\right]\u2254\left\{1,\dots ,n\right\}$, ${F}_{2}\in {\mathbb{R}}^{n\times {n}^{2}}$, ${F}_{1}\in {\mathbb{R}}^{n\times n}$ are time-independent matrices, and the inhomogeneity ${F}_{0}\left(t\right)\in {\mathbb{R}}^{n}$ is a ${C}^{1}$ continuous function of $t$.

$$\frac{\mathrm{d}u}{\mathrm{d}t}={F}_{2}{u}^{\otimes 2}+{F}_{1}u+{F}_{0}\left(t\right),\hspace{1em}\hspace{1em}u\left(0\right)={u}_{\mathrm{i}\mathrm{n}}.$$

[1]

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*${F}_{2}$, ${F}_{1}$,

*and*${F}_{0}$

*are*$s$-

*sparse*(

*i.e.*,

*have at most*$\mathrm{s}$

*nonzero entries in each row and column*), ${F}_{1}$

*is diagonalizable*,

*and the eigenvalues*${\lambda}_{j}$

*of*${F}_{1}$

*satisfy*$Re\left({\lambda}_{n}\right)\le \dots \le Re\left({\lambda}_{1}\right)<0$.

*We parametrize the problem in terms of the quantity*

$$\text{R}\u2254\frac{1}{|Re\left({\lambda}_{1}\right)|}\left(\Vert {u}_{\mathrm{i}\mathrm{n}}\Vert \Vert {F}_{2}\Vert +\frac{\Vert {F}_{0}\Vert}{\Vert {u}_{\mathrm{i}\mathrm{n}}\Vert}\right).$$

[2]

*For some given*$T>0$,

*we assume the values*$Re\left({\lambda}_{1}\right)$, $\Vert {F}_{2}\Vert $, $\Vert {F}_{1}\Vert $, $\Vert {F}_{0}\left(t\right)\Vert $

*for each*$t\in \left[0,T\right]$,

*and*$\Vert {F}_{0}\Vert \u2254\underset{t\in \left[0,T\right]}{max}\Vert {F}_{0}\left(t\right)\Vert $, $\Vert {F}_{0}^{\prime}\Vert \u2254\underset{t\in \left[0,T\right]}{max}\Vert {F}_{0}^{\prime}\left(t\right)\Vert $

*are known*,

*and that we are given oracles*${O}_{{F}_{2}}$, ${O}_{{F}_{1}}$,

*and*${O}_{{F}_{0}}$

*that provide the locations and values of the nonzero entries of*${F}_{2}$, ${F}_{1}$,

*and*${F}_{0}\left(t\right)$

*for any specified*$t$,

*respectively*,

*for any desired row or column*.

*We are also given the value*$\Vert {u}_{\mathrm{i}\mathrm{n}}\Vert $

*and an oracle*${O}_{x}$

*that maps*$\left|00\dots 0\right.\u27e9\in {\mathbb{C}}^{n}$

*to a quantum state proportional to*${u}_{\mathrm{i}\mathrm{n}}$.

*Our goal is to produce a quantum state proportional to*$u\left(T\right)$

*for some given*$T>0$

*within some prescribed*${\ell}^{2}$

*error tolerance*$\u03f5>0$.

When ${F}_{0}\left(t\right)=0$ (i.e., the ODE is homogeneous), the quantity $\text{R}=\frac{\Vert {u}_{\mathrm{i}\mathrm{n}}\Vert \Vert {F}_{2}\Vert}{|Re\left({\lambda}_{1}\right)|}$ 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 $\text{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 (20–22). 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 [with the tridiagonal block structurewhere ${\u0177}_{j}={u}^{\otimes j}\in {\mathbb{R}}^{{n}^{j}}$, ${\u0177}_{\mathrm{i}\mathrm{n}}=\left[{u}_{\mathrm{i}\mathrm{n}};{u}_{\mathrm{i}\mathrm{n}}^{\otimes 2};\dots ;{u}_{\mathrm{i}\mathrm{n}}^{\otimes N}\right]$, and ${A}_{j+1}^{j}\in {\mathbb{R}}^{{n}^{j}\times {n}^{j+1}}$, ${A}_{j}^{j}\in {\mathbb{R}}^{{n}^{j}\times {n}^{j}}$, ${A}_{j-1}^{j}\in {\mathbb{R}}^{{n}^{j}\times {n}^{j-1}}$ for $j\in \left[N\right]$ satisfyingTo construct a system of linear equations, we next divide $\left[0,T\right]$ into $m=T/h$ time steps and apply the forward Euler method on [where ${y}^{k}\in {\mathbb{R}}^{n+{n}^{2}+\cdots +{n}^{N}}$ approximates $\u0177\left(kh\right)$ for each $k\in {\left[m+1\right]}_{0}\u2254\left\{\mathrm{0,1},\dots ,m\right\}$, with ${y}^{0}={y}_{\mathrm{i}\mathrm{n}}\u2254\u0177\left(0\right)={\u0177}_{\mathrm{i}\mathrm{n}}$, and letting all ${y}^{k}$ be equal for $k\in {\left[m+p+1\right]}_{0}\backslash {\left[m+1\right]}_{0}$, for some sufficiently large integer $p$. (It is unclear whether another discretization could improve performance, as discussed further in

**1**], we apply the Carleman procedure to obtain the system of linear ODEs$$\frac{\mathrm{d}\u0177}{\mathrm{d}t}=A\left(t\right)\u0177+b\left(t\right),\hspace{1em}\hspace{1em}\u0177\left(0\right)={\u0177}_{\mathrm{i}\mathrm{n}}$$

[3]

$$\frac{\mathrm{d}}{\mathrm{d}t}\left(\begin{array}{c}\hfill {\u0177}_{1}\hfill \\ \hfill {\u0177}_{2}\hfill \\ \hfill \vdots \hfill \\ \hfill {\u0177}_{N}\hfill \end{array}\right)=\left(\begin{array}{ccc}\hfill {A}_{1}^{1}\hfill & \hfill {A}_{2}^{1}\hfill & \hfill \hfill \\ \hfill {A}_{1}^{2}\hfill & \hfill {A}_{2}^{2}\hfill & \hfill \ddots \hfill \\ \hfill \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill {A}_{N}^{N-1}\hfill \\ \hfill \hfill & \hfill \hfill & \hfill {A}_{N-1}^{N}\hfill & \hfill {A}_{N}^{N}\hfill \end{array}\right)\left(\begin{array}{c}\hfill {\u0177}_{1}\hfill \\ \hfill {\u0177}_{2}\hfill \\ \hfill \vdots \hfill \\ \hfill {\u0177}_{N}\hfill \end{array}\right)+\left(\begin{array}{c}\hfill {F}_{0}\left(t\right)\hfill \\ \hfill 0\hfill \\ \hfill \vdots \hfill \\ \hfill 0\hfill \end{array}\right),$$

[4]

$${A}_{j+1}^{j}={F}_{2}\otimes {I}^{\otimes j-1}\hspace{0.17em}+\hspace{0.17em}I\otimes {F}_{2}\otimes {I}^{\otimes j-2}\hspace{0.17em}+\hspace{0.17em}\dots \hspace{0.17em}+\hspace{0.17em}{I}^{\otimes j-1}\otimes {F}_{2},$$

[5]

$${A}_{j}^{j}={F}_{1}\otimes {I}^{\otimes j-1}\hspace{0.17em}+\hspace{0.17em}I\otimes {F}_{1}\otimes {I}^{\otimes j-2}\hspace{0.17em}+\hspace{0.17em}\dots \hspace{0.17em}+\hspace{0.17em}{I}^{\otimes j-1}\otimes {F}_{1},$$

[6]

$${A}_{j+1}^{j}={F}_{0}\left(t\right)\hspace{0.17em}\otimes \hspace{0.17em}{I}^{\otimes j-1}\hspace{0.17em}+\hspace{0.17em}I\hspace{0.17em}\otimes \hspace{0.17em}{F}_{0}\left(t\right)\hspace{0.17em}\otimes \hspace{0.17em}{I}^{\otimes j-2}\hspace{0.17em}+\hspace{0.17em}\dots \hspace{0.17em}+\hspace{0.17em}{I}^{\otimes j-1}\hspace{0.17em}\otimes \hspace{0.17em}{F}_{0}\left(t\right).$$

[7]

**3**], letting$${y}^{k+1}=\left[I+A\left(kh\right)h\right]{y}^{k}+b\left(kh\right)$$

[8]

*Discussion*.)In the above system, the first $n$ components of ${y}^{k}$ for $k\in {\left[m+p+1\right]}_{0}$ (i.e., ${y}_{1}^{k}$) is $\u03f5$-close to the exact solution $u\left(T\right)$, up to normalization. We apply the high-precision QLSA (23) to [

**8**] and postselect on $k$ to produce ${y}_{1}^{k}/\Vert {y}_{1}^{k}\Vert $ for some $k\in {\left[m+p+1\right]}_{0}\backslash {\left[m\right]}_{0}$.*SI Appendix*, Section 2, gives more details about the system of linear equations and the measurement procedure.We control the additive approximation error $\u03f5$ 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*$\text{R}<1$.

*Let*$q\u2254\Vert {u}_{\mathrm{i}\mathrm{n}}\Vert /\Vert u\left(T\right)\Vert $.

*There exists a quantum algorithm producing a state that approximates*$u\left(T\right)/\Vert u\left(T\right)\Vert $

*with error at most*$\u03f5\le 1$,

*succeeding with probability*$\mathrm{\Omega}\left(1\right)$,

*with query and gate complexity*

$$\frac{s{T}^{2}q}{\u03f5}\cdot \text{poly}\left(\mathrm{log}T,\mathrm{log}n,\mathrm{log}1/\u03f5\right).$$

[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 $\text{R}\ge \sqrt{2}$. 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 $\text{R}\ge \sqrt{2}$, this procedure distinguishes states with overlap $1-\u03f5$ in time $\text{poly}\left(\mathrm{log}\left(1/\u03f5\right)\right)$. 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*$\text{R}\ge \sqrt{2}$.

*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\left(T\right)/\Vert u\left(T\right)\Vert $

*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 ${F}_{1}$ 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*${\lambda}_{j}$

*of*${F}_{1}$

*satisfy*$Re\left({\lambda}_{n}\right)\le \dots \le Re\left({\lambda}_{1}\right)<0$.

*Then for any*$j\in \left[N\right]$,

*the error*${\eta}_{j}\left(t\right)\u2254{u}^{\otimes j}\left(t\right)-{\u0177}_{j}\left(t\right)$

*satisfies*

$$\Vert {\eta}_{j}\left(t\right)\Vert \le tN\Vert {F}_{2}\Vert {\Vert {u}_{\mathrm{i}\mathrm{n}}\Vert}^{N+1}$$

[10]

*where we require*$\text{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*$\text{R}<1$

*as defined in*[

**2**].

*Choose a sufficiently small time step*$h$.

*Suppose the error from Carleman linearization satisfies*$\Vert \eta \left(t\right)\Vert \le \Vert u\left(T\right)\Vert /4$.

*Then the global error from the forward Euler method*[

**8**]

*on the interval*$\left[0,T\right]$

*for*[

**3**]

*satisfies*

$$\Vert {\u0177}_{1}\left(T\right)-{y}_{1}\left(T\right)\Vert \le 3{N}^{2.5}Th\left[{\left(\Vert {F}_{2}\Vert +\Vert {F}_{1}\Vert +\Vert {F}_{0}\Vert \right)}^{2}+\Vert {F}_{0}^{\prime}\Vert \right].$$

[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*

$$\kappa \le 3\left(m+p+1\right).$$

[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*$\Vert {u}_{\mathrm{i}\mathrm{n}}\Vert $,

*and let*${O}_{x}$

*be an oracle that maps*$\left|00\dots 0\right.\u27e9\in {\mathbb{C}}^{n}$

*to a normalized quantum state*$\left|{u}_{\mathrm{i}\mathrm{n}}\right.\u27e9$

*proportional to*${u}_{\mathrm{i}\mathrm{n}}$.

*Assume we are also given the values*$\Vert {F}_{0}\left(t\right)\Vert $

*for each*$t\in \left[0,T\right]$,

*and let*${O}_{{F}_{0}}$

*be an oracle that provides the locations and values of the nonzero entries of*${F}_{0}\left(t\right)$

*for any specified*$t$.

*Then a quantum state encoding the right-hand side of the linear system can be prepared using*$O\left(N\right)$

*queries to*${O}_{x}$

*and*$O\left(m\right)$

*queries to*${O}_{{F}_{0}}$,

*with gate complexity larger by a*$\text{poly}\left(\mathrm{log}N,\mathrm{log}n\right)$

*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*$\Vert \eta \left(t\right)\Vert \le \Vert u\left(T\right)\Vert /4$,

*and the error from Euler method satisfies*$\Vert \u0177\left(T\right)-{y}^{m}\Vert \le \Vert u\left(T\right)\Vert /4$.

*Then the probability of measuring a state*$\left|{y}_{1}^{k}\right.\u27e9$

*for*$k={\left[m+p+1\right]}_{0}\backslash {\left[m+1\right]}_{0}$

*satisfies*

$${P}_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{e}}\ge \frac{p+1}{9\left(m+p+1\right)N{q}^{2}}.$$

[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 $\text{R}$ defined in [

**2**]. Whereas there is an efficient quantum algorithm for $\text{R}<1$ (as shown in*Theorem 1*), we show here that the problem is intractable for $\text{R}\ge \sqrt{2}$. 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*$|\psi \rangle ,|\varphi \rangle $

*be states of a qubit with*$\left|\u27e8\psi |\varphi \u27e9\right|=1-\u03f5$.

*Suppose we are given either a black box that prepares*$\left|\psi \right.\u27e9$

*or a black box that prepares*$\left|\varphi \right.\u27e9$.

*Then any bounded-error protocol for determining whether the state is*$\left|\psi \right.\u27e9$

*or*$\left|\varphi \right.\u27e9$

*must use*$\mathrm{\Omega}\left(1/\u03f5\right)$

*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*$\text{R}\ge \sqrt{2}$,

*and two states of a qubit with overlap*$1-\u03f5$ (

*for*$0<\u03f5<1-3/\sqrt{10}$)

*as possible initial conditions*,

*such that the two final states after evolution time*$T=O\left(\mathrm{log}\left(1/\u03f5\right)\right)$

*have an overlap no larger than*$3/\sqrt{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 modeldescribes the population of $P$ individuals composed of four components: susceptible (${P}_{S}$), exposed (${P}_{E}$), infected (${P}_{I}$), and recovered (${P}_{R}$). Here ${r}_{\text{tra}}$ is the rate of transmission from an infected to a susceptible person, ${r}_{\text{vac}}$ is the vaccination rate, ${T}_{\text{lat}}$ is the latent time until an exposed person becomes infectious, ${T}_{\text{inf}}$ is the infectious time that an infectious person can infect others, and $\mathrm{\Lambda}$ 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 $\text{R}<1$ (as discussed in more detail in

$$\frac{\mathrm{d}{P}_{S}}{\mathrm{d}t}=-\mathrm{\Lambda}\frac{{P}_{S}}{P}-{r}_{\text{vac}}{P}_{S}+\mathrm{\Lambda}-{r}_{\text{tra}}{P}_{S}\frac{{P}_{I}}{P}$$

[14]

$$\frac{\mathrm{d}{P}_{E}}{\mathrm{d}t}=-\mathrm{\Lambda}\frac{{P}_{E}}{P}-\frac{{P}_{E}}{{T}_{\text{lat}}}+{r}_{\text{tra}}{P}_{S}\frac{{P}_{I}}{P}$$

[15]

$$\frac{\mathrm{d}{P}_{I}}{\mathrm{d}t}=-\mathrm{\Lambda}\frac{{P}_{I}}{P}+\frac{{P}_{E}}{{T}_{\text{lat}}}-\frac{{P}_{I}}{{T}_{\text{inf}}}$$

[16]

$$\frac{\mathrm{d}{P}_{R}}{\mathrm{d}t}=-\mathrm{\Lambda}\frac{{P}_{R}}{P}+{r}_{\text{vac}}{P}_{S}+\frac{{P}_{I}}{{T}_{\text{inf}}}$$

[17]

*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 equationwith linear damping and a forcing term. Equations of the form [

$${\partial}_{t}\mathbf{u}+\left(\mathbf{u}\cdot \nabla \right)\mathbf{u}+\beta \mathbf{u}=\nu {\nabla}^{2}\mathbf{u}+\mathbf{f}$$

[18]

**18**] are ubiquitous in fluid mechanics (27), and related models such as those studied in refs. 28–30 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, $\nu =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 equationwhich is the one-dimensional case of Eq.

$${\partial}_{t}u+u{\partial}_{x}u=\nu {\partial}_{x}^{2}u+f,$$

[19]

**18**with $\beta =0$. This is often used as a simple model of convective flow (35). For concreteness, let the initial condition be $u\left(x,0\right)={U}_{0}\mathrm{sin}\left(2\pi x/{L}_{0}\right)$ on the domain $x\in \left[-{L}_{0}/2,{L}_{0}/2\right]$, and use Dirichlet boundary conditions $u\left(-{L}_{0}/\mathrm{2,0}\right)=u\left({L}_{0}/\mathrm{2,0}\right)=0$. We force this equation using a localized off-center Gaussian with a sinusoidal time dependence, given by $f\left(x,t\right)={U}_{0}\mathrm{exp}\left(-\frac{{\left(x-{L}_{0}/4\right)}^{2}}{2{\left({L}_{0}/32\right)}^{2}}\right)\mathrm{cos}\left(2\pi t\right)$. To solve this equation using the Carleman method, we discretize the spatial domain into ${n}_{x}$ 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 $\text{R}\approx 44$, which violates the requirement $\text{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 $\text{R}>1$.Fig. 1.

## 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 $\text{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 $\text{R}<1$. Furthermore, we established a lower bound showing that for general quadratic differential equations with $\text{R}\ge \sqrt{2}$, 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 $\text{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 $\mathrm{\Theta}\left({n}^{N}\right)$. 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 $\Vert {\lambda}_{1}\Vert T\ll 1$ but let the nonlinearity and forcing be even smaller such that $\text{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 $\text{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 $1\le \text{R}<\sqrt{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 $\text{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, ${T}^{2}\text{poly}\left(\mathrm{log}T\right)$. 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 $\mathrm{B}\mathrm{Q}\mathrm{P}=\mathrm{P}\mathrm{P}$ (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 $\text{R}<1$ is not a strict requirement for the viscous Burgers equation, since we observe convergence even though $\text{R}\approx 44$. 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 $\text{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)

- Download
- 384.02 KB

## 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

#### Classifications

#### Copyright

© 2021. Published under the PNAS license.

#### 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

#### 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

#### Competing Interests

The authors declare no competing interest.

## Metrics & Citations

### Metrics

#### 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.