Solving Burgers’ equation with quantum computing

Computational fluid dynamics (CFD) simulations are a vital part of the design process in the aerospace industry. Although reliable CFD results can be obtained with turbulence models, direct numerical simulation of complex bodies in three spatial dimensions (3D) is impracticable due to the massive amount of computational elements. For instance, a 3D direct numerical simulation of a turbulent boundary-layer over the wing of a commercial jetliner that resolves all relevant length scales using a serial CFD solver on a modern digital computer would take approximately 750 million years or roughly 20% of the earth’s age. Over the past 25 years, quantum computers have become the object of great interest worldwide as powerful quantum algorithms have been constructed for several important, computationally challenging problems that provide enormous speed-up over the best-known classical algorithms. In this paper, we adapt a recently introduced quantum algorithm for partial differential equations to Burgers’ equation and develop a quantum CFD solver that determines its solutions. We used our quantum CFD solver to verify the quantum Burgers’ equation algorithm to find the flow solution when a shockwave is and is not present. The quantum simulation results were compared to: (i) an exact analytical solution for a flow without a shockwave; and (ii) the results of a classical CFD solver for flows with and without a shockwave. Excellent agreement was found in both cases, and the error of the quantum CFD solver was comparable to that of the classical CFD solver.


Introduction
Computational fluid dynamics (CFD) plays an essential role in the design of flight vehicles. It provides detailed, accurate predictions of the airflow about such vehicles, and it is a cost-effective alternative to expensive wind tunnel testing. However, highresolution CFD simulations of complex flows can take weeks or more to complete with existing supercomputers. Finding ways to speedup CFD simulations is a perennial problem.
Over the past 25 years, there has been growing interest in standing up a technology that will allow the construction of a robust, scalable quantum computer. The engineering challenge for quantum computer hardware is to reliably generate quantum entanglement and quantum state superposition in a scalable manner, while protecting these effects from decoherence, so that they can be exploited by quantum algorithms. If this can be done, quantum algorithms exist which, when run on a quantum computer, can speedup many important, computationally challenging problems. Perhaps the best-known example is Shor's quantum factoring algorithm [20] which provides an exponential speedup over the best existing classical factoring algorithms and, consequently, has greatly impacted the state of the art in cryptography.
It is natural to ask whether a quantum computer might allow a speedup of CFD simulations. Recently, a quantum algorithm for solving partial differential equations (PDE) has been introduced which was applied to the Navier-Stokes (NS) equations of fluid dynamics [8]. The resulting quantum algorithm was tested on a steady-state, inviscid, compressible nozzle flow problem which allows for shockwave formation and for which an exact solution is known. Excellent agreement with the exact solution was found, both when a shockwave was and was not present. The quantum NS algorithm was shown to provide a quadratic (exponential) speedup over classical random (deterministic) algorithms for non-smooth/turbulent flows.
Prior to Reference [8], Yepez used a quantum lattice-gas model to simulate the flow of a NS fluid [23][24][25], and to determine solutions of Burgers' equation [26]. The quantum lattice-gas model was shown to provide a significant advantage over a classical lattice-gas model. Steijl [21,22] used a hybrid quantum/classical approach to solve Poisson's equation which arises (inter alia) in incompressible NS flows. References [4,5,16,19] also examined fluid flow problems through quantum algorithms.
In this paper, we apply the quantum PDE algorithm of Reference [8] to Burgers' equation. We numerically simulate its application to flows in which a shockwave is and is not present. We find excellent agreement between our simulation results and: (i) an exact solution of a flow without a shockwave; and (ii) the results of a standard/classical CFD simulation of flows with and without a shockwave. The outline of this paper is as follows. In Sect. 2 we describe the application of the quantum PDE algorithm to Burgers' equation, while in Sect. 3 we present the results of our numerical simulation of the quantum Burgers' equation algorithm, and compare them to an exact solution and a classical CFD simulation as described above. Finally, we make closing remarks in Sect. 4.

Governing equations
Burgers' equation [2] (BE) is an important PDE that is often used in CFD as a simplified model for the Navier-Stokes dynamics. It contains both nonlinear and viscous terms and is widely used to test new techniques and algorithms. In one spatial dimension it is: Here u(x, t) is the local flow velocity and ν is the kinematic viscosity. It proves convenient to rewrite Eq. (1) in terms of non-dimensional variables: x * = x/L; u * = u/u ∞ ; and t * = tu ∞ /L, where L is a characteristic length and u ∞ is a characteristic velocity. The resulting equation is: where Re = u ∞ L/ν is the Reynolds number. In the remainder of this paper, we only work with non-dimensional variables and so will suppress the asterisk-superscript. For inviscid flow the RHS of Eq. (2) vanishes, and in conservation form it becomes:

Quantum BE algorithm
Following the prescription of Reference [8], we begin by discretizing space: The spatial boundary points correspond to the grid-points x 1 and x m . It is important to note that time t remains a continuous parameter. We also replace the spatial derivative by a first-order upwind scheme. This reduces the PDE, Eq.(3), to a coupled set of ordinary differential equations (ODEs): where, for simplicity, the lattice spacing x = x j − x j−1 is assumed to be constant.
We see that BE has been replaced by a coupled set of m − 2 ODEs. We postpone discussion of the initial-and boundary-conditions to Sect. 3.
To solve the BE system of ODEs (Eq. (5)) we use a quantum algorithm for solving systems of nonlinear ODEs introduced by Kacewicz [14]. Specifically, we look for a bounded function A( j, t) which approximates the exact solution u( j, t) over the time interval 0 ≤ t ≤ T . We require both u( j, t) and A( j, t) to satisfy the initial condition: The driver function f (u) is assumed to have continuous, bounded derivatives to order r , with the r th derivative satisfying the Hölder condition: Here H > 0 and 0 < ρ ≤ 1. The driver function's smoothness is parameterized by q = r + ρ, with q 1 (q ∼ 1) corresponding to smooth (non-smooth) functions. Functions satisfying these conditions are known as Hölder class functions [7,9] and are elements of the Hölder space F r ,ρ .
Kacewicz' quantum ODE algorithm begins by partitioning the time interval , and the approximate solution within T i,m by A i,m ( j, t). Taylor's method [1,11,17] is used to write A i,m ( j, t) as a truncated Taylor series about t i,m : For Hölder class functions f ∈ F r ,ρ , the parameter r is given. For a quasi-smooth driver function f (u), the parameter r is chosen so the error As noted earlier, we required that the {y i ( j)} provide the initial condition for the approximate solution A i ( j, t) for the i th primary subinterval T i . Thus, at t = t i ≡ t i,0 , we require: We see that once the parameters n, k, and are assigned values, the above construction determines the approximate solution A( j, t). How the parameters n and k are chosen is discussed in Sect. 3. We now explain how the {y i ( j)} are chosen.
To that end, Eq. (5) is integrated over T i : which is exact (the second term has been added and subtracted). To obtain an equation that relates the ; discards the third term on the RHS as it is O(h r +1 ); and writes τ =hz so that Eq. (8) becomes: for 0 ≤ i ≤ n − 2. Equation (9) determines y i+1 ( j) from y i ( j) and the Taylor poly- throughout the primary subinterval T 0 = [0, t 1 ] as described above. Equation (9) then determines y 1 ( j) from y 0 ( j), once the integral on the RHS is evaluated. To that end, Kacewicz introduces N k knot times {z m, p } in each secondary subinterval T i,m and approximates the integral by its average value over the knot times: The Quantum Amplitude Estimation Algorithm [3] (QAEA) is used to estimate the average value of f appearing on the RHS of Eq. (10). Note that this is the only task in Kacewicz' quantum algorithm that requires a quantum computer. Finally, before the QAEA can be used to evaluate Eq. (10), the summand f must be shifted and rescaled so that the new summand takes values in the range [0, 1]. Novak [18] and Heinrich [10] showed how the QAEA could be used to evaluate a function average, and Reference [8] explains how the shift and rescaling is implemented. In this way the {y 1 ( j)} are determined. They, in turn, determine the {A 1,m ( j, t)} throughout T 1 = [t 1 , t 2 ] as described above. This allows the RHS of Eq. (9) to be evaluated (using the QAEA to approximate the integral) giving the {y 2 ( j)}. Iterating this procedure over the remaining primary subintervals T i gives the approximate solution A( j, t), where A( j, t) = A i ( j, t) for t ∈ T i and 0 ≤ i ≤ n − 1. Reference [14] shows that for Hölder class functions the solution error satisfies (for n ≥ 5): Here α k = k(q + 1) − 1 and q = r + ρ is the driver function smoothness parameter. To reach this level of performance Kacewicz requires the upper bound 1 on the QAEA estimate of the function average and the probability 1 − δ 1 that this bound is satisfied [14] be given by 1 = 1/n k−1 and 1 − δ 1 = (1 − δ) 1/n k . In Sect. 3 we discuss how 1 and δ are assigned values.

Complexity analysis
Both quantum and classical algorithms for BE must discretize the spatial continuum, and in the interests of an apples-to-apples comparison, we assume both algorithms use the same discretization procedure. This guarantees the discretization costs for both are the same, and so comparison of the complexity of the quantum and classical BE algorithms reduces to the relative complexity of their respective ODE solvers. Numerical solution of an ODE requires multiple evaluations of its driver function f (u). We assume that an oracle is available. For a classical algorithm, the oracle is a black-box function/subroutine that returns f (u), while for a quantum algorithm, the oracle is an operator whose action encodes f(u) in the quantum state. The quantum oracle used in the QAEA is described in [8,18].
The computational cost of an algorithm A over a family of driver functions F is the maximum number of oracle calls needed by A to compute an approximate solution A( j, t) over all driver functions f ∈ F. It is thus the number of oracle calls needed to solve the hardest driver function in F. For a given > 0, the -complexity comp t (F, ) for an algorithm of type t = {deterministic, random, quantum} is the minimum (computational) cost over all algorithms A of type t whose error satisfies e t (A, F) < . It is thus the cost of the best algorithm of type t on driver functions f ∈ F.
Over a twenty year period, Kacewicz [12][13][14][15] determined upper and lower bounds on the -complexity for quantum, classical random, and classical deterministic ODE algorithms applied to Hölder class driver functions. We reproduce his results in Table 1.
The bounds depend on the error tolerance , and on the parameters (introduced in Sect. 2.1): (i) 0 < γ ≤ 1; (ii) 0 < ρ ≤ 1; and (iii) the smoothness parameter q = r + ρ, where r is the highest order derivative kept in the truncated Taylor series (see Eq. (7)). The upper bounds for the quantum and classical random ODE algorithms are for the algorithms introduced in [14]. The quoted lower bounds for quantum and classical random algorithms apply to all respective algorithms in these two classes, as shown in [12,13]. The algorithms in [14] are thus (almost) optimal as the exponents in the respective upper and lower bounds agree up to a small parameter γ , and in the quantum case, to within a logarithmic factor which we have suppressed. Kacewicz' classical deterministic ODE algorithm is seen to be optimal.
The -complexity is seen to be exponentially sensitive to the smoothness parameter q = r + ρ. For smooth driver functions, r is large (derivatives exist to high order) and so q 1. For non-smooth driver functions, r , q = O(1). In the most extreme case of a continuous, but non-differentiable function, r = 0 and so q = ρ ≤ 1.
For smooth functions (q 1), Table 1 gives Thus, for smooth driver functions, the complexity of all three types of ODE algorithm are the same slowly varying function of , and so there is no quantum speed-up in this case. However, for non-smooth driver functions, there is a quantum speed-up. This is apparent from the upper bounds in Table 1. We see that the exponent in the quantum upper bound is smaller than the corresponding exponent for the classical random and classical deterministic algorithms. Its complexity is thus smaller, corresponding to less oracle calls, and so to a shorter runtime. Thus for non-smooth driver functions there is a quantum speed-up. The degree of speed-up increases with decreasing q, being largest for q, γ 1. In this limit, The quantum algorithm thus has a square-root reduction in its -complexity (oracle calls) over classical random algorithms, and so a quadratic speed-up in runtime. The quantum speed-up is more pronounced when compared to classical deterministic algorithms where the speed-up is exponential in 1/q 1. To summarize: the quantum ODE algorithm provides a substantial quantum speedup over classical random and classical deterministic ODE algorithms for a non-smooth driver function, which is the computationally most challenging case. The largest speedup occurs in the regime q, γ 1. As noted at the beginning of this subsection, the same conclusion applies to the quantum BE algorithm.

Results
In this section, we present the results of a simulation of the quantum BE algorithm applied to two BE flows. In Sect. 3.1 we examine a smooth flow problem for which an exact analytical solution can be found, while in Sect. 3.2 we examine a flow which contains a travelling shockwave discontinuity. We also carried out a classical CFD simulation of the BE dynamics whose results we compare with the results of our quantum BE simulation.
Before entering into a presentation of our results, we specify the input parameters for our simulation. These parameter values were used in both simulations. We consider flows along the x-axis with 0 ≤ x ≤ 3. Recall that we are working with non-dimensional flow velocities, positions, and times. We discretize space by introducing a lattice with m = 61 grid-points. We chose the error bound 1 = 0.005, and the probability δ = 0.005 which is the probability that the quantum BE algorithm returns a solution that violates Eq. (11). We noted in Sect. 2.1 that Kacewicz requires 1 = 1/n k . Solving this expression for k gives: where z denotes the smallest integer greater than z. To obtain a second relation between n and k, we require that the durationh of a secondary subinterval be consistent with the Courant-Friedrichs-Lewy (CFL) stability condition. [6] This condition introduces a local time increment t( j) at each grid-point j, with C < 1. We chose C = 0.1, and the CFL time-increment is then t C F L = min {2≤ j≤m−1} t( j). The duration time of the simulation is then T = N t t C F L , where N t = 1000 was chosen. As discussed in Sect. 2.1, T =hn k . Setting these two expressions for T equal to each other (with minor rearrangement) gives: To insure our time partition is consistent with the CFL stability condition we require that n and k be chosen so thath/ t C F L < 1. This requires N t n k < 1.
Equations (14) and (17) allow n and k to be determined iteratively. Begin by choosing a value n 0 for n. Then use Eq. (14) to determine k 0 . Next evaluate N t /n k 0 0 and test whether it is less than 1. If yes, then set n = n 0 and k = k 0 and stop. Otherwise, set n 1 = n 0 + 1 and evaluate k 1 using Eq. (14). Then, test whether N t /n k 1 1 < 1. If yes, stop with n = n 1 and k = k 1 . Otherwise, continue to iterate in this fashion until first

Smooth flow
Here we look for a solution of BE for a smooth inviscid flow. The boundary condition at x 1 = 0 is u(x 1 , t) = 0. At x 61 = 3, because the flow is out of the simulation region, the boundary condition is a floating boundary condition. The flow velocity u(x 61 , t) is found by linear extrapolation from the flow velocity at x 59 and x 60 : u(x 61 , t) = 2u(x 60 , t) − u(x 59 , t). Finally, we impose the initial condition u(x, 0) = x, which satisfies both boundary conditions. An exact solution for this flow problem can be found using separation of variables: u(x, t) = X (x)T (t). Inserting this ansatz into BE, carrying out the separation, and imposing the initial condition gives X (x) = x and T (t) = 1/(1 + t). The exact analytic solution is then: We plot the exact solution in Fig. 1 at three different times: 0 ≤ t 1 < t 2 < t 3 ≤ T . The flow is seen to go asymptotically to zero everywhere with increasing time. To examine the performance of the quantum BE algorithm we wrote a quantum CFD solver (Q-CFD) to simulate the application of the quantum BE algorithm to this problem. The solver's output u q ( j, t) was compared with the exact analytical solution u a ( j, t). To assess the accuracy of the quantum BE algorithm solution we calculated the percent-error: at the end of the simulation: T = 1.667. We plot the percent-error in Fig. 2. We see that the quantum BE algorithm result is in excellent agreement with the exact solution. The largest percent-error occurs near x = 0 which is to be expected because of our use of an upwind scheme. The percent-error in that region is of the order of a few tenths of a percent. For x > 0.5, the percent-error drops to a few hundredths of a percent. Figure 2 also plots the percent-error for a classical CFD solver (C-CFD) applied to this problem. We see that that the quantum and classical solvers produce comparable errors.

Travelling shockwave
Here we examine the solution of BE when the initial condition contains a shockwave discontinuity: The boundary conditions are u(0, t) = 1 and u(3, t) = 0. These conditions require the initial shockwave to travel in the direction of increasing x. We ran the simulation until the shockwave discontinuity arrived at x = 3.0. This insures that the boundary condition at x = 3.0 is obeyed throughout the simulation. For this flow problem, no exact analytical solution is known. Thus we compared the results of our Q-CFD solver to that of a C-CFD solver. We plot the results for both solvers in Fig. 3 for t = 0, 0.83, 1.67, and 2.5. We see that: (i) the quantum BE algorithm successfully converged even though the flow contained a travelling shockwave discontinuity; and (ii) the Q-CFD results are in excellent agreement with that of the C-CFD solver. Finally, notice that both solvers slightly rounded the kinks in the flow solution by comparable amounts.

Conclusion
In this paper we have applied the quantum PDE algorithm of Reference [8] to BE in one-spatial dimension. We numerically simulated the quantum algorithm's application to inviscid flows with and without a shockwave. The results found were compared to an exact solution which is known when no shockwave is present in the flow, and to a standard CFD simulation when a shockwave is and is not present. Agreement was excellent, and the quantum algorithm's error was seen to be of the same order as that of the classical CFD simulation. We underscore that the quantum algorithm successfully converged to the solution even when a shockwave was present. We also discussed the computational complexity and speedup of the quantum BE algorithm which was found to be identical to that of the quantum PDE [8] and ODE [14] algorithms. It was shown that there is a quadratic (exponential) speedup over classical random (deterministic) algorithms for non-smooth driver functions which is the computationally most challenging case. It is important to note that this speedup is only expected to occur when the quantum algorithm is run on a quantum computer. Numerical simulation of the quantum algorithm on a classical computer is not expected to show a quantum speedup as classical computers do not generate the quantum entanglement or state superposition that underlie the speedup.
To close, we list a few useful directions for future work: -Application of the quantum PDE, NS, and BE algorithms to problems in two or more spatial dimensions; -Determining the quantum circuit implementation of the quantum PDE algorithm; and -Applying the quantum NS and BE algorithms to more complex flows.