Efficient quantum algorithm for
dissipative nonlinear differential equations
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 -dimensional ordinary differential equations. Assuming , where is a parameter characterizing the ratio of the nonlinearity and forcing to the linear dissipation, this algorithm has complexity , where is the evolution time, is the allowed error, and measures decay of the solution. This is an exponential improvement over the best previous quantum algorithms, whose complexity is exponential in . 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 novel 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 . Finally, we discuss potential applications, showing that the 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 .
1 Introduction
Models governed by both ordinary differential equations (ODEs) and partial differential equations (PDEs) arise extensively in natural and social science, 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. Note that polynomials of degree higher than two, and even more general nonlinearities, can be reduced to the quadratic case by introducing additional variables [35, 30]. The quadratic case also directly includes many archetypal models, such as the logistic equation in biology, the Lorenz system in atmospheric dynamics, and the Navier–Stokes equations in fluid dynamics.
Quantum algorithms offer the prospect of rapidly characterizing the solutions of high-dimensional systems of linear ODEs [7, 9, 19] and PDEs [22, 15, 43, 23, 20, 28, 40]. Such algorithms can produce a quantum state proportional to the solution of a sparse (or block-encoded) -dimensional system of linear differential equations in time using the quantum linear system algorithm [32].
Early work on quantum algorithms for differential equations already considered the nonlinear case [39]. It 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 due to exponentially increasing resources used to maintain sufficiently many copies of the solution to represent the nonlinearity throughout the evolution.
Recently, heuristic quantum algorithms for nonlinear ODEs have been studied. Reference [34] explores a linearization technique known as the Koopman–von Neumann method that might be amenable to the quantum linear system algorithm. In [27], the authors provide a high-level description of how linearization can help solve nonlinear equations on a quantum computer. However, neither paper makes precise statements about concrete implementations or running times of quantum algorithms. The recent preprint [41] 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 [21] could be avoided in their approach.
While quantum mechanics is described by linear dynamics, possible nonlinear modifications of the theory have been widely studied. Generically, such modifications enable quickly solving hard computational problems (e.g., solving unstructured search among items in time ), making nonlinear dynamics exponentially difficult to simulate in general [2, 1, 21]. Therefore, constructing efficient quantum algorithms for general classes of nonlinear dynamics has been considered largely out of reach.
In this article, we design and analyze a quantum algorithm that overcomes this limitation using Carleman linearization [16, 37, 30]. This approach embeds polynomial nonlinearities into an infinite-dimensional system of linear ODEs, and then truncates it to obtain a finite-dimensional linear approximation. The Carleman method has previously been used in the analysis of dynamical systems [49, 3, 42] and the design of control systems [13, 46, 31], but to the best of our knowledge it has not been employed in the context of quantum algorithms. We discretize the finite ODE system in time using the forward Euler method and solve the resulting linear equations with the quantum linear system algorithm [32, 18]. We control the approximation error of this approach by combining a novel convergence theorem with a bound for the global error of the Euler method. Furthermore, we provide an upper bound for the condition number of the linear system and lower bound the success probability of the final measurement. Subject to the condition , where the quantity (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 algorithm is , where is the sparsity, is the evolution time, quantifies the decay of the final solution relative to the initial condition, is the dimension, and is the allowed error (see Theorem 1). In the regime , this is an exponential improvement over [39], which has complexity exponential in .
Note that the solution cannot decay exponentially in for the algorithm to be efficient, as captured by the dependence of the complexity on —a known limitation of quantum ODE algorithms [9]. For homogeneous ODEs with , the solution necessarily decays exponentially in time (see equation (4.11)), so the algorithm is not asymptotically efficient. However, even for solutions with exponential decay, we still find an improvement over the best previous result [39] 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 2), or can decay slowly (i.e., can be ). Inhomogeneous equations arise in many applications, including for example the discretization of PDEs with nontrivial boundary conditions.
We also provide a quantum lower bound for the worst-case complexity of simulating strongly nonlinear dynamics, showing that the algorithm’s condition cannot be significantly improved in general (Theorem 2). Following the approach of [2, 21], we construct a protocol for distinguishing two states of a qubit driven by a certain quadratic ODE. Provided , this procedure distinguishes states with overlap in time . Since nonorthogonal quantum states are hard to distinguish, this implies a lower bound on the complexity of the quantum ODE problem.
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 provided the Reynolds number is sufficiently small. We also note that while the formal validity of our arguments assumes , we find in one numerical experiment that our proposed approach remains valid for larger (see Section 6).
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 Section 7.
The remainder of this paper is structured as follows. Section 2 introduces the quantum quadratic ODE problem. Section 3 presents the Carleman linearization procedure and describes its performance. Section 4 gives a detailed analysis of the quantum Carleman linearization algorithm. Section 5 establishes a quantum lower bound for simulating quadratic ODEs. Section 6 describes how our approach could be applied to several well-known ODEs and PDEs and presents numerical results for the case of the viscous Burgers equation. Finally, we conclude with a discussion of the results and some possible future directions in Section 7.
2 Quadratic ODEs
We focus on an initial value problem described by the -dimensional quadratic ODE
| (2.1) |
Here , , each is a function of on the interval for , , are time-independent matrices, and the inhomogeneity is a continuous function of . We let denote the spectral norm.
The main computational problem we consider is as follows.
Problem 1.
In the quantum quadratic ODE problem, we consider an -dimensional quadratic ODE as in (2.1). We assume , , and are -sparse (i.e., have at most nonzero entries in each row and column), is diagonalizable, and that the eigenvalues of satisfy . We parametrize the problem in terms of the quantity
| (2.2) |
For some given , we assume the values , , , for each , and , are known, and that we are given oracles , , and that provide the locations and values of the nonzero entries of , , and for any specified , respectively, for any desired row or column. We are also given the value and an oracle that maps to a quantum state proportional to . Our goal is to produce a quantum state proportional to for some given within some prescribed error tolerance .
When (i.e., the ODE is homogeneous), the quantity is qualitatively similar to the Reynolds number, which characterizes the ratio of the (nonlinear) convective forces to the (linear) viscous forces within a fluid [14, 38]. More generally, quantifies the combined strength of the nonlinearity and the inhomogeneity relative to dissipation.
Note that without loss of generality, given a quadratic ODE satisfying (2.1) with , we can modify it by rescaling with a suitable constant to satisfy
| (2.3) |
and
| (2.4) |
with left unchanged by the rescaling. We use this rescaling in our algorithm and its analysis. With this rescaling, a small implies both small and small relative to .
3 Quantum Carleman linearization
As discussed in Section 1, it is challenging to directly simulate quadratic ODEs using quantum computers, and indeed the complexity of the best known quantum algorithm is exponential in the evolution time [39]. However, for a dissipative nonlinear ODE without a source, any quadratic nonlinear effect will only be significant for a finite time because of the dissipation. To exploit this, we can create a linear system that approximates the initial nonlinear evolution within some controllable error. After the nonlinear effects are no longer important, the linear system properly captures the almost-linear evolution from then on.
We develop such a quantum algorithm using the concept of Carleman linearization [16, 37, 30]. 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 equations, where the truncation level depends on the allowed approximation error, giving a finite linear ODE system.
Let us describe the Carleman linearization procedure in more detail. Given a system of quadratic ODEs (2.1), we apply the Carleman procedure to obtain the system of linear ODEs
| (3.1) |
with the tri-diagonal block structure
| (3.2) |
where , , and , , for satisfying
| (3.3) | ||||
| (3.4) | ||||
| (3.5) |
Note that is a -sparse matrix. The dimension of (3.1) is
| (3.6) |
To construct a system of linear equations, we next divide into time steps and apply the forward Euler method on (3.1), letting
| (3.7) |
where approximates for each , with , and letting all be equal for , for some sufficiently large integer . (It is unclear whether another discretization could improve performance, as discussed further in Section 7.) This gives an linear system
| (3.8) |
that encodes (3.7) and uses it to produce a numerical solution at time , where
| (3.9) |
and
| (3.10) |
with a normalizing factor . Observe that the system (3.8) has the lower triangular structure
| (3.11) |
In the above system, the first components of for (i.e., ) approximate the exact solution , up to normalization. We apply the high-precision quantum linear system algorithm (QLSA) [18] to (3.8) and postselect on to produce for some . The resulting error is at most
| (3.12) |
This error includes contributions from both Carleman linearization and the forward Euler method. (The QLSA also introduces error, which we bound separately. Note that we could instead apply the original QLSA [32] instead of its subsequent improvement [18], but this would slightly complicate the error analysis and might perform worse in practice.)
Now we state our main algorithmic result.
Theorem 1 (Quantum Carleman linearization algorithm).
Consider an instance of the quantum quadratic ODE problem as defined in Problem 1, with its Carleman linearization as defined in (3.1). Assume . Let
| (3.13) |
There exists a quantum algorithm producing a state that approximates with error at most , succeeding with probability , with a flag indicating success, using at most
| (3.14) |
queries to the oracles , , and . The gate complexity is larger than the query complexity by a factor of . Furthermore, if the eigenvalues of are all real, the query complexity is
| (3.15) |
and the gate complexity is larger by a factor of .
4 Algorithm analysis
In this section we establish several lemmas and use them to prove Theorem 1.
4.1 Solution error
The solution error has three contributions: the error from applying Carleman linearization to (2.1), the error in the time discretization of (3.1) by the forward Euler method, and the error from the QLSA. Since the QLSA produces a solution with error at most with complexity [18], we focus on bounding the first two contributions.
4.1.1 Error from Carleman linearization
First, we provide an upper bound for the error from Carleman linearization for arbitrary evolution time . To the best of our knowledge, the first and only explicit bound on the error of Carleman linearization appears in [30]. However, they only consider homogeneous quadratic ODEs; and furthermore, to bound the error for arbitrary , they assume the logarithmic norm of is negative (see Theorems 4.2 and 4.3 of [30]), which is too strong for our case. Instead, we give a novel analysis under milder conditions, providing the first convergence guarantee for general inhomogeneous quadratic ODEs.
We begin with a lemma that describes the decay of the solution of (2.1).
Lemma 1.
Proof.
Consider the derivative of . We have
| (4.2) |
If , then
| (4.3) |
Letting , , and with , we consider a -dimensional quadratic ODE
| (4.4) |
Since , the discriminant satisfies
| (4.5) |
Thus, defined in (4.1) are distinct real roots of . Since and , we have . We can rewrite the ODE as
| (4.6) |
Letting , we obtain an associated homogeneous quadratic ODE
| (4.7) |
Since the homogeneous equation has the closed-form solution
| (4.8) |
the solution of the inhomogeneous equation can be obtained as
| (4.9) |
Therefore we have
| (4.10) |
Since , is located between the two roots and , and thus . This implies in (4.10) decreases from , so we have for any . ∎
We remark that since and , so approaches to a stationary point of the right-hand side of (2.1) (called an attractor in the theory of dynamical systems).
Note that for a homogeneous equation (i.e., ), this shows that the dissipation inevitably leads to exponential decay. In this case we have , so (4.10) gives
| (4.11) |
which decays exponentially in .
On the other hand, as mentioned in the introduction, the solution of a dissipative inhomogeneous quadratic ODE can remain asymptotically nonzero. Here we present an example of this. Consider a time-independent uncoupled system with , , with , , , , and . We see that each decreases from to , with . Hence, the norm of is bounded as for any . In general, it is hard to lower bound , but the above example shows that a nonzero inhomogeneity can prevent the norm of the solution from decreasing to zero.
We now give an upper bound on the error of Carleman linearization.
Lemma 2.
Proof.
The exact solution of the original quadratic ODE (2.1) satisfies
| (4.13) |
and the approximated solution satisfies (3.2). Comparing these equations, we have
| (4.14) |
with the tri-diagonal block structure
| (4.15) |
Consider the derivative of . We have
| (4.16) |
For , we bound each term as
| (4.17) | ||||
Define a matrix with nonzero entries , , and . Then
| (4.18) |
Since , is strictly diagonally dominant and thus the eigenvalues of satisfy . Thus we have
| (4.19) |
For , we have
| (4.20) |
Since , and , we have
| (4.21) |
Note that (4.25) can be bounded alternatively by
| (4.27) |
and thus . We select (4.25) because it avoids including an additional parameter .
We also give an improved analysis that works for homogeneous quadratic ODEs () under milder conditions. This analysis follows the proof in [30] closely.
Corollary 1.
Proof.
We again consider satisfying (4.14). Since , (4.14) reduces to a time-independent ODE with an upper triangular block structure,
| (4.30) |
and
| (4.31) |
We proceed by backward substitution. Since , we have
| (4.32) |
For , (3.4) gives and (3.3) gives . By Lemma 1, . We can therefore upper bound (4.32) by
| (4.33) | ||||
For , (4.30) gives
| (4.34) |
Again, since , we have
| (4.35) |
which has the upper bound
| (4.36) | ||||
where we used (4.33) in the last step. Iterating this procedure for , we find
| (4.37) | ||||
While Problem 1 makes some strong assumptions about the system of differential equations, they appear to be necessary for our analysis. Specifically, the conditions and are required to ensure arbitrary-time convergence.
Since the Euler method for (3.1) is unstable if [7, 25], we only consider the case . If , (4.40) reduces to
| (4.42) |
giving the error bound
| (4.43) |
instead of (4.41). Then the error bound can be made arbitrarily small for a finite time by increasing , but after , the error bound diverges.
Furthermore, if , Bernoulli’s inequality gives
| (4.44) |
where the right-hand side upper bounds as in (4.41). Assuming is smaller than , we require
| (4.45) |
In other words, to apply (4.29) for the Carleman linearization procedure, the truncation order given by Lemma 2 must grow exponentially with .
4.1.2 Error from forward Euler method
Next, we provide an upper bound for the error incurred by approximating (3.1) with the forward Euler method. This problem has been well studied for general ODEs. Given an ODE on with an inhomogenity that is an -Lipschitz continuous function of , the global error of the solution is upper bounded by , although in most cases this bound overestimates the actual error [4]. To remove the exponential dependence on in our case, we derive a tighter bound for time discretization of (3.1) in Lemma 3 below. This lemma is potentially useful for other ODEs as well and can be straightforwardly adapted to other problems.
Lemma 3.
Consider an instance of the quantum quadratic ODE problem as defined in Problem 1, with as defined in (2.2). Choose a time step
| (4.46) |
in general, or
| (4.47) |
if the eigenvalues of are all real. Suppose the error from Carleman linearization as defined in Lemma 2 is bounded by
| (4.48) |
where is defined in (3.13). Then the global error from the forward Euler method (3.7) on the interval for (3.1) satisfies
| (4.49) |
Proof.
We define a linear system that locally approximates the initial value problem (3.1) on for as
| (4.50) |
where is the exact solution of (3.1). For , we denote the local truncation error by
| (4.51) |
and the global error by
| (4.52) |
where in (3.7) is the numerical solution. Note that .
For the local truncation error, we Taylor expand with a second-order Lagrange remainder, giving
| (4.53) |
for some . Since by (3.1), we have
| (4.54) |
and thus the local error (4.51) can be bounded as
| (4.55) |
where .
By the triangle inequality, the global error (4.52) can therefore be bounded as
| (4.56) |
Since and are obtained by the same linear system with different right-hand sides, we have the upper bound
| (4.57) | ||||
In order to provide an upper bound for for all , we write
| (4.58) |
where
| (4.59) | ||||
| (4.60) | ||||
| (4.61) |
We provide upper bounds separately for , , and , and use the bound .
The eigenvalues of consist of all -term sums of the eigenvalues of . More precisely, they are where are the eigenvalues of and is a -tuple of indices. The eigenvalues of are thus . With , we have
| (4.62) | ||||
where we used by (4.46). Therefore
| (4.63) |
We also have
| (4.64) |
and
| (4.65) |
Using the bounds (4.63) and (4.64), we aim to select the value of to ensure
| (4.66) |
The assumption (4.46) implies
| (4.67) |
(note that the denominator is non-zero due to (2.3)). Then we have
| (4.68) | ||||||
so as claimed.
The choice (4.46) can be improved if an upper bound on is known. In particular, if , (4.67) simplifies to
| (4.69) |
which is satisfied by (4.47) using .
Using this in (4.57), we have
| (4.70) |
Plugging (4.70) into (4.56) iteratively, we find
| (4.71) |
Using (4.55), this shows that the global error from the forward Euler method is bounded by
| (4.72) |
and when , ,
| (4.73) |
Finally, we remove the dependence on . Since
| (4.74) | ||||
we have
| (4.75) |
Maximizing each term for , we have
| (4.76) |
| (4.77) |
| (4.78) |
| (4.79) |
and using , by (4.48), and , we have
| (4.80) | ||||
for all . Therefore, substituting the bounds (4.76)–(4.80) into (4.75), we find
| (4.81) | ||||
Thus, (4.73) gives
| (4.82) |
and when , ,
| (4.83) |
as claimed. ∎
4.2 Condition number
Now we upper bound the condition number of the linear system.
Lemma 4.
Proof.
We begin by upper bounding . We write
| (4.85) |
where
| (4.86) | ||||
| (4.87) | ||||
| (4.88) |
Clearly . Furthermore, by (4.66), which follows from the choice of time step (4.46). Therefore,
| (4.89) |
Next we upper bound
| (4.90) |
We express as
| (4.91) |
where satisfies
| (4.92) |
Given any for , we define
| (4.93) |
where . We first upper bound , and then use this to upper bound .
We consider two cases. First, for fixed , we directly calculate for each by the forward Euler method (3.7), giving
| (4.94) |
| (4.95) | ||||
4.3 State preparation
We now describe a procedure for preparing the right-hand side of the linear system (3.8), whose entries are composed of and for .
The initial vector is a direct sum over spaces of different dimensions, which is cumbersome to prepare. Instead, we prepare an equivalent state that has a convenient tensor product structure. Specifically, we embed into a slightly larger space and prepare the normalized version of
| (4.102) |
where is some standard vector (for simplicity, we take ). If lives in a vector space of dimension , then lives in a space of dimension while lives in a slightly smaller space of dimension . Using standard techniques, all the operations we would otherwise apply to can be applied instead to , with the same effect.
Lemma 5.
Assume we are given the value , and let be an oracle that maps to a normalized quantum state proportional to . Assume we are also given the values for each , and let be an oracle that provides the locations and values of the nonzero entries of for any specified . Then the quantum state defined in (3.10) (with replaced by ) can be prepared using queries to and queries to , with gate complexity larger by a factor.
Proof.
We first show how to prepare the state
| (4.103) |
where
| (4.104) |
This state can be prepared using queries to the initial state oracle applied in superposition to the intermediate state
| (4.105) |
To efficiently prepare , notice that
| (4.106) |
where (assuming for simplicity that is a power of ) and is the -bit binary expansion of . Observe that
| (4.107) |
where
| (4.108) |
(Notice that .) Each tensor factor in (4.107) is a qubit that can be produced in constant time. Overall, we prepare these qubit states and then apply times.
We now discuss how to prepare the state
| (4.109) |
in which we replace by in (3.10), and define
| (4.110) |
This state can be prepared using the above procedure for and queries to with that implement for , applied in superposition to the intermediate state
| (4.111) |
Here the queries are applied conditionally upon the value in the first register: we prepare if the first register is and if the first register is for . We can prepare (i.e., perform a unitary transform mapping ) in time complexity [48] using the known values of and .
Overall, we use queries to and queries to to prepare . The gate complexity is larger by a factor. ∎
4.4 Measurement success probability
After applying the QLSA to (3.8), we perform a measurement to extract a final state of the desired form. We now consider the probability of this measurement succeeding.
Lemma 6.
Consider an instance of the quantum quadratic ODE problem defined in Problem 1, with the QLSA applied to the linear system (3.8) using the forward Euler method (3.7) with time step (4.46). Suppose the error from Carleman linearization satisfies as in (4.48), and the global error from the forward Euler method as defined in Lemma 3 is bounded by
| (4.112) |
where is defined in (3.13). Then the probability of measuring a state for satisfies
| (4.113) |
where is also defined in (3.13).
Proof.
The idealized quantum state produced by the QLSA applied to (3.8) has the form
| (4.114) |
where the states and for and are subnormalized to ensure .
We decompose the state as
| (4.115) |
where
| (4.116) | ||||
Note that for all . We lower bound
| (4.117) |
by lower bounding the terms of the product
| (4.118) |
First, according to (4.48) and (4.112), the exact solution and the approximate solution defined in (3.8) satisfy
| (4.119) |
Since , using (4.119), we have
| (4.120) |
Choosing , we have . Using amplitude amplification, iterations suffice to succeed with constant probability.
4.5 Proof of Theorem 1
Proof.
We first present the quantum Carleman linearization (QCL) algorithm and then analyze its complexity.
The QCL algorithm.
We start by rescaling the system to satisfy (2.3) and (2.4). Given a quadratic ODE (2.1) satisfying (where is defined in (2.2)), we define a scaling factor , and rescale to
| (4.125) |
Replacing by in (2.1), we have
| (4.126) | ||||
We let , , and so that
| (4.127) | ||||
Note that is invariant under this rescaling, so still holds for the rescaled equation.
Concretely, we take111In fact, one can show that any suffices to satisfy (2.3) and (2.4).
| (4.128) |
After rescaling, the new quadratic ODE satisfies . Since by Lemma 1, we have , so (2.4) holds. Furthermore, is located between the two roots and , which implies as shown in Lemma 1, so (2.3) holds for the rescaled problem.
Having performed this rescaling, we henceforth assume that (2.3) and (2.4) are satisfied. We then introduce the choice of parameters as follows. Given and an error bound , we define
| (4.129) |
Given , , and , we choose
| (4.130) |
Since by (4.129) and , Lemma 2 gives
| (4.131) |
Thus, (4.48) holds since .
Now we discuss the choice of . On the one hand, must satisfy (4.46) to satisfy the conditions of Lemma 3 and Lemma 4. On the other hand, Lemma 3 gives the upper bound
| (4.132) |
This also ensures that (4.112) holds since . Thus, we choose
| (4.133) | ||||
Combining (4.131) with (4.132), we have
| (4.134) |
Thus, (4.119) holds since . Using
| (4.135) |
and (4.134), we obtain
| (4.136) |
i.e., the normalized output state is -close to .
We follow the procedure in Lemma 5 to prepare the initial state . We apply the QLSA [18] to the linear system (3.8) with , giving a solution . We then perform a measurement to obtain a normalized state of for some and . By Lemma 6, the probability of obtaining a state for some , giving the normalized vector , is
| (4.137) |
By amplitude amplification, we can achieve success probability with repetitions of the above procedure.
Analysis of the complexity.
By Lemma 5, the right-hand side in (3.8) can be prepared with queries to and queries to , with gate complexity larger by a factor. The matrix in (3.8) is an matrix with nonzero entries in any row or column. By Lemma 4 and our choice of parameters, the condition number of is at most
| (4.138) | ||||
Here we use and
| (4.139) | ||||
The first inequality above is from the sum of and , where . Consequently, by Theorem 5 of [18], the QLSA produces the state with
| (4.140) | ||||
queries to the oracles , , , and . Using steps of amplitude amplification to achieve success probability , the overall query complexity of our algorithm is
| (4.141) | ||||
and the gate complexity exceeds this by a factor of .
If the eigenvalues of are all real, by (4.47), the condition number of is at most
| (4.142) | ||||
Similarly, the QLSA produces the state with
| (4.143) |
queries to the oracles , , , and . Using amplitude amplification to achieve success probability , the overall query complexity of the algorithm is
| (4.144) |
and the gate complexity is larger by a factor of as claimed. ∎
5 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 defined in (2.2). Whereas there is an efficient quantum algorithm for (as shown in Theorem 1), we show here that the problem is intractable for .
Theorem 2.
Assume . 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 with bounded error must have worst-case time complexity exponential in .
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. Note that since our algorithm only approximates the quantum state corresponding to the solution, we must lower bound the query complexity of approximating the solution of a quadratic ODE.
5.1 Hardness of state discrimination
Previous work on the computational power of nonlinear quantum mechanics shows that the ability to distinguish non-orthogonal states can be applied to solve unstructured search (and other hard computational problems) [2, 1, 21]. Here we show a similar limitation using an information-theoretic argument.
Lemma 7.
Let be states of a qubit with . Suppose we are either given 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 queries.
Proof.
Using the black box times, we prepare states with overlap . By the well-known relationship between fidelity and trace distance, these states have trace distance at most . Therefore, by the Helstrom bound (which states that the advantage over random guessing for the best measurement to distinguish two quantum states is given by their trace distance [33]), we need to distinguish the states with bounded error. ∎
5.2 State discrimination with nonlinear dynamics
Lemma 7 can be used to establish limitations on the ability of quantum computers to simulate 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 [2, 1] and the Gross-Pitaevskii equation [21]), here we aim to show the difficulty of efficiently simulating more general nonlinear ODEs—in particular, quadratic ODEs with dissipation—using quantum algorithms.
Lemma 8.
There exists an instance of the quantum quadratic ODE problem as defined in Problem 1 with , and two states of a qubit with overlap (for ) as possible initial conditions, such that the two final states after evolution time have an overlap no larger than .
Proof.
Consider a -dimensional system of the form
| (5.1) | ||||
for some , with an initial condition satisfying . According to the definition of in (2.2), we have , so henceforth we write this parameter as . The analytic solution of (5.1) is
| (5.2) | ||||
When , is finite within the domain
| (5.3) |
when , we have for all ; and when , goes to as . The behavior of depends similarly on .
Without loss of generality, we assume . For , both and are finite within the domain (5.3), which we consider as the domain of .
Now we consider 1-qubit states that provide inputs to (5.1). Given a sufficiently small , we first define by
| (5.4) |
We then construct two 1-qubit states with overlap , namely
| (5.5) |
and
| (5.6) |
Then the overlap between the two initial states is
| (5.7) |
The initial overlap (5.7) is larger than the target overlap in Lemma 8 provided . For simplicity, we denote
| (5.8) | |||
and let and denote solutions of (5.1) with initial conditions and , respectively. Since , we see that increases with , satisfying
| (5.9) |
and
| (5.10) |
for any time , whatever the behavior of .
We now study the outputs of our problem. For the state , the initial condition for (5.1) is . Thus, the output for any is
| (5.11) |
For the state , the initial condition for (5.1) is . We now discuss how to select a terminal time to give a useful output state . For simplicity, we denote the ratio of and by
| (5.12) |
Noticing that goes to infinity as approaches , while remains finite within (5.3), there exists a terminal time such that222More concretely, we take that upper bounds on the domain , in which is a finite value since . Then there exists a terminal time such that , and hence .
| (5.13) |
The normalized output state at this time is
| (5.14) |
Finally, we estimate the evolution time , which is implicitly defined by (5.13). We can upper bound its value by . According to (5.3), we have
| (5.16) |
since the function decreases monotonically with for . Using (5.7) to rewrite this expression in terms of , we have
| (5.17) |
which scales like as . Therefore as claimed. ∎
5.3 Proof of Theorem 2
We now establish our main lower bound result.
Proof.
As introduced in the proof of Lemma 8, consider the quadratic ODE (5.1); the two initial states of a qubit and defined in (5.5) and (5.6), respectively; and the terminal time defined in (5.13).
Suppose we have a quantum algorithm that, given a black box to prepare a state that is either or , can produce quantum states or that are within distance of and , respectively. Since by Lemma 8, and have constant overlap, the overlap between and is also constant for sufficiently small . More precisely, we have
| (5.18) |
by (5.15), which implies
| (5.19) |
We also have
| (5.20) |
and similarly for . These three inequalities give us
| (5.21) |
which is at least a constant for (say) .
Lemma 7 therefore shows that preparing the states and requires time , as these states can be used to distinguish the two possibilities with bounded error. By Lemma 8, this time is . This shows that we need at least exponential simulation time to approximate the solution of arbitrary quadratic ODEs to within sufficiently small bounded error when . ∎
Note that exponential time is achievable since our QCL algorithm can solve the problem by taking to be exponential in , where is the truncation level of Carleman linearization. (The algorithm of Leyton and Osborne also solves quadratic differential equations with complexity exponential in , but requires the additional assumptions that the quadratic polynomial is measure-preserving and Lipschitz continuous [39].)
6 Applications
Due to the assumptions of our analysis, our quantum Carleman linearization algorithm can only be applied to problems with certain properties. First, there are two requirements to guarantee convergence of the inhomogeneous Carleman linearization: the system must have linear dissipation, manifested by ; and the dissipation must be sufficiently stronger than both the nonlinear and the forcing terms, so that . Dissipation typically leads to an exponentially decaying solution, but for the dependency on and in (4.141) to allow an efficient algorithm, the solution cannot exponentially approach zero.
However, this issue does not arise if the forcing term resists the exponential decay towards zero, instead causing the solution to decay towards some non-zero (possibly time-dependent) state. The norm of the state that is exponentially approached can possibly decay towards zero, but this decay itself must happen slower than exponentially for the algorithm to be efficient.333Also note that the QCL algorithm might provide an advantage over classical computation for homogeneous equations in cases where only evolution for a short time is of interest.
We now investigate possible applications that satisfy these conditions. First we present an application governed by ordinary differential equations, and then we present possible applications in partial differential equations.
Several physical systems can be represented in terms of quadratic ODEs. Examples include models of interacting populations of predators and prey [51], dynamics of chemical reactions [44, 17], and the spread of an epidemic [11]. We now give an example of the latter, based on the epidemiological model used in [52] to describe the early spread of the COVID-19 virus.
The so-called SEIR model divides a population of individuals into four components: susceptible (), exposed (), infected (), and recovered (). We denote the rate of transmission from an infected to a susceptible person by , the typical time until an exposed person becomes infectious by the latent time , and the typical time an infectious person can infect others by the infectious time . Furthermore we assume that there is a flux of individuals constantly refreshing the population. This flux corresponds to susceptible individuals moving into, and individuals of all components moving out of, the population, in such a way that the total population remains constant.
To ensure that there is sufficiently strong linear decay to guarantee Carleman convergence, we also add a vaccination term to the component. We choose an approach similar to that of [53], but denote the vaccination rate, which is approximately equal to the fraction of susceptible individuals vaccinated each day, by . The model is then
| (6.1) | ||||
| (6.2) | ||||
| (6.3) | ||||
| (6.4) |
The sum of equations (6.1)–(6.4) shows that is a constant. Hence we do not need to include the equation for in our analysis, which is crucial since the component would have introduced positive eigenvalues. The maangces corresponding to (2.1) are then
| (6.5) | ||||
| (6.6) |
Since is a triangular matrix, its eigenvalues are located on its diagonal, so . Furthermore we can bound , , and , so
| (6.7) |
We see that the condition for guaranteed convergence of Carleman linearization is . Essentially, the Carleman method only converges if the (nonlinear) transmission is sufficiently slower than the (linear) decay.
To assess how restrictive this assumption is, we consider the SEIR parameters used in [52]. Note that they also included separate components for asymptomatic and hospitalized persons, but to simplify the analysis we include both of these components in the component. In their work, they considered a city with approximately inhabitants. In a specific period, they estimated a travel flux444Since we require that does not approach exponentially, we can assume that the travel flux is some non-zero, but small, value, e.g., individual per day. of individuals per day, latent time days, infectious time days, and transmission rate . We let the initial condition be dominated by the susceptible component so that , and we assume555This example arguably corresponds to quite rapid vaccination, and is chosen here such that remains smaller than one, as required to formally guarantee convergence of the Carleman method. However, as shown in the upcoming example of the Burgers equation, larger values of might still allow for convergence in practice, suggesting that our algorithm might handle lower values of the vaccination rate. that . With the stated parameters, a direct calculation gives , showing that the assumptions of our algorithm can correspond to some real-world problems that are only moderately nonlinear.
While the example discussed above has only a constant number of variables, 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, similar to what is done in [10] and [45].
Other examples of high-dimensional ODEs arise from the discretization of certain PDEs. Consider, for example, equations for of the type
| (6.8) |
with the forcing term being a function of both space and time. This equation can be cast in the form of (2.1) by using standard discretizations of space and time. Equations of the form (6.8) can represent Navier–Stokes-type equations, which are ubiquitous in fluid mechanics [38], and related models such as those studied in [36, 47, 50] to describe the formation of large-scale structure in the universe. Similar equations also appear in models of magnetohydrodynamics (e.g., [26]), or the motion of free particles that stick to each other upon collision [12]. In the inviscid case, , the resulting Euler-type equations with linear damping are also of interest, both for modeling micromechanical devices [5] and for their intimate connection with viscous models [24].
As a specific example, consider the one-dimensional forced viscous Burgers equation
| (6.9) |
which is the one-dimensional case of equation (6.8) with . Equation (6.9) is often used as a simple model of convective flow [14]. For concreteness, let the initial condition be on the domain , and use Dirichlet boundary conditions . We force this equation using a localized off-center Gaussian with a sinusoidal time dependence,666Note that this forcing does not satisfy the general conditions for efficient implementation of our algorithm since it is not sparse. However, we expect that the algorithm can still be implemented efficiently for a structured non-sparse forcing term such as in this example. given by . To solve this equation using the Carleman method, we discretize the spatial domain into points and use central differences for the derivatives to get
| (6.10) |
with . This equation is of the form (2.1) and can thus generate the Carleman system (3.2). The resulting linear ODE can then be integrated using the forward Euler method, as shown in Figure 1. In this example, the viscosity is defined such that the Reynolds number , and spatial discretization points were sufficient to resolve the solution. The figure compares the Carleman solution with the solution obtained via direct integration of (6.10) with the forward Euler method (i.e., without Carleman linearization). By inserting the matrices , , and corresponding to equation (6.10) into the definition of (2.2), we find that is indeed negative as required, given Dirichlet boundary conditions, but the parameters used in this example result in . Even though this does not satisfy the requirement of the QCL algorithm, we see from the absolute error plot in Figure 1 that the maximum absolute error over time decreases exponentially as the truncation level is incremented (in this example, the maximum Carleman truncation level considered is ). Surprisingly, this suggests that in this example, the error of the classical Carleman method converges exponentially with , even though .
7 Discussion
In this paper we have presented a quantum Carleman linearization (QCL) algorithm for a class of quadratic nonlinear differential equations. Compared to the previous approach of [39], our algorithm improves the complexity from an exponential dependence on to a nearly quadratic dependence, under the condition as defined in (2.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 . Furthermore, we established a lower bound showing that for general quadratic differential equations with , quantum algorithms must have worst-case complexity exponential in . 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 makes differential equations classically tractable. Clearly a naive integration of the truncated Carleman system (3.2) is not efficient on a classical computer since the system size is . But furthermore, it is unlikely that any classical algorithm for this problem can run in time polylogarithmic in . If we consider Problem 1 with dissipation that is small compared to the total evolution time, but let the nonlinearity and forcing be even smaller such that , 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 non-dissipative linear differential equations, which is a BQP-hard problem even when the dynamics are unitary [29]. 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 , 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 the problem can be solved efficiently by quantum computers in the worst case.
Furthermore, the complexity of our algorithm has nearly quadratic dependence on , namely . It is unknown whether the complexity for quadratic ODEs must be at least linear or quadratic in . Note that sublinear complexity is impossible in general because of the no-fast-forwarding theorem [8]. 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 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 (see Section 8 of [9] 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.
It is possible that variations of the Carleman linearization procedure could increase the accuracy of the result. For instance, instead of using just tensor powers of as auxiliary variables, one could use other nonlinear functions. Several previous papers on Carleman linearization have suggested using multidimensional orthogonal polynomials [6, 30]. They also discuss approximating higher-order terms with lower-order ones in (3.2) instead of simply dropping them, possibly improving accuracy. Such changes would however alter the structure of the resulting linear ODE, which could affect the quantum implementation.
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 3, which states the global error increases at most linearly with . To establish this result for the Euler method, it suffices to choose the time step (4.46) to ensure , and then estimate the growth of global error by (4.73). 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 Section 6 (see in particular Figure 1) suggest that the condition is not a strict requirement for the viscous Burgers equation, since we observe convergence even though . This suggests that some property of equation (6.9) makes it more amenable to Carleman linearization than our current analysis predicts. 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 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 in 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 [37], 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.
Acknowledgments
We thank Paola Cappellaro for valuable discussions and comments. JPL thanks Aaron Ostrander for inspiring discussions. HØK thanks Bernhard Paus Græsdal for introducing him to the Carleman method. HKK and NFL thank Seth Lloyd for preliminary discussions on nonlinear equations and the quantum linear systems algorithm. We also think Dong An and anonymous referees for pointing out the exponential decay of solutions to dissipative homogeneous equations. AMC and JPL did part of this work while visiting the Simons Institute for the Theory of Computing in Berkeley and gratefully acknowledge its hospitality.
AMC and JPL 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 National Science Foundation (CCF-1813814). HØK, HKK, and NFL were partially funded by award no. DE-SC0020264 from the Department of Energy.
References
- [1] Scott Aaronson, NP-complete problems and physical reality, ACM SIGACT News 36 (2005), no. 1, 30–52, arXiv:quant-ph/0502072.
- [2] Daniel S. Abrams and Seth Lloyd, Nonlinear quantum mechanics implies polynomial-time solution for NP-complete and #P problems, Physical Review Letters 81 (1998), no. 18, 3992, arXiv:quant-ph/9801041.
- [3] Roberto F. S. Andrade, Carleman embedding and Lyapunov exponents, Journal of Mathematical Physics 23 (1982), no. 12, 2271–2275.
- [4] Kendall E. Atkinson, An introduction to numerical analysis, John Wiley & Sons, 2008.
- [5] Min-Hang Bao, Micro mechanical transducers: pressure sensors, accelerometers and gyroscopes, Elsevier, 2000.
- [6] Richard Bellman and John M. Richardson, On some questions arising in the approximate solution of nonlinear differential equations, Quarterly of Applied Mathematics 20 (1963), 333–339.
- [7] Dominic W. Berry, High-order quantum algorithm for solving linear differential equations, Journal of Physics A: Mathematical and Theoretical 47 (2014), no. 10, 105301, arXiv:1010.2745.
- [8] Dominic W. Berry, Graeme Ahokas, Richard Cleve, and Barry C. Sanders, Efficient quantum algorithms for simulating sparse Hamiltonians, Communications in Mathematical Physics 270 (2007), 359–371, arXiv:quant-ph/0508139.
- [9] Dominic W. Berry, Andrew M. Childs, Aaron Ostrander, and Guoming Wang, Quantum algorithm for linear differential equations with exponentially improved dependence on precision, Communications in Mathematical Physics 356 (2017), no. 3, 1057–1081, arXiv:1701.03684.
- [10] Derdei Bichara, Yun Kang, Carlos Castillo-Chávez, Richard Horan, and Charles Perrings, SIS and SIR epidemic models under virtual dispersal, 03 2015.
- [11] Fred Brauer and Carlos Castillo-Chavez, Mathematical models in population biology and epidemiology, vol. 2, Springer, 2012.
- [12] Yann Brenier and Emmanuel Grenier, Sticky particles and scalar conservation laws, SIAM Journal on Numerical Analysis 35 (1998), no. 6, 2317–2328.
- [13] Roger Brockett, The early days of geometric nonlinear control, Automatica 50 (2014), no. 9, 2203–2224.
- [14] Johannes M. Burgers, A mathematical model illustrating the theory of turbulence, Advances in Applied Mechanics, vol. 1, Elsevier, 1948, pp. 171–199.
- [15] Yudong Cao, Anargyros Papageorgiou, Iasonas Petras, Joseph Traub, and Sabre Kais, Quantum algorithm and circuit design solving the Poisson equation, New Journal of Physics 15 (2013), no. 1, 013021, arXiv:1207.2485.
- [16] Torsten Carleman, Application de la théorie des équations intégrales linéaires aux systèmes d’équations différentielles non linéaires, Acta Mathematica 59 (1932), no. 1, 63–87.
- [17] Alessandro Ceccato, Paolo Nicolini, and Diego Frezzato, Recasting the mass-action rate equations of open chemical reaction networks into a universal quadratic format, Journal of Mathematical Chemistry 57 (2019), 1001–1018.
- [18] Andrew M. Childs, Robin Kothari, and Rolando D. Somma, Quantum algorithm for systems of linear equations with exponentially improved dependence on precision, SIAM Journal on Computing 46 (2017), no. 6, 1920–1950, arXiv:1511.02306.
- [19] Andrew M. Childs and Jin-Peng Liu, Quantum spectral methods for differential equations, Communications in Mathematical Physics 375 (2020), 1427–1457, arXiv:1901.00961.
- [20] Andrew M. Childs, Jin-Peng Liu, and Aaron Ostrander, High-precision quantum algorithms for partial differential equations, 2020, arXiv:2002.07868.
- [21] Andrew M. Childs and Joshua Young, Optimal state discrimination and unstructured search in nonlinear quantum mechanics, Physical Review A 93 (2016), no. 2, 022314, arXiv:1507.06334.
- [22] B. David Clader, Bryan C. Jacobs, and Chad R. Sprouse, Preconditioned quantum linear system algorithm, Physical Review Letters 110 (2013), no. 25, 250504, arXiv:1301.2340.
- [23] Pedro C. S. Costa, Stephen Jordan, and Aaron Ostrander, Quantum algorithm for simulating the wave equation, Physical Review A 99 (2019), no. 1, 012323, arXiv:1711.05394.
- [24] Constantine M. Dafermos, Hyperbolic conservation laws in continuum physics, vol. 3, Springer, 2005.
- [25] Germund G. Dahlquist, A special stability problem for linear multistep methods, BIT Numerical Mathematics 3 (1963), no. 1, 27–43.
- [26] P. A. Davidson, An introduction to magnetohydrodynamics, Cambridge Texts in Applied Mathematics, Cambridge University Press, 2001.
- [27] Ilya Y. Dodin and Edward A. Startsev, On applications of quantum computing to plasma simulations, 2020, arXiv:2005.14369.
- [28] Alexander Engel, Graeme Smith, and Scott E. Parker, Quantum algorithm for the Vlasov equation, Physical Review A 100 (2019), no. 6, 062315, arXiv:1907.09418.
- [29] Richard P. Feynman, Quantum mechanical computers, Optics News 11 (1985), 11–20.
- [30] Marcelo Forets and Amaury Pouly, Explicit error bounds for Carleman linearization, 2017, arXiv:1711.02552.
- [31] Alfredo Germani, Costanzo Manes, and Pasquale Palumbo, Filtering of differential nonlinear systems via a Carleman approximation approach, Proceedings of the 44th IEEE Conference on Decision and Control, pp. 5917–5922, 2005.
- [32] Aram W. Harrow, Avinatan Hassidim, and Seth Lloyd, Quantum algorithm for linear systems of equations, Physical Review Letters 103 (2009), no. 15, 150502, arXiv:0811.3171.
- [33] Carl W. Helstrom, Quantum detection and estimation theory, Journal of Statistical Physics 1 (1969), 231–252.
- [34] Ilon Joseph, Koopman-von Neumann approach to quantum simulation of nonlinear classical dynamics, 2020, arXiv:2003.09980.
- [35] Edward H. Kerner, Universal formats for nonlinear ordinary differential systems, Journal of Mathematical Physics 22 (1981), no. 7, 1366–1371.
- [36] Lev Kofman, Dmitri Pogosian, and Sergei Shandarin, Structure of the universe in the two-dimensional model of adhesion, Monthly Notices of the Royal Astronomical Society 242 (1990), no. 2, 200–208.
- [37] Krzysztof Kowalski and Willi-Hans Steeb, Nonlinear dynamical systems and Carleman linearization, World Scientific, 1991.
- [38] Pierre Gilles Lemarié-Rieusset, The Navier-Stokes problem in the 21st century, CRC Press, 2018.
- [39] Sarah K. Leyton and Tobias J. Osborne, A quantum algorithm to solve nonlinear differential equations, 2008, arXiv:0812.4423.
- [40] Noah Linden, Ashley Montanaro, and Changpeng Shao, Quantum vs. classical algorithms for solving the heat equation, arXiv:2004.06516.
- [41] Seth Lloyd, Giacomo De Palma, Can Gokler, Bobak Kiani, Zi-Wen Liu, Milad Marvian, Felix Tennie, and Tim Palmer, Quantum algorithm for nonlinear differential equations, 2020, arXiv:2011.06571.
- [42] Gerasimos Lyberatos and Christos A. Tsiligiannis, A linear algebraic method for analysing Hopf bifurcation of chemical reaction systems, Chemical Engineering Science 42 (1987), no. 5, 1242–1244.
- [43] Ashley Montanaro and Sam Pallister, Quantum algorithms and the finite element method, Physical Review A 93 (2016), no. 3, 032324, arXiv:1512.05903.
- [44] Elliott W. Montroll, On coupled rate equations with quadratic nonlinearities, Proceedings of the National Academy of Sciences 69 (1972), no. 9, 2532–2536.
- [45] Rany Qurratu Aini, Deden Aldila, and Kiki Sugeng, Basic reproduction number of a multi-patch SVI model represented as a star graph topology, 10 2018, p. 020237.
- [46] Andreas Rauh, Johanna Minisini, and Harald Aschemann, Carleman linearization for control and for state and disturbance estimation of nonlinear dynamical processes, IFAC Proceedings Volumes 42 (2009), no. 13, 455–460.
- [47] Sergei F. Shandarin and Yakov B. Zeldovich, The large-scale structure of the universe: turbulence, intermittency, structures in a self-gravitating medium, Reviews of Modern Physics 61 (1989), no. 2, 185–220.
- [48] Vivek V. Shende, Stephen S. Bullock, and Igor L. Markov, Synthesis of quantum-logic circuits, IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems 25 (2006), no. 6, 1000–1010, arXiv:quant-ph/0406176.
- [49] Willi-Hans Steeb and F. Wilhelm, Non-linear autonomous systems of differential equations and Carleman linearization procedure, Journal of Mathematical Analysis and Applications 77 (1980), no. 2, 601–611.
- [50] Massimo Vergassola, Bérengère Dubrulle, Uriel Frisch, and Alain Noullez, Burgers’ equation, devil’s staircases and the mass distribution for large-scale structures, Astronomy and Astrophysics 289 (1994), 325–356.
- [51] Paul Waltman, Competition models in population biology, SIAM, 1983.
- [52] Chaolong Wang, Li Liu, Xingjie Hao, Huan Guo, Qi Wang, Jiao Huang, Na He, Hongjie Yu, Xihong Lin, An Pan, Sheng Wei, and Tangchun Wu, Evolving epidemiology and impact of non-pharmaceutical interventions on the outbreak of coronavirus disease 2019 in Wuhan, China, Journal of the American Medical Association 323 (2020), no. 19, 1915–1923.
- [53] Gul Zaman, Yong Han Kang, and Il Hyo Jung, Stability analysis and optimal vaccination of an SIR epidemic model, Biosystems 93 (2008), no. 3, 240–249.