Quantum Algorithms for Solving Ordinary Differential Equations via Classical Integration Methods

Identifying computational tasks suitable for (future) quantum computers is an active field of research. Here we explore utilizing quantum computers for the purpose of solving differential equations. We consider two approaches: (i) basis encoding and fixed-point arithmetic on a digital quantum computer, and (ii) representing and solving high-order Runge-Kutta methods as optimization problems on quantum annealers. As realizations applied to two-dimensional linear ordinary differential equations, we devise and simulate corresponding digital quantum circuits, and implement and run a 6$^{\mathrm{th}}$ order Gauss-Legendre collocation method on a D-Wave 2000Q system, showing good agreement with the reference solution. We find that the quantum annealing approach exhibits the largest potential for high-order implicit integration methods. As promising future scenario, the digital arithmetic method could be employed as an"oracle"within quantum search algorithms for inverse problems.


Introduction
Solving differential equations (DEs) is a ubiquitous task in the scientific and engineering community. Traditional numerical algorithms need to integrate the time steps sequentially due to their interdependence, which is at odds with the trend towards parallelization in modern highly-parallel high-performance computing (HPC) architectures. This trend already motivated research in various directions, such as parallelin-time algorithms [13]. An alternative route could emerge from quantum computers, where we could benefit from better asymptotic scaling properties compared to classical computers when applied to certain tasks [2,15,31,33]. Here we explore approaches of utilizing quantum computers for the purpose of solving DEs: (i) using basis encoding to describe the Benjamin Zanger: benjamin.zanger@tum.de Christian B. Mendl: christian.mendl@tum.de Martin Schulz: schulzm@in.tum.de Martin Schreiber: martin.schreiber@tum.de state of the differential equation, and mimicking arithmetic operations of classical algorithms on a digital quantum computer (Section 3); and (ii) reformulating the time integration as a quantum annealing problem (Section 4).
To realize (i), we will first design quantum circuits for arithmetic operations in numerical fixed-point representation, and then apply these to solve linear ordinary differential equations (ODEs) in two dimensions.
Our second approach (ii) reformulates Runge-Kutta integration methods (both explicit and implicit) as minimization problems suitable for a quantum annealing framework. We additionally design a flexible number representation to reach high accuracy. We then run the annealing task (based on a method of order six) on a D-Wave 2000Q system, and are able to demonstrate a good agreement with the reference solution obtained on a classical computer.
Related work includes proposals for quantum algorithms that are capable of solving both linear partial differential equations (PDEs) [7,18,39] and nonlinear PDEs [22,25,26]. While for linear differential equations an exponential advantage in the resources is known for some time, exponential advantages for nonlinear differential equations has only been found recently. Both algorithms, for the linear and nonlinear case, assume that the state of the PDE is amplitude encoded (i.e., via the amplitudes of the quantum wavefunction), which allows a logarithmic scaling of quantum resources for an increase in the dimension of the state space. However, this also requires state preparation for the input and quantum state tomography for measuring the output, which is imposing scaling issues. There have been efforts to improve these scaling for certain PDEs, e.g. for the wave equation in [7]. Liu et al. [22] approached this problem by embedding the initial state of nonlinear PDEs into a higher dimensional state and to derive a measurement success probability of their algorithm, enabling a better scaling for the readout with the use of amplitude amplification. Unfortunately, their algorithm only works for dissipative differential equations. To our knowledge, there currently does not exist any algorithm which solves the state preparation and readout problem for general PDEs.
In contrast, our work does not provide a logarithmic scaling of quantum resources with the state space but does not obey the state preparation and readout problem. For our second approach (ii), the time complexity for a single time step in our algorithm does not scale with the state space. Potential applications include all kind of linear and nonlinear ODEs which are expressible with low order polynomial terms, including problems from biology and fluid dynamics.

Classical Time Integration Methods
This section recapitulates the numerical time integration schemes that are used in the quantum algorithms discussed later. Specifically, we will discuss the Runge-Kutta methods and how they can be phrased as optimization problems. Since we can only cover a few aspects of time integration, we solely target single-step methods and refer interested readers to Durran [9], Hairer et al. [16,17] and Shu and Osher [32].
Concretely, we consider a system of ordinary differential equations written as with u(t) ∈ R N for all t ≥ 0 and the initial condition The Runge-Kutta (RK) formulation unifies various kinds of explicit and implicit time integration schemes. Each particular RK method can be described by a Butcher table, consisting of a matrix A ∈ R s×s and two vectors, b, c ∈ R s , where s is called the number of stages of the method [4]. The RK formulation for a time step ∆t is then given by the equations (i ∈ {1, . . . , s}) withũ(t+∆t) being the approximated solution at the next time step. Explicit methods correspond to a strictly lowerdiagonal matrix A, hence avoiding any implicit dependencies. Otherwise, the method is implicit, with coefficients, e.g., given by the collocation method [17] that leads to a dense matrix A. Details are skipped here for sake of brevity, but we point out the large challenge to cope with multiple stages depending implicitly on each other and various associated iterative algorithms [10,11,28].
Given the RK formulation and the demand for higher-order implicit time integration methods, efficient solvers are required to cope with the implicit dependencies in these equations. To provide one example, a class of prominent solvers are based on spectral deferred corrections [10], which avoid the implicit dependency of multiple stages. Here, the underlying idea is to iteratively correct an approximation of the solution, while each correction only requires evaluating forward and backward Euler time steps. However, the underlying procedure is still iterative, requiring multiple iterations to gain higher-order accuracy, and hence is computationally demanding. Quantum computers could provide a new approach to solve such problems efficiently as a combinatorial problem, as we will elaborate in Section 4.

Digital Quantum Circuit Time Integration
In this section, we investigate time integration by utilizing qubit-based arithmetic operations on digital quantum computers. Realizing integer operations via quantum circuits is well-known in the literature [8,36]. Our contribution is a generalization to fixed-point number representations.

Arithmetic on Digital Quantum Computers
Guiding principles for the following identities can be gleaned from analogies to continuous variable quantum computing [3,24,38]. In this framework, a quantum register |x stores a real number x ∈ R. This register is acted on by the position operatorX , formally defined asX = and its conjugate momentum operatorP. They obey the canonical commutation relations [X ,P] = i and are related viaP =F †XF , whereF is the Fourier transformation.
The conjugated variables can be used to realize addition and subtraction [23]. Applying the Hamilto-nianĤ =P for a time c leads to a shift in the variable by c, since (setting = 1) Thus the addition can be expressed as time evolution. Given a continuous variable quantum state |a , a ∈ R, we get e −icP |a = |a + c .
Our goal is to approximate this system using n qubits, i.e., 2 n available states. As a first step, we describe integer arithmetic and then investigate the required modifications for a fixed-point representation. Given a ∈ {0, . . . , 2 n − 1}, the quantum register state |a is canonically identified by the corresponding tensor product of single qubit states based on the binary representation of a, i.e., |a = |a n−1 · · · |a 1 |a 0 for a = a n−1 · · · a 1 a 0 . We will encounter the discrete Fourier transform at several occasions, for which we use the convention e −2πijk/2 n |k ∀j = 0, . . . , 2 n − 1.
(8) F andF † can be efficiently implemented using the quantum Fourier transform (QFT) algorithm.
With these preparations, we define a discretized position operator byX The discrete momentum operator is related toX bŷ analogous to the continuous variable case.
Eq. (7) holds analogously for discrete registers [37]: given integers a and c, To realize this operation on a digital quantum computer, first note that, based on Eq. (10), e −2πicX/2 n can be evaluated as follows, using thatX is a diagonal matrix with respect to the computational basis: The unitary matrix (qubit gate) representation of the inner operation is R n−k (c) † , with the definition for any integer m ≥ 0 and ϑ ∈ R. In the following, we will abbreviate R m (1) by R m . Assume we want to add two integers a, b ∈ {0, . . . , 2 n − 1}, which are both stored in quantum registers, i.e., the initial state is |a, b : according to Verdon et al. [37], a "von Neumann measurement" of |a combined with (11) can be used to increment the second register by a: We can decompose the operator on the left as e −2πiX⊗P /2 n = I ⊗F † e −2πiX⊗X/2 n I ⊗F , (16) and rewrite the inner operator as From this representation, one observes that R m gates controlled by |a and applied to the |b register in case a j = 0, b k = 0 are sufficient; see Fig. 1 for an illustration. Figure 1: An addition circuit for two quantum registers with n = 2, implementing Eq. (15); also compare with [8].
For general n, the operator in Eq. (17) requires 1 2 n(n + 1) controlled-R m gates. Since the number of gates for the quantum Fourier transform has an O(n 2 ) scaling as well, the overall cost for the addition in Eq. (15) is O(n 2 ) two-qubit gates (controlled phase rotations).
Analogous to the addition circuits, subtraction can be realized by letting c → −c mod 2 n , and noting that similar to Eq. (15). Thus it suffices to take the adjoint of Eqs. (16) and (17) to implement subtraction. As a remark, the principle behind (15) can be generalized as follows: given a map g : Then e −2πiŶg⊗P /2 n |a, b = |a, g(a) + b mod 2 n . (20) The final operation we discuss here is multiplication: given a, b, c ∈ {0, . . . , 2 n − 1}, it holds that e −2πiX⊗X⊗P /2 n |a, b, c = |a, b, ab + c mod 2 n , (21) again analogous to Eq. (15). We can decompose and represent the inner operator as All the exponential functions with j + k + ≥ n, a j = 0, b k = 0, or c = 0 evaluate to 1 and need not be taken into account explicitly. It suffices to act with R m gates defined in Eq. (14) on the |c register, controlled by both |a and |b ; Fig. 2 illustrates this procedure for n = 2. The overall number of doubly-controlled R m gates is 1 6 n(2+3n+n 2 ), when taking the condition j+k+ < n into account. This O(n 3 ) scaling, together with the required doubly-controlled gates, which have to be emulated by single and two qubit gates on current digital quantum computers, pose a considerable practical limitation on the present approach.
The advertised fixed-point representation with scaling factor 2 −q for an integer q ≥ 0 amounts to reinterpreting a quantum state |a as representing the number 2 −q a, written in binary representation as 2 −q a = a n−1 · · · a q . a q−1 · · · a 0 (24) (q digits after the dot). The two's complement of binary numbers is suitable for including negative numbers as well. Then the quantum state |a for a ∈ {0, . . . , 2 n − 1} represents 2 −q −2 n−1 a n−1 + 2 n−2 a n−2 + · · · + a 0 on the logical level. The representable numbers are thus (26) and the arithmetic operations are understood modulo 2 n−q . It is important to note that the operators and circuits for addition and subtraction derived so far remain exactly the same. Concerning multiplication, the exact product of two numbers from D n,q requires (in general) 2q digits after the dot in binary format. To cast this back into an element of D n,q , one disregards the trailing q digits, i.e., rounds the number. In binary representation, a bit shift to the right by q places realizes this operation for non-negative numbers. To implement this procedure using quantum circuits, we first introduce the following operator (acting on two n-qubit registers) where denotes the (logical) right bit shift. Then As before, we can diagonalize the operator by Fourier transformation applied to the |c register. The resulting diagonal operator then reads (compare with Eq. (23)): The bit shift corresponds to the 2 −q factor in the exponent, and the rounding to the condition j + k ≥ q. Finally, we describe how integer division by 2 and subsequent rounding towards −∞ is achievable as quantum circuit using the two's complement representation. First note that 1 2 −2 n−1 a n−1 + 2 n−2 a n−2 + · · · + a 0 = −2 n−1 a n−1 + 2 n−2 a n−1 + 2 n−3 a n−2 + · · · + 2 −1 a 0 .
a n−1 now appears twice, and rounding corresponds to dropping the term 2 −1 a 0 . The quantum circuit shown in Fig. 3 realizes this procedure, illustrated for n = 4.

Demonstration and Results
We employ the explicit Euler scheme as a demonstration for digital quantum circuit time integration. In principle, the approach also works for other higherorder explicit Runge-Kutta schemes since solely elementary arithmetic operations are required, and these operations can be directly mapped to quantum gates. Specifically, we consider the coupled linear differential equations (N = 2 dimensions) d dt with t ≥ 0, u 0 = (x 1 , x 2 ). We use this system as a representation of a semi-discrete hyperbolic PDE.
The analytical solution of this equation is given by To implement Eq. (31), we need a circuit with quantum registers initialized toũ 1 (t) andũ 2 (t). The output of the circuit are registers storingũ 1 (t + ∆t) and u 2 (t + ∆t). We also need two "ancilla" temporary registers for evaluating f (ũ(t)). The overall circuit is depicted in Fig. 4. Note that applying a sequence of time steps requires a reinitialization of the ancilla registers to |0 ; alternatively, this could be achieved by "uncomputing" them.
To implement the circuit, ∆t has to be a power of 1 2 , which allows us to save qubits during the simulation. This is solely for efficiency reasons and other ∆t values could be realized as well. Instead of multiplying by ∆t, we can divide multiple times by 2. Note, that ∆t does not need to be stored in a quantum register. Fig. 5 shows the circuits implementing the right side f (u) of the differential equation (31).
We simulate the quantum circuits using Qiskit [1], without noise or decoherence. The time step is chosen as ∆t = 1 2 . Each number register consists of 4 qubits, and the fixed-point arithmetic scaling factor is set to q = 1.
As verification, we implement an explicit Euler also on a classical computer and compare the results. In all cases, the results match in a numerical sense, as expected.
We first study the stationary case (du/dt = 0) at the equilibrium point u(t) = (0, 0), see Fig. 6a. Although this seems to be trivial, such stationary cases play an important role for, e.g., more complex cases like geostrophic balance in atmospheric simulations [34], and therefore are included here as proof-of-concept. Indeed the DQC solution remains constant as well, as expected.
Next, we simulate the same differential equation, but using initial values u 0 = (0, −1); the results are visualized in Fig. 6b. The simulated solution amplifies the real one, as expected for the explicit Euler scheme [9]: since the digital quantum computer mimics the operations of a classical solver, the achievable accuracy is likewise bounded by the numerical errors of the classical method. Our demonstrations mainly serve as proof-of-principle.

Quantum Annealing Time Integration
Annealing is a method for solving optimization problems. The solution of a problem is encoded into the global minimum ("ground state") of an energy function. One possible encoding, relevant for combinatorial optimization problems, is quadratic unconstrained binary optimization (QUBO), see, e.g., Glover et al. [14]. In this formulation, the goal is to minimize the binary variables σ i ∈ {0, 1} of the energy function (Ising-type Hamiltonian) for given parameters J i,j , h i ∈ R. While thermal annealing uses thermal fluctuations to overcome local minima of the target function, quantum annealing uses quantum tunneling effects for that purpose, potentially exhibiting a faster convergence [20]. Nevertheless, neither thermal nor quantum annealing are guaranteed to find the ground state of the optimization problem.

Time Integration as an Optimization Problem
In the following, we consider a system of autonomous linear differential equations and M is the highest polynomial degree of the differential equation. An extension to other nonlinear equations is possible with our framework, as well, as long as the action of f can be formulated as part of a QUBO problem. A time step of collocation-based implicit Runge-Kutta time integration [17] can be represented in the form ofũ (t + ∆t) =ũ(t) + ∆tSf ũ(t + ∆t) Mul ∆t calculate f 2 (ũ) |∆tf 2 (ũ(t)) ũ 1 (t + ∆t) u 2 (t + ∆t) Figure 4: Circuit implementation of the explicit Euler scheme for Eq. (31). The input of the circuit are two quantum registers, initialized with the values ofũ1(t) andũ2(t) in fixed-point arithmetic. The circuit calculates one iteration of the explicit Euler scheme and outputsũ1(t + ∆t) andũ2(t + ∆t). The green region marks ancilla registers. withũ(0) a suitably chosen vector with the initial conditions, and S the spectral integration matrix given by the particular quadrature method used in the collocation formulation. In this formulation, the vector u also stores all intermediate state solutions at the quadrature points. We can interpret (35) as an optimization problem (see, e.g., [11]) of the form In the following, we will phrase this optimization problem in the context of a quantum annealer, to obtainũ(t + ∆t) for a givenũ(t). A similar approach for solving linear systems of equations has already been introduced by O'Malley and Vesselinov in [29]. Our approach for the extension to polynomial equations is mostly inspired by Cheng Chang et al [6]. Based on the Runge-Kutta methods (see Eqs. (2) and (3)) and using Eq. (34), we derivẽ The variables to optimize areũ j (t m+1 ) and K oj for o = 1, . . . , s, j = 1, . . . , N . The vectorũ j (t m ) remains constant in Eq. (39) and, therefore, is not part of the optimization. Unfortunately, we cannot directly map this equation to the Ising model because the parameters to optimize are real. To overcome this issue, we develop (another) tailored approximation using binary variables, as shown below.

Binary Number Representation
Different from the fixed-point representation described in Section 3.1, here our goal is an adaptive version (similar to the floating-point format). We approximate a real number g ∈ R by where d, k ∈ R are parameters and σ i (for i = 0, . . . , n − 1) the available binary variables. During an optimization step on a quantum computer, only the binary variables (σ i ) are modified. The detailed procedure will be described below.

Variational Approach
One limitation of current quantum annealers is the connectivity between qubits available in hardware. In order to use arbitrary couplings in the QUBO model, embeddings have to be used. An embedding maps a QUBO problem onto a slightly different one. This new model solves the same task using fewer connections per qubit, but a larger overall number of qubits [21]. Thus, in general, an important goal is the design of algorithms with sparse connectivity and small qubit counts. Typically, one chooses a variational approach to achieve this [27]. The key idea is to only run a part of the problem on the quantum device and perform post-processing on a classical computer to optimize parameters. This procedure is repeated until the output converges to the final solution. However, convergence using the post-processing technique is only guaranteed in the convex case, which in general is only satisfied for linear actions f . In the following, we analyze the number of qubits and connectivity required by our algorithm to achieve a given numerical accuracy. Subsequently, we introduce a variational formulation.

Numerical Accuracy and Qubit Connectivity
We assume that n qubits are sufficient to approximate the range of values needed to represent an entry of the solution vector. By choosing k = n and d = 0 in Eq. (40), this interval is [0, 1], and the accuracy for each simulated number is given by 2 −n . Furthermore, we assume that the solution of the DEQ remains in this interval.
The overall number of required qubits for s stages in the Runge-Kutta method is n(s + 1), since an additional n qubits store the resultũ j (t + ∆t) at the next time step. Additionally, auxiliary qubits are needed to reduce cubic and higher terms to quadratic ones. To achieve this, we use reduction by substitution which was first presented by Rosenberg [30]. For this method, one auxiliary qubit per reduction is used. The amount of auxiliary qubits rises exponentially with the order M , but is also limited by the number of stages s and the system size of the DEQ. For the remainder of this analysis, we will only assume linear DEQ (M = 1) for which no auxiliary qubits are needed.
Next, we determine the required connectivity between qubits. According to Eq. (39), a multiplication between two registers requires n connections between the qubits and a multiplication with itself n − 1 qubits. Eq. (39) also shows that the connectivity depends both on the density of the Runge-Kutta table A, as well as the DEQ matrix L. For our analysis, we consider the worst-case scenario of dense matrices.
There are three sets of terms in Eq. (39) in which two different registers are multiplied. The first term is the multiplication ofũ j (t m+1 ) with itself (directly after the second equal sign in (39)), for which n − 1 connections are needed. In the second set of terms, the entries of K are multiplied with each other (several occasions after the second equal sign). Since there are sN of these entries, where N is the system size, we count sN (sN −1)/2 multiplications in total. However, sN of these multiplications are with the same registers. Therefore, there are sN (sN − 1)/2 − sN = s 2 N 2 /2 − 3sN/2 multiplications with other and sN multiplications with the same register. As a consequence, the total number of connections needed for the second scenario is s 2 N 2 n/2 − 3sN n/2 + sN (n − 1) = s 2 N 2 n/2 − sN n/2 − sN . In the third set of terms, all entries of K are multiplied with the entries ofũ j (t m+1 ). In total, these are sN multiplications, accounting for sN n connections. All of the three sets of terms create new connections, and therefore the connection counts must be summed up, leading to a total number of connections required for this algorithm (in the worst case) of (connect.) = n − 1 + s 2 N 2 n/2 + sN n/2 − sN. (41)

Variational Adaptive Number Representation
According to the s 2 N 2 n/2 term in Eq. (41), increasing the number of qubits n for approximating a real number will demand a much higher connectivity on the quantum device. To mitigate this issue, we switch to an adaptive number representation, as already anticipated in Eq. (40). Specifically, for each individual number g appearing in the algorithm, we iteratively refine k and d in Eq. (40), denoted g i , k i and d i for the i th iteration. As update step, we set with a fixed constant c > 0. The idea, reminiscent of nested intervals, is to shrink the range of values covered by 2 −k σ (which is optimized by the quantum computer) in Eq. (40), while absorbing the best candidate for g into the classical offset parameter d. The shift by 2 n−ki+1−1 in Eq. (43) ensures that both positive and negative corrections are possible. We remark that, as our experiments in the next section shows, c must not be chosen too large on real quantum annealers, and should ideally correspond to the convergence rate of the iteration.

Time Complexity
In the following, we analyze the time complexity of our algorithm. Since, to our knowledge, there is neither any evidence that quantum annealers have an advantage compared to classical computers, nor we are aware of any method to determine the time complexity of an algorithm on a quantum annealer, we study the time complexity of the algorithm on an adiabatic quantum computer. In this type of quantum computing, the solution of a problem is found by evolving a time-dependent Hamiltonian where the solution is encoded into the ground state of the final Hamiltonian.
The system is first initialized to the ground state of an initial Hamiltonian. If the evolution is slow enough, at the end of the evolution the system will remain in the state of the desired solution. The time the evolution takes scales with T = O 1 g 2 min , where g min is the smallest energy gap between the ground state and another state during the evolution process [12].
We choose an initial Hamiltonian H init with a sufficient high energy gap, so that it does not limit the time complexity of our algorithm. The timedependent Hamiltonian we use is given bŷ whereĤ QUBO is given by the QUBO Hamiltonian of our minimization function. Although the real minimal gap energy gap g min is found by determining the energy difference of the ground and first excited state for allĤ(t) ∀t, we determine this gap for the QUBO HamiltonianĤ QUBO in order to get an idea how this gap changes with parameters we chose in the algorithm. Afterwards, we compare our findings with numerical results of the minimal gap inĤ(t).
For linear ODEs the QUBO Hamiltonian is strictly convex. This makes it possible to calculate boundaries on the energy gap under the assumption that the lowest energy is 0. We first show the calculation in detail for the explicit Euler of a one dimensional ODE and later deal with the generalized Runge-Kutta case of linear ODEs. For the explicit Euler of a linear ODE with L (2) = Λ, a time step is given as In our Hamiltonian, u(t) is exactly represented, while u(t + ∆t) is discretized and represented by different quantum states. Therefore, we use the notation u(t + ∆t) = u n+1 ∆x, where u n+1 ∈ Z and ∆x is the smallest number resolution in our algorithm. The solution for u n+1 which minimizes Eq. (39) is given by We assume that u * n+1 ∈ Z, which means that a state for the energy of 0 exists. Because of the construction, we know that this is the lowest energy state and that the minimization function is symmetric. Therefore, the second lowest energy level is given by To find the energy gap g min , we can simply subtract the energy of u n+1 by the energy of u * n+1 . Our finding is that Next, we compare this result to the scaling of the minimal gap inĤ(t), which we computed numerically using and initial state ofĤ init = −0.
The comparison is depicted in Fig. 7. For sufficient small ∆x our previously determined scaling and numerical results perfectly match. Another important parameter to choose is ∆t. According to our derivation forĤ QUBO , a change in ∆t should not have any influence in g min . The comparison of the numerical and expected results confirms that changes in ∆t do not have an influence on g min . Our assumption that the scaling of the energy gap of the ground and first excited state ofĤ QUBO is comparable to the overall minimal energy gap g min agrees with numerical results. Therefore, we conclude that the time complexity of the explicit Euler using our algorithm approximately is Next, we determine the complexity of a generalized Runge-Kutta method. In this case, we have to discretize the stages as well, We refer to the minimizer of Eq. (39) using the * symbol and assume that these are contained in Z. Let us define the deviation from a solution K o,n and u n+1 to these minimizers as ∆K o,n = K o,n − K * o,n and ∆u n+1 = u n+1 − u * n+1 . With these definitions, the minimal energy gap in the QUBO Hamiltonian is calculated by where ∆u n+1 + s o=1 ∆K o,n ≥ 1. First, note that ∆u n+1 , ∆K 1,n , ..., ∆K s,n ∈ {0, 1} in order to minimize Eq. (50). In the following, we want to give a lower bound on the energy gap using only the first term of Eq. (50). Under the assumption, that ∆t < 1, this term is either minimized if all variables deviate or only a single of the stages. Therefore, where b = min{b i for 1 ≤ i ≤ s}. We verify that the bound also holds forĤ(t) by evaluating the minimal gap numerically. As a Runge-Kutta scheme, we use Crank-Nicolson. The scaling in ∆x for ∆t = 0.1 and ∆t = 10 −6 is depicted in Fig. 8. For this specific case, it can be seen that the minimal gap has a better scaling, g min ∝ O (∆x), compared to the lower bound on g QUBO we found. For our concrete example, we can also see that higher time steps result in a lower energy gap.
To get a scaling including the amount of stages, we limit the Runge-Kutta scheme on the collocation method using Chebyshev-Gauss quadrature in the following. This allows b = 1/s for increasing amount of stages. From the boundary we found in Eq. (51), we conclude that the time complexity of the Runge-Kutta method, using Chebyshev-Gauss collocation, is approximately better than For both the explicit Euler case, as well as the generalized Runge-Kutta case, it is possible to generalize the results to systems of equations of arbitrary sizes. This is since the term we used in the minimization function is independent of any mutual interference for other variables in a system of equations. Compared to a classical computer, where so far no parallel processing without an additional time scaling is known, we don't have that in this quantum algorithm. Instead, a larger system of equations only allocates more hardware/qubit resources.

Results of the QA Time Integration
We perform the quantum annealing task on the D-Wave 2000Q system, which features 2041 qubits. To encode the QUBO model onto the annealer, the values of the J matrix and the h vector are first normalized to be within the interval [−2, 1] for the J matrix and [−2, 2] for the h vector. For our problem, we directly map the values of J and h onto the interval [−1, 1]. Furthermore, as arbitrary connectivity between all qubits is not possible, the concrete connections available on the target annealing hardware have to be taken into account and the QUBO model has to be adapted to those hardware characteristics. If the QUBO problem needs more connections than the hardware provides, an additional embedding has to be found. We solve both of these tasks using the library minorminer, which uses a method given by Cai et al. [5]. We perform 100 reads on the annealer per iteration. For the last example, we additionally used the library PyQUBO [35] which maps higher order polynomials to quadratic ones using the reduction by substitution method.   9 shows the analytical and annealed result of the linear coupled differential equation given in Eq. (31) with the initial conditions u 1 (0) = 1 and u 2 (0) = 0. As a specific integration method, we choose the Gauss-Legendre collocation method of order 6 with a timestep of ∆t = 0.5. For every grid point, 3 qubits are used. We perform 15 iterations per timestep, starting with k 0 = 1 and c = 0.5. The numerical result obtained on a digital computer, using the same numerical method as for the annealing result, is visually indistinguishable from the exact solution. Therefore, the solution of the quantum annealer does not perfectly match the numerical solution of a classical computer. This is both due to errors in the quantum hardware as well as not enough iterations per timestep in our algorithm. We expect that the errors of future quantum annealers will be significantly lower.
Next, we analyze the convergence of our variational approach by performing a single Runge-Kutta step on a quantum annealer and compare it to numerical results using the same Runge-Kutta method. shows these errors for the Gauss-Legendre collocation method of order 6. The dashed line results are obtained on the D-Wave 2000Q system, while the solid lines are simulated on a digital computer by an exact solver, which finds the exact minimum of the combinatorial optimization function. We use two different exponent shift constants, given in Eq. (42), for the variational approach: The blue lines use c = 0.3, the green lines use c = 0.8. For the higher exponent shift, c = 0.8, first, the annealing solution follows the exact solution closely and converges faster than the lower convergence case. Nevertheless, after an error spike at the tenth iteration, the convergence stops and the error of the annealing solution is not further reduced.
For c = 0.3, the results of the quantum annealer follows the exact solver over all iterations, except of some spikes, which are caused by errors in the computing hardware. However, due to the lower exponent shift, the variational approach is still able to recover the same result even with errors present. Last, we solve a Ricatti equation, given by using our algorithm. This DEQ has two equilibrium points, an unstable one at u = 0 and a stable one at u = 1. As the integration method we used the Crank-Nicolson method with 2 qubits per register. For this case, only two auxiliary qubits are needed, making it 8 qubits in total, excluding those qubits which are needed for the embedding on a quantum annealer. Although the variational approach is not guaranteed to converge, we try the variational approach with c = 0.8 and 10 iterations. For this nonlinear case, we only present the results from an annealing simulator, not on the actual D-Wave 2000Q system. This is because of the non convex structure of the problem and high noise on the D-Wave system the results are hard to interpret. First, we tested the unstable equilibrium point and verified that the solution stays at the unstable equilibrium point, u = 0. Due to the high noise on current quantum annealers, this result is not achieved on a physical device. However, it shows that our algorithm is able to evaluate such solutions correctly. Next, we evaluate a dynamic case, starting at  u(0) = 0.1. In this case, the solution is expected to converge to the stable point u(∞) = 1. The results of the exact annealing are depicted in Fig. 11. Using the same integration method and time step, the annealed solution does not match the Newton method. We verified that the results by the Newton's algorithm are global minima of Eq. (39). The variational approach does not guarantee to converge to a global minima for non convex problems. Without the variational approach, with a sufficient amount of qubits for each register, and Eq. (39) being continuous and only having a single global minimum, we believe that our annealing algorithm is superior to the Newton's method. Namely, Newton's method is likely prone to local minima if the initial estimate for the method is not chosen correctly. A sufficient amount of qubits in our annealing algorithm ensures that a configuration close to the global minimum is a better minimizer to Eq. (39) compared to any other configuration close to a local minimum. However, this assumes hardware which is free of noise and a large amount of qubits.
Our findings serve as proof-of-concept and demonstrates the feasibility of our algorithm on current quantum annealers. The second example demonstrates that the exponent shift parameter has to be chosen wisely, dependent on the error of the quantum hardware. Finding best suitable exponent shifts for specific hardware is not covered in this paper and additional research is needed. The third example shows that with our algorithm the Runge-Kutta method can also be applied to nonlinear ODEs. However, for the nonlinear example studied, the annealed results do not match the expected results. A discussion about conditions, which we believe are necessary for nonlinear ODEs, is provided.

Discussion and Outlook
We have explored the usage of quantum computers for the task of solving differential equations, a ubiquitous problem in engineering and scientific computing. Our focus lies on the technique instead of a benchmark comparison, as the latter is not (yet) reasonable given the small problem size that had to studied here.
Concerning the time integration based on arithmetic operations on a digital quantum computer (Section 3), one might ask how this could possibly provide an advantage, given that the operations of a classical computer are directly mimicked. As scenario for exploiting the inherent parallelism of quantum computers, the circuits developed here could be used as building blocks within encompassing quantum algorithms, e.g., as an "oracle" within a quantum search method. Specifically, in the context of ordinary differential equations, this approach could solve inverse problems directly, like finding initial conditions given a desired state at the final time. From a broader perspective, we imagine that quantum computers could serve as parallel accelerators of subroutines within classical algorithms (like time integration).
Looking at the quantum annealing approach (Section 4), we see the largest potential for high-order implicit time integration methods, as the associated computational costs scales more favorably with integration order when using quantum annealing (compared to classical computers). In addition, we were able to decouple the number of required qubits from the accuracy of the number representation for the quantum annealing approach, and the arithmetic steps are implemented implicitly into the QUBO model, reducing the qubit count for temporary variables. Even though this decoupling is only guaranteed to work for linear ODEs, we applied it to a nonlinear ODE where our results did not match the expected results. We provided a discussion about the conditions, which we believe are necessary to solve nonlinear ODEs with our algorithm. We also gave an upper bound for the time complexity of our annealing algorithm if it is applied to an adiabatic quantum computer in the linear case. We showed that this complexity does not depend on the size of the problem itself, but rather on the discretization size of the variables and the time step size.
As an interesting alternative direction, we point out that the inherent dynamics of quantum systems could be used for solving nonlinear differential equations via the Koopman-von Neumann formulation of classical mechanics [19].