Quantum-Solving Algorithm for d’Alembert Solutions of the Wave Equation

Abstract

When faced with a quantum-solving problem for partial differential equations, people usually transform such problems into Hamiltonian simulation problems or quantum-solving problems for linear equation systems. In this paper, we propose a third approach to solving partial differential equations that differs from the two approaches. By using the duality quantum algorithm, we construct a quantum-solving algorithm for solving the first-order wave equation, which represents a typical class of partial differential equations. Numerical results of the quantum circuit have high precision consistency with the theoretical d’Alembert solution. Then the routine is applied to the wave equation with either a dissipation or dispersion term. As shown by complexity analysis for all these cases of the wave equation, our algorithm has a quadratic acceleration for each iteration compared to the classical algorithm.

1. Introduction

Most scientific problems can be solved by studying the laws governing the evolution of physical quantities in space and time. Therefore, partial differential equations undoubtedly play an extremely important role in the field of natural sciences. However, the problem of solving partial differential equations is extremely difficult. While if quantum algorithms are introduced and the problems of partial differential equations are solved on a quantum computer, it can achieve accelerated characteristics compared to classical algorithms.

The usual quantum algorithm for solving partial differential equations proceeds as follows. First, discretize the space so that the function $f\left(x,t\right)$ becomes a vector $\mathbf{f}\left(t\right)$ and map its normalized components to the quantum state components, i.e., $|\mathbf{x}\left(t\right)\rangle ={\sum }_i{f}_i^{\prime }\left(t\right)|x_i\rangle$, where $f_i^{\prime }\left(t\right)$ is the i-th component of the vector $\mathbf{f}\left(t\right)$ after normalization. Next, the vectors encoded onto the quantum states are mapped onto a fixed model. Most quantum algorithms for solving partial differential equations rely on Hamiltonian simulations [1,2,3,4] or a linear equation system-solving algorithm (HHL algorithm) [5].

In the following, the main ideas of the above two solution methods will be briefly reviewed with examples. The solution method based on Hamiltonian simulation [6,7,8] that maps partial differential equations to the Schrödinger equation will be introduced first. This method maps equations with a similar structure to the Schrödinger equation to the Schrödinger equation and transforms the equation solving problem into a Hamiltonian simulation problem. For example, solving the Black–Scholes equation [9]

$$\frac{\partial f}{\partial t}=af+b\frac{\partial f}{\partial x}-c\frac{{\partial }^{2}f}{\partial x^2}.$$

(1)

The equation can be written in the following form

$$\frac{\partial f}{\partial t}=Af.$$

(2)

It is obvious that the equation is formally similar to the Schrödinger equation. Thus the A operator can be mapped to the Hamiltonian in Schrödinger’s equation in such a way that $A=ib\widehat{p}+\left(a\mathrm{I}+c{\widehat{p}}^2\right),\widehat{p}=-i{\partial }_x$. One can split A into Hermitian and anti-Hermitian parts, i.e., $A=A_H+A_{aH}$, where

$$A_{aH}=ib\widehat{p},A_H=a\mathrm{I}+c{\widehat{p}}^2.$$

(3)

The vector $\mathbf{f}\left(t\right)$, obtained by discretizing the function, is encoded onto the state vector $\left|\mathbf{x}\right(ϵ)\rangle$, using the Trotter product formula

$$|\mathbf{x}\left(ϵ\right)\rangle =e^{Aϵ}|{\mathbf{x}}_0\rangle \approx e^{A_Hϵ}{e}^{A_{aH}ϵ}|{\mathbf{x}}_0\rangle .$$

(4)

The problem of solving the partial differential equation is transformed into the problem of a Hamiltonian simulation. The process of simulating the action of the above Hamiltonian, i.e., the quantum state $|{\mathbf{x}}_0\rangle$, evolves under the designed Hamiltonian to obtain the final state. The solution of the original equation at different moments can be obtained by measuring the final state for different iterations.

In fact, it is efficient to use a Hamiltonian simulation to construct quantum algorithms for solving partial differential equations, which can solve first-order partial differential equations (requiring that the Hermitian and anti-Hermitian parts of the matrix A decomposition commute with each other) and second-order partial differential equations such as the wave equation. However, not all partial differential equations have the algebraic structure of Schrödinger’s equation. The quantum algorithm [7,10,11,12] for solving partial differential equations is presented below using the HHL algorithm. For the partial differential equation with the following structure after spatial discretization

$$\dot{\mathbf{x}}=A\mathbf{x}+\mathbf{b}.$$

(5)

Using Euler’s method to discretize time gives

$$\frac{\mathbf{x}\left(t_{j+1}\right)-\mathbf{x}\left(t_j\right)}{h}\approx A\mathbf{x}\left(t_j\right)+\mathbf{b}.$$

(6)

Let ${\mathbf{x}}_j=\mathbf{x}\left(t_j\right)$, the partial differential equation can be transformed into the following linear equation system; as an example, only the result of $j\le 2$ is given here,

$$\left[\begin{array}{ccc}\mathrm{I}& 0& 0\\ -(\mathrm{I}+Ah)& \mathrm{I}& 0\\ 0& -(\mathrm{I}+Ah)& \mathrm{I}\end{array}\right]\left[\begin{array}{c}{\mathbf{x}}_0\hfill \\ {\mathbf{x}}_1\hfill \\ {\mathbf{x}}_2\hfill \end{array}\right]=\left[\begin{array}{c}{\mathbf{x}}_{\mathrm{in}}\\ \mathbf{b}h\\ \mathbf{b}h\end{array}\right].$$

(7)

This system of equations is then solved using the HHL algorithm to obtain the following quantum states

$$|\mathbf{x}\rangle =\sum _{j=0}^{N_t}\left|t_j\right.〉\left|{\mathbf{x}}_j\right.〉.$$

(8)

The quantum state contains the solution of the partial differential equation $t_0$ to $t_j$ moments. Therefore, are there any other efficient algorithms for solving partial differential equations besides the above two methods? In this article, a third method different from the above two, the duality quantum algorithm [13,14,15,16] with amplitude amplification [17,18,19,20], is used to construct a quantum algorithm for solving the partial differential equation. The duality quantum algorithm also brings a speed-up effect compared to the classical algorithm [21,22,23].

This paper is organized as follows. First of all, the duality quantum algorithm will be used to solve the first-order wave equation with the d’Alembert solution. In the second and third parts of this paper, we will use the duality quantum algorithm to construct a solution algorithm to the wave equations with dissipation and dispersion terms. In these three parts of the paper, for these three problems, we use numerical simulations and present the results of the quantum algorithm solutions in the form of pictures for comparison with the theoretical values. At the end of this paper, we will analyze the complexity of our algorithm for solving the wave equations.

2. Duality Quantum Algorithm for Solving the First-Order Wave Equation

When talking about the wave equation, people must first think of the second-order linear hyperbolic type equation

$$\frac{{\partial }^{2}u}{\partial t^2}+k^2\frac{{\partial }^{2}u}{\partial x^2}=0.$$

(9)

The general solution can be written as $u\left(x,t\right)=f\left(x-kt\right)+g\left(x+kt\right)$, where f, g are two arbitrary functions. $f\left(x-kt\right)$ and $g\left(x+kt\right)$ represent waves passing along the x-axis with constant velocity to the right and to the left. Since Equation (9) is a linear homogeneous equations, its solutions are superposed. Therefore, f and g are two traveling waves that propagate independently without interfering with each other. If one focuses on only one of these two waves, Equation (9) degenerates to a linear hyperbolic equation of the first-order

$$\frac{\partial u}{\partial t}+k\frac{\partial u}{\partial x}=0.$$

(10)

In this paper, we discretize the continuous independent variable x into N points, i.e., $\mathbf{x}=\left(x_0,x_1,\cdots ,x_{N-1}\right)$. Then the spatial part of the function $u\left(x,t_j\right)$ at the moment $t_j$ is discretized into the vector

$$\mathbf{u}\left(\mathbf{x},t_j\right)=\left(u\left(x_0,t_j\right),u\left(x_1,t_j\right),\cdots ,u\left(x_{N-1},t_j\right)\right),$$

(11)

encode it onto the computational basis and define the quantum state ${|\psi \rangle }_j$ as

$${|\psi \rangle }_j=\frac{{\sum }_{i=0}^{N-1}u\left(x_i,t_j\right)|i\rangle }{\sqrt{{\sum }_{i=0}^{N-1}{u}^2\left(x_i,t_j\right)}}.$$

(12)

In the following, we will give the quantum algorithm for solving Equation (10) based on the non-unitary evolution of the quantum system. First, the Taylor expansion for each order partial differential term of Equation (10) is

$$\begin{array}{cc}\hfill & \frac{\partial u}{\partial t}=\frac{u\left(x,t+\tau \right)-u\left(x,t\right)}{\tau }+o\left(\tau \right),\hfill \\ \hfill & \frac{\partial u}{\partial x}=\frac{u\left(x+h,t\right)-u\left(x,t\right)}{h}+o\left(h\right).\hfill \end{array}$$

(13)

Pluging Equation (13) into Equation (10), the difference equation form of Equation (10) is obtained as

$$\frac{u\left(x_i,t_{j+1}\right)-u\left(x_i,t_j\right)}{\tau }+k\frac{u\left(x_{i+1},t_j\right)-u\left(x_i,t_j\right)}{h}=0.$$

(14)

Its local truncation error is $o(\tau +h)$. When $\tau$, $h\to 0$, Equation (14) approximates the original Equation (10). Organizing Equation (14) leads to

$$u\left(x_i,t_{j+1}\right)=\frac{\tau k\left[u\left(x_i,t_j\right)-u\left(x_{i+1},t_j\right)\right]}{h}+u\left(x_i,t_j\right)$$

(15)

Let $\mathsf{\Delta }=\frac{\tau }{h}$, then the following iterative relation can be obtained from Equation (15).

$$u\left(x_i,t_{j+1}\right)=\left(1+\mathsf{\Delta }k\right)u\left(x_i,t_j\right)-\mathsf{\Delta }ku\left(x_{i+1},t_j\right)$$

(16)

Taking the periodic boundary condition that $u\left(x_N,t\right)=\left(x_0,t\right)$, the equation describing the whole system can be written in the following form

$$\left[\begin{array}{c}u\left(x_0,t_{j+1}\right)\\ u\left(x_1,t_{j+1}\right)\\ \dots \\ u\left(x_{N-1},t_{j+1}\right)\end{array}\right]=A\left[\begin{array}{c}u\left(x_0,t_j\right)\\ u\left(x_1,t_j\right)\\ \\ \dots \\ u\left(x_{N-1},t_j\right)\end{array}\right]$$

(17)

where

$$A=\left[\begin{array}{ccccccc}1+\mathsf{\Delta }k& -\mathsf{\Delta }k& 0& & \dots & & 0\\ 0& 1+\mathsf{\Delta }k& -\mathsf{\Delta }k& 0& \dots & & 0\\ & \dots & & & & & \\ -\mathsf{\Delta }k& 0& & & \dots & 0& 1+\mathsf{\Delta }k\end{array}\right]$$

(18)

Then the state ${|\psi \rangle }_{j+1}$ of the system at the next moment, i.e., the moment $t_{j+1}$, can be expressed as ${A|\psi \rangle }_j$. It is obvious that the A-matrix is not an unitary matrix, so there is no way to achieve it directly by the product of quantum logic gates. Instead, the A-matrix has to be split into linear combinations of the unitary operators by the duality model of quantum computation, i.e., $A=(1+\mathsf{\Delta }k)A_0-\mathsf{\Delta }k{A}_1$, where $A_0$ is a unitary matrix of order N and

$$A_1={\left[\begin{array}{ccccc}0& 1& 0& & 0\\ 0& 0& 1& \dots & 0\\ & \cdots & & & \\ 0& 0& 0& & 1\\ 1& 0& 0& 0& 0\end{array}\right]}_{N\otimes N}.$$

(19)

By introducing an auxiliary qubit, the operation of the linear combination of unitary operators can be realized and thus equivalently the non-unitary evolution, i.e., ${A|\psi \rangle }_j$. Its quantum circuit is shown in Figure 1.

Where the matrix $A_1$ can be decomposed into ${\mathrm{C}}^n\left(X\right)$ gates as well as X gates with $O\left({log}_{2}N\right)$. The specific quantum circuit that implements the $A_1$ operation is shown in Figure 2.

According to Lemma 5.5 and Lemma 7.1 in article [24], a total of $4N-5{log}_{2}N-4$ CNOT gates and single-qubit rotating gates are needed if the $A_1$ operation continues to be disassembled. The following will explain the duality quantum algorithm for the solution of the first-order wave equation according to Figure 1, where first the auxiliary qubit passes through the $W_0$ gate, which has the following effect

$$W_0:|0\rangle \to \frac{\left(1+\mathsf{\Delta }k\right)|0\rangle -\mathsf{\Delta }k|1\rangle }{\sqrt{{\left(1+\mathsf{\Delta }k\right)}^2+{\left(\mathsf{\Delta }k\right)}^2}}.$$

(20)

Next, after two controlled quantum gates $|0\rangle \langle 0|\otimes A_0$ and $|1\rangle \langle 1|\otimes A_1$, the quantum state evolves as

$$\frac{\left(1+\mathsf{\Delta }k\right)|0\rangle A_0{|\psi \rangle }_j-\mathsf{\Delta }k|1\rangle A_1{|\psi \rangle }_j}{\sqrt{{\left(1+\mathsf{\Delta }k\right)}^2+{\left(\mathsf{\Delta }k\right)}^2}}.$$

(21)

Then the quantum state after the Hadamard transformation is

$$\frac{1}{\sqrt{2}}|0\rangle \left[\frac{\left(1+\mathsf{\Delta }k\right)A_0{|\psi \rangle }_j-\mathsf{\Delta }k{A}_1{|\psi \rangle }_j}{\sqrt{{\left(1+\mathsf{\Delta }k\right)}^2+{\left(\mathsf{\Delta }k\right)}^2}}\right]+\frac{1}{\sqrt{2}}|1\rangle \left[\frac{\left(1+\mathsf{\Delta }k\right)A_0{|\psi \rangle }_j+\mathsf{\Delta }k{A}_1{|\psi \rangle }_j}{\sqrt{{\left(1+\mathsf{\Delta }k\right)}^2+{\left(\mathsf{\Delta }k\right)}^2}}\right].$$

(22)

Finally, after the measurement to select the state of the auxiliary qubit as 0, the state of the working qubits at this time is ${|\psi \rangle }_{j+1}$, that is

$$\frac{\left(1+\mathsf{\Delta }k\right)A_0{|\psi \rangle }_j-\mathsf{\Delta }k{A}_1{|\psi \rangle }_j}{\sqrt{{\left(1+\mathsf{\Delta }k\right)}^2+{\left(\mathsf{\Delta }k\right)}^2}}.$$

(23)

Define the coefficients $C_j$ as

$$\sqrt{\sum _{i=0}^{N-1}{u}^2\left(x_i,t_j\right)\left[{\left(1+\mathsf{\Delta }k\right)}^2+{\left(\mathsf{\Delta }k\right)}^2\right]}.$$

(24)

The amplitude under the computational basis of the quantum state ${|\psi \rangle }_{j+1}$ is enlarged by a factor of $C_j$ to obtain the column vector $\mathbf{u}(x_i,t_{j+1})$, which is the state of the system at the moment $t_{j+1}$. The analysis yields that the computational complexity of this algorithm is $O\left(N\right)$ per iteration, while the complexity of the classical algorithm is $O\left(N^2\right)$. The specific calculation of the complexity is presented at the end of this paper.

The following equation will be used as an example

$$\left\{\begin{array}{c}\frac{\partial u}{\partial t}+2\frac{\partial u}{\partial x}=0\hfill \\ u\left(x,0\right)=-sin2\pi x+\frac{sin 4\pi x}{2}-\frac{sin 6\pi x}{3}\hfill \end{array}\right.$$

(25)

to show the duality quantum algorithm for the solution of the first-order wave equation. First of all, one period of the function, i.e., $[0,1]$, is chosen, and this interval is discretized into 32 points, i.e., $h=0.03125$ in Equation (15). Thus, the function value of 32 discrete points can be encoded using 5 qubits. Choose $\mathsf{\Delta }=0.1$ in Equation (16). Then $\tau =0.003125$, which represents the time interval for each iteration of the system evolution. According to Equation (20), the effect of the action of $W_0$ can be determined as

$$W_0:|0\rangle \to \frac{6\sqrt{37}}{37}|0\rangle -\frac{\sqrt{37}}{37}|1\rangle$$

(26)

This gives $W_0$ as $R_y(-0.33)$. Thus far, the quantum circuit for each iteration of the solution to Equation (25) can be given in Figure 3.

Numerical simulation of the first 10 iterations of this quantum circuit, whose results are shown in Figure 4. Where the orange curve represents the theoretical value. The blue points represent the results given by the numerical simulation of the quantum solution algorithm.

3. Duality Quantum Algorithm for the Solution of the Traveling Wave Dissipation Problem

By adding the dissipation term to Equation (10), the traveling wave equation with dissipation is obtained as

$$\frac{\partial u}{\partial t}+k\frac{\partial u}{\partial x}-\alpha \frac{{\partial }^{2}u}{\partial x^2}=0.$$

(27)

The Taylor expansion for each term of Equation (27) is

$$\begin{array}{cc}\hfill \frac{\partial u}{\partial t}& =\frac{u\left(x,t+\tau \right)-u\left(x,t\right)}{\tau }+o\left(\tau \right),\hfill \\ \hfill \frac{\partial u}{\partial x}& =\frac{u\left(x+h,t\right)-u\left(x,t\right)}{h}+o\left(h\right),\hfill \\ \hfill \frac{{\partial }^{2}u}{\partial x^2}& =\frac{u(x-h,t)-2u(x,t)+u(x+h,t)}{h^2}+o\left(h^2\right).\hfill \end{array}$$

(28)

Differentiating Equation (27) yields

$$\begin{array}{cc}\hfill \frac{u\left(x_i,t_{j+1}\right)-u\left(x_i,t_j\right)}{\tau }+& k\frac{u\left(x_{i+1},t_j\right)-u\left(x_i,t_j\right)}{h}\hfill \\ \hfill -& \alpha \frac{u\left(x_{i-1},t_j\right)-2u\left(x_i,t_j\right)+u\left(x_{i+1},t_j\right)}{h^2}=0.\hfill \end{array}$$

(29)

The collation leads to

$$\begin{array}{cc}\hfill u\left(x_i,t_{j+1}\right)=& \tau k\frac{u\left(x_i,t_j\right)-u\left(x_{i+1},t_j\right)}{h}\hfill \\ \hfill & +\tau \alpha \frac{u\left(x_{i-1},t_j\right)-2u\left(x_i,t_j\right)+u\left(x_{i+1},t_j\right)}{h^2}+u\left(x_i,t_j\right).\hfill \end{array}$$

(30)

Let $\mathsf{\Delta }=\frac{\tau }{h}$ and ${\mathsf{\Delta }}_1=\alpha \frac{\mathsf{\Delta }}{h}$, then the following iterative relation can be obtained

$$u\left(x_i,t_{j+1}\right)={\mathsf{\Delta }}_{1}u\left(x_{i-1},t_j\right)+\left(k\mathsf{\Delta }-2{\mathsf{\Delta }}_1+1\right)u\left(x_i,t_j\right)+\left({\mathsf{\Delta }}_1-k\mathsf{\Delta }\right)u\left(x_{i+1},t_j\right).$$

(31)

Take the periodic boundary condition that $u\left(x_N,t\right)=u\left(x_0,t\right)$. Then, the equation describing the whole system can be written in the following form

$$\left[\begin{array}{c}u\left(x_0,t_{j+1}\right)\\ u\left(x_1,t_{j+1}\right)\\ \dots \\ u\left(x_{N-1},t_{j+1}\right)\end{array}\right]=A\left[\begin{array}{c}u\left(x_0,t_j\right)\\ u\left(x_1,t_j\right)\\ \dots \\ u\left(x_{N-1},t_j\right)\end{array}\right].$$

(32)

In Equation (32)

$$A=\left[\begin{array}{ccccccccc}a& b& 0& & \dots & & & 0& c\\ c& a& b& 0& \dots & & & & 0\\ & \dots & & & & & & \\ 0& & & & \dots & 0& c& a& b\\ b& 0& & & \dots & & 0& c& a\end{array}\right],$$

(33)

in which

$$\begin{array}{cc}\hfill & a=k\mathsf{\Delta }-2{\mathsf{\Delta }}_1+1,\hfill \\ \hfill & b={\mathsf{\Delta }}_1-k\mathsf{\Delta },\hfill \\ \hfill & c={\mathsf{\Delta }}_1.\hfill \end{array}$$

(34)

Following the encoding method of Equation (12), then the state ${|\psi \rangle }_{j+1}$ of the system at the next moment, i.e., the moment $t_{j+1}$, can be expressed as ${A|\psi \rangle }_j$. The matrix $A=\left(k\mathsf{\Delta }-2{\mathsf{\Delta }}_1+1\right)A_0+\left({\mathsf{\Delta }}_1-k\mathsf{\Delta }\right)A_1+{\mathsf{\Delta }}_1{A}_2$, where $A_2=A_1^{†}$. Thus the operation of a linear combination of unitary operators can be equivalently implemented by introducing two auxiliary qubits, whose quantum circuit is shown in Figure 5.

In the following, the duality quantum algorithm for solving the dissipation problem of the first-order wave equation is explained in conjunction with the quantum circuit (Figure 5), where the first the auxiliary qubits pass through the $W_1$ gate, which has the following effect

$$W_1:|00\rangle \to \frac{\left(k\mathsf{\Delta }-2{\mathsf{\Delta }}_1+1\right)|00\rangle +\left({\mathsf{\Delta }}_1-k\mathsf{\Delta }\right)|01\rangle +{\mathsf{\Delta }}_1|10\rangle }{\sqrt{{\left(k\mathsf{\Delta }-2{\mathsf{\Delta }}_1+1\right)}^2+{\left({\mathsf{\Delta }}_1-k\mathsf{\Delta }\right)}^2+{\mathsf{\Delta }}_1^2}}.$$

(35)

Next, after three controlled quantum gates $|00\rangle \langle 00|\otimes A_0$, $|01\rangle \langle 01|\otimes A_1$ and $|10\rangle \langle 10|\otimes A_2$ the quantum state evolves as

$$\frac{\left(k\mathsf{\Delta }-2{\mathsf{\Delta }}_1+1\right)|00\rangle A_0{|\psi \rangle }_j+\left({\mathsf{\Delta }}_1-k\mathsf{\Delta }\right)|01\rangle A_1{|\psi \rangle }_j+{\mathsf{\Delta }}_1|10\rangle A_2{|\psi \rangle }_j}{\sqrt{{\left(k\mathsf{\Delta }-2{\mathsf{\Delta }}_1+1\right)}^2+{\left({\mathsf{\Delta }}_1-k\mathsf{\Delta }\right)}^2+{\mathsf{\Delta }}_1^2}}.$$

(36)

Then, the quantum state evolves after the Hadamard transformation of two auxiliary qubits as

$$\begin{array}{cc}\hfill \frac{1}{2}& |00\rangle \left[\frac{\left(k\mathsf{\Delta }-2{\mathsf{\Delta }}_1+1\right)A_0{|\psi \rangle }_j+\left({\mathsf{\Delta }}_1-k\mathsf{\Delta }\right)A_1{|\psi \rangle }_j+{\mathsf{\Delta }}_1{A}_2{|\psi \rangle }_j}{\sqrt{{\left(k\mathsf{\Delta }-2{\mathsf{\Delta }}_1+1\right)}^2+{\left({\mathsf{\Delta }}_1-k\mathsf{\Delta }\right)}^2+{\mathsf{\Delta }}_1^2}}\right]\hfill \\ \hfill & +\frac{1}{2}|01\rangle \left[\frac{\left(k\mathsf{\Delta }-2{\mathsf{\Delta }}_1+1\right)A_0{|\psi \rangle }_j-\left({\mathsf{\Delta }}_1-k\mathsf{\Delta }\right)A_1{|\psi \rangle }_j+{\mathsf{\Delta }}_1{A}_2{|\psi \rangle }_j}{\sqrt{{\left(k\mathsf{\Delta }-2{\mathsf{\Delta }}_1+1\right)}^2+{\left({\mathsf{\Delta }}_1-k\mathsf{\Delta }\right)}^2+{\mathsf{\Delta }}_1^2}}\right]\hfill \\ \hfill & +\frac{1}{2}|10\rangle \left[\frac{\left(k\mathsf{\Delta }-2{\mathsf{\Delta }}_1+1\right)A_0{|\psi \rangle }_j+\left({\mathsf{\Delta }}_1-k\mathsf{\Delta }\right)A_1{|\psi \rangle }_j-{\mathsf{\Delta }}_1{A}_2{|\psi \rangle }_j}{\sqrt{{\left(k\mathsf{\Delta }-2{\mathsf{\Delta }}_1+1\right)}^2+{\left({\mathsf{\Delta }}_1-k\mathsf{\Delta }\right)}^2+{\mathsf{\Delta }}_1^2}}\right]\hfill \\ \hfill & +\frac{1}{2}|11\rangle \left[\frac{\left(k\mathsf{\Delta }-2{\mathsf{\Delta }}_1+1\right)A_0{|\psi \rangle }_j-\left({\mathsf{\Delta }}_1-k\mathsf{\Delta }\right)A_1{|\psi \rangle }_j-{\mathsf{\Delta }}_1{A}_2{|\psi \rangle }_j}{\sqrt{{\left(k\mathsf{\Delta }-2{\mathsf{\Delta }}_1+1\right)}^2+{\left({\mathsf{\Delta }}_1-k\mathsf{\Delta }\right)}^2+{\mathsf{\Delta }}_1^2}}\right].\hfill \end{array}$$

(37)

Finally, the auxiliary qubits are measured to select the state of the auxiliary qubit as 00; then the state of the working qubit at this time is ${|\psi \rangle }_{j+1}$, that is

$$\frac{\left(k\mathsf{\Delta }-2{\mathsf{\Delta }}_1+1\right)A_0{|\psi \rangle }_j+\left({\mathsf{\Delta }}_1-k\mathsf{\Delta }\right)A_1{|\psi \rangle }_j+{\mathsf{\Delta }}_1{A}_2{|\psi \rangle }_j}{\sqrt{{\left(k\mathsf{\Delta }-2{\mathsf{\Delta }}_1+1\right)}^2+{\left({\mathsf{\Delta }}_1-k\mathsf{\Delta }\right)}^2+{\mathsf{\Delta }}_1^2}}.$$

(38)

Define the coefficient $C_j$ as

$$\sqrt{\sum _{i=0}^{N-1}{u}^2\left(x_i,t_j\right)\left[{\left(k\mathsf{\Delta }-2{\mathsf{\Delta }}_1+1\right)}^2+{\left({\mathsf{\Delta }}_1-k\mathsf{\Delta }\right)}^2+{\mathsf{\Delta }}_1^2\right]}.$$

(39)

The column vector $\mathbf{u}(x_i,t_{j+1})$, which is the state of the system at $t_{j+1}$ moments, is obtained by amplifying the amplitude under the computational basis of the quantum state ${|\psi \rangle }_{j+1}$ by a factor of $C_j$.

According to Lemma 5.5 and Lemma 7.1 in article [24] and combined with Figure 5, the computational complexity of this algorithm per iteration is $O\left(N\right)$, while the complexity of the classical algorithm is $O\left(N^2\right)$. Thus, the present algorithm has the property of speeding up in each iteration compared to the classical algorithm. The specific calculation of the complexity is presented at the end of this paper.

The following equation will be used as an example

$$\{\begin{array}{c}\frac{\partial u}{\partial t}+2\frac{\partial u}{\partial x}-0.1\frac{{\partial }^{2}u}{\partial x^2}=0\hfill \\ u(x,0)=-sin 2\pi x+\frac{sin 4\pi x}{2}-\frac{sin 6\pi x}{3}\hfill \end{array}$$

(40)

to show the duality quantum algorithm for the solution of the dissipation problem of the first-order wave equation. First, one period of the function is chosen, i.e., $[0,1]$, and this interval is discretized into 32 points, i.e., $h=0.03125$ in Equation (28). Thus, the function value of 32 discrete points can be encoded using 5 qubits. Choose $\mathsf{\Delta }=0.2$ and ${\mathsf{\Delta }}_1=1.28$ in Equation (31). Then $\tau =0.00625$, which represents the time interval for each iteration of system evolution. According to Equation (35), the effect of the action of $W_1$ can be determined as

$$W_1:|00\rangle \to -\frac{\sqrt{29}}{9}|00\rangle +\frac{22\sqrt{29}}{261}|01\rangle +\frac{32\sqrt{29}}{261}|10\rangle .$$

(41)

It is constructed as shown in Figure A1 with Equations (A1) and (A2). The revolving gate $R_n\left(\theta \right)$ of the first auxiliary qubit is constructed according to Equation (A1), such that

$$R_n\left(\theta \right)|00\rangle \to \left(\frac{5\sqrt{1537}}{261}|0\rangle +\frac{32\sqrt{29}}{261}|1\rangle \right)|0\rangle .$$

(42)

It is obtained that $R_n\left(\theta \right)$ is $R_y\left(1.442\right)$. According to Equation (A3), the controlled operator $U_1$ is to achieve the following action

$$U_1|0\rangle \to \frac{29\sqrt{53}}{265}|0\rangle +\frac{22\sqrt{53}}{265}|1\rangle .$$

(43)

The controlled operator $U_1$ can be obtained as $R_y\left(1.298\right)$. At this point, we can give the quantum circuit for each iteration of the solution Equation (40), as shown in Figure 6

Numerical simulation of the first 10 iterations of this quantum circuit results in Figure 7. The orange curve represents the resolved theoretical value. The blue points represent the values solved by the numerical simulation of the quantum algorithm.

4. Duality Quantum Algorithm for Solving Traveling Wave Dispersion Problems

By adding the dispersion term to Equation (10), the traveling wave equation with dispersion is obtained as

$$\frac{\partial u}{\partial t}+k\frac{\partial u}{\partial x}+\beta \frac{{\partial }^{3}u}{\partial x^3}=0.$$

(44)

Taylor expansion of the terms of Equation (44)

$$\begin{array}{cc}\hfill & \frac{\partial u}{\partial t}=\frac{u\left(x,t+\tau \right)-u\left(x,t\right)}{\tau }+o\left(\tau \right),\hfill \\ \hfill & \frac{\partial u}{\partial x}=\frac{u\left(x+h,t\right)-u\left(x,t\right)}{h}+o\left(h\right),\hfill \\ \hfill & \frac{{\partial }^{3}u}{\partial x^3}=\frac{-u(x-2h,t)+3u(x-h,t)-3u(x,t)+u(x+h,t)}{h^3}+o\left(h^3\right).\hfill \end{array}$$

(45)

Differentiating Equation (44) yields

$$\begin{array}{cc}\hfill & \frac{u\left(x_i,t_{j+1}\right)-u\left(x_i,t_j\right)}{\tau }+k\frac{u\left(x_{i+1},t_j\right)-u\left(x_i,t_j\right)}{h}\hfill \\ \hfill & +\beta \frac{-u\left(x_{i-2},t_j\right)+3u\left(x_{i-1},t_j\right)-3u\left(x_i,t_j\right)+u\left(x_{i+1},t_j\right)}{h^3}=0.\hfill \end{array}$$

(46)

Let $\mathsf{\Delta }=\frac{\tau }{h}$, ${\mathsf{\Delta }}_2=\beta \frac{\mathsf{\Delta }}{h^2}$, and the following iterative relation can be obtained

$$\begin{array}{cc}\hfill u\left(x_i,t_{j+1}\right)={\mathsf{\Delta }}_{2}u\left(x_{i-2},t_j\right)-& 3{\mathsf{\Delta }}_{2}u\left(x_{i-1},t_j\right)\hfill \\ \hfill & +\left(k\mathsf{\Delta }+3{\mathsf{\Delta }}_2+1\right)u\left(x_i,t_j\right)-\left(k\mathsf{\Delta }+{\mathsf{\Delta }}_2\right)u\left(x_{i+1},t_j\right).\hfill \end{array}$$

(47)

Take the periodic boundary condition that $u\left(x_N,t\right)=u\left(x_0,t\right)$. Then the equation describing the whole system can be written in the following form

$$\left[\begin{array}{c}u\left(x_0,t_{j+1}\right)\\ u\left(x_1,t_{j+1}\right)\\ \dots \\ u\left(x_{N-1},t_{j+1}\right)\end{array}\right]=A\left[\begin{array}{c}u\left(x_0,t_j\right)\\ u\left(x_1,t_j\right)\\ \dots \\ u\left(x_{N-1},t_j\right)\end{array}\right],$$

(48)

where

$$A=\left[\begin{array}{ccccccccccccc}a& b& 0& & & & \dots & & \dots & & 0& d& c\\ c& a& b& 0& & & \dots & & \dots & & & 0& d\\ d& c& a& b& 0& & \dots & & \dots & & & & 0\\ 0& d& c& a& b& 0& \dots & & \dots & & & & 0\\ & & \dots & & & & \dots & & \dots & & & & \\ 0& & & \dots & & & \dots & 0& d& c& a& b& 0\\ 0& & & \dots & & & \dots & & 0& d& c& a& b\\ b& 0& & \dots & & & \dots & & & 0& d& c& a\end{array}\right],$$

(49)

in which

$$\begin{array}{cc}\hfill & a=k\mathsf{\Delta }+3{\mathsf{\Delta }}_2+1,\hfill \\ \hfill & b=-\left(k\mathsf{\Delta }+{\mathsf{\Delta }}_2\right),\hfill \\ \hfill & c=-3{\mathsf{\Delta }}_2,\hfill \\ \hfill & d={\mathsf{\Delta }}_2.\hfill \end{array}$$

(50)

It can be seen that $A=\left(k\mathsf{\Delta }+3{\mathsf{\Delta }}_2+1\right)A_0-\left(k\mathsf{\Delta }+{\mathsf{\Delta }}_2\right)A_1-3{\mathsf{\Delta }}_2{A}_2+{\mathsf{\Delta }}_2{A}_3$, where

$$A_3={\left[\begin{array}{ccccccccc}0& 0& & & & 0& 0& 1& 0\\ 0& 0& & & & 0& 0& 0& 1\\ 1& 0& & & \dots & 0& 0& 0& 0\\ 0& 1& & & & 0& 0& 0& 0\\ & & \dots & & & & & & \\ & \dots & & 1& 0& 0& 0& 0& 0\\ & & & 0& 1& 0& 0& 0& 0\\ 0& 0& & 0& 0& 1& 0& 0& 0\\ 0& 0& & 0& 0& 0& 1& 0& 0\end{array}\right]}_{N\otimes N}.$$

(51)

Therefore, the operation of the linear combination of unitary operators can be equivalently implemented by introducing two auxiliary qubits, whose quantum circuit is shown in Figure 8.

It is not difficult to find $A_3=A_2^2$, and the following will be combined with the quantum circuit (Figure 8) to explain the duality quantum algorithm for the solution of the dispersion problem of the first-order wave equation. First, the auxiliary qubits pass through the $W_2$ gate, the effect of which is as follows

$$W_2:|00\rangle \to \frac{\left(k\mathsf{\Delta }+3{\mathsf{\Delta }}_2+1\right)|00\rangle +\left(k\mathsf{\Delta }+{\mathsf{\Delta }}_2\right)|01\rangle -3{\mathsf{\Delta }}_2|10\rangle +{\mathsf{\Delta }}_2|11\rangle }{\sqrt{{\left(k\mathsf{\Delta }+3{\mathsf{\Delta }}_2+1\right)}^2+{\left(k\mathsf{\Delta }+{\mathsf{\Delta }}_2\right)}^2+10{\mathsf{\Delta }}_2^2}}.$$

(52)

Next, after four controlled quantum gates, $|00\rangle \langle 00|\otimes A_0$, $|01\rangle \langle 01|\otimes A_1$, $|10\rangle \langle 10|\otimes A_2$ and $|11\rangle \langle 11|\otimes A_3$, the quantum state evolves as

$$\frac{\left(k\mathsf{\Delta }+3{\mathsf{\Delta }}_2+1\right)|00\rangle A_0{|\psi \rangle }_j+\left(k\mathsf{\Delta }+{\mathsf{\Delta }}_2\right)|01\rangle A_1{|\psi \rangle }_j-3{\mathsf{\Delta }}_2|10\rangle A_2{|\psi \rangle }_j+{\mathsf{\Delta }}_2|11\rangle A_3{|\psi \rangle }_j}{\sqrt{{\left(k\mathsf{\Delta }+3{\mathsf{\Delta }}_2+1\right)}^2+{\left(k\mathsf{\Delta }+{\mathsf{\Delta }}_2\right)}^2+10{\mathsf{\Delta }}_2^2}}.$$

(53)

Then, the Hadamard transform is performed for the two auxiliary qubits, and the quantum state evolves as

$$\begin{array}{cc}\hfill \frac{1}{2}& |00\rangle \left[\frac{\left(k\mathsf{\Delta }+3{\mathsf{\Delta }}_2+1\right)A_0{|\psi \rangle }_j+\left(k\mathsf{\Delta }+{\mathsf{\Delta }}_2\right)A_1{|\psi \rangle }_j-3{\mathsf{\Delta }}_2{A}_2{|\psi \rangle }_j+{\mathsf{\Delta }}_2{A}_3{|\psi \rangle }_j}{\sqrt{{\left(k\mathsf{\Delta }+3{\mathsf{\Delta }}_2+1\right)}^2+{\left(k\mathsf{\Delta }+{\mathsf{\Delta }}_2\right)}^2+10{\mathsf{\Delta }}_2^2}}\right]\hfill \\ \hfill +& \frac{1}{2}|01\rangle \left[\frac{\left(k\mathsf{\Delta }+3{\mathsf{\Delta }}_2+1\right)A_0{|\psi \rangle }_j-\left(k\mathsf{\Delta }+{\mathsf{\Delta }}_2\right)A_1{|\psi \rangle }_j-3{\mathsf{\Delta }}_2{A}_2{|\psi \rangle }_j-{\mathsf{\Delta }}_2{A}_3{|\psi \rangle }_j}{\sqrt{{\left(k\mathsf{\Delta }+3{\mathsf{\Delta }}_2+1\right)}^2+{\left(k\mathsf{\Delta }+{\mathsf{\Delta }}_2\right)}^2+10{\mathsf{\Delta }}_2^2}}\right]\hfill \\ \hfill +& \frac{1}{2}|10\rangle \left[\frac{\left(k\mathsf{\Delta }+3{\mathsf{\Delta }}_2+1\right)A_0{|\psi \rangle }_j+\left(k\mathsf{\Delta }+{\mathsf{\Delta }}_2\right)A_1{|\psi \rangle }_j+3{\mathsf{\Delta }}_2{A}_2{|\psi \rangle }_j-{\mathsf{\Delta }}_2{A}_3{|\psi \rangle }_j}{\sqrt{{\left(k\mathsf{\Delta }+3{\mathsf{\Delta }}_2+1\right)}^2+{\left(k\mathsf{\Delta }+{\mathsf{\Delta }}_2\right)}^2+10{\mathsf{\Delta }}_2^2}}\right]\hfill \\ \hfill +& \frac{1}{2}|11\rangle \left[\frac{\left(k\mathsf{\Delta }+3{\mathsf{\Delta }}_2+1\right)A_0{|\psi \rangle }_j-\left(k\mathsf{\Delta }+{\mathsf{\Delta }}_2\right)A_1{|\psi \rangle }_j+3{\mathsf{\Delta }}_2{A}_2{|\psi \rangle }_j+{\mathsf{\Delta }}_2{A}_3{|\psi \rangle }_j}{\sqrt{{\left(k\mathsf{\Delta }+3{\mathsf{\Delta }}_2+1\right)}^2+{\left(k\mathsf{\Delta }+{\mathsf{\Delta }}_2\right)}^2+10{\mathsf{\Delta }}_2^2}}\right].\hfill \end{array}$$

(54)

Finally, the state of the auxiliary qubit is selected as 00 after measurement, then the state of the working qubits at this time is ${|\psi \rangle }_{j+1}$, that is

$$\frac{\left(k\mathsf{\Delta }+3{\mathsf{\Delta }}_2+1\right)A_0{|\psi \rangle }_j+\left(k\mathsf{\Delta }+{\mathsf{\Delta }}_2\right)A_1{|\psi \rangle }_j-3{\mathsf{\Delta }}_2{A}_2{|\psi \rangle }_j+{\mathsf{\Delta }}_2{A}_3{|\psi \rangle }_j}{\sqrt{{\left(k\mathsf{\Delta }+3{\mathsf{\Delta }}_2+1\right)}^2+{\left(k\mathsf{\Delta }+{\mathsf{\Delta }}_2\right)}^2+10{\mathsf{\Delta }}_2^2}}.$$

(55)

Define the coefficient $C_j$ as

$$\sqrt{\sum _{i=0}^{N-1}{u}^2\left(x_i,t_j\right)\left[{\left(k\mathsf{\Delta }+3{\mathsf{\Delta }}_2+1\right)}^2+{\left(k\mathsf{\Delta }+{\mathsf{\Delta }}_2\right)}^2+10{\mathsf{\Delta }}_2^2\right]}.$$

(56)

The amplitude under the computational basis of the quantum state ${|\psi \rangle }_{j+1}$ is enlarged $C_j$ times to obtain the column vector $\mathbf{u}\left(x_i,t_{j+1}\right)$, which is the state of the system at the moment $t_{j+1}$. According to Lemma 5.5 and Lemma 7.1 in article [24], and combined with the analysis of Figure 8, we can obtain that the computational complexity of this algorithm for each iteration is $O\left(N\right)$, while the complexity of the classical algorithm is $O\left(N^2\right)$. The specific calculation of the complexity is presented at the end of this paper.

The following equation will be used as an example

$$\left\{\begin{array}{c}\frac{\partial u}{\partial t}+2\frac{\partial u}{\partial x}+0.01\frac{{\partial }^{3}u}{\partial x^3}=0\hfill \\ u\left(x,0\right)=-sin2\pi x+\frac{sin 4\pi x}{2}-\frac{sin 6\pi x}{3}\hfill \end{array}\right.$$

(57)

to show the duality quantum algorithm for the solution of the dispersion problem of the first-order wave equation. Firstly, one period of the function is chosen, i.e., $[0,1]$, and this interval is discretized into 32 points, i.e., $h=0.03125$ in Equation (47). Thus, the function value of 32 discrete points can be encoded using 5 qubits. Choose $\mathsf{\Delta }=0.05$ and ${\mathsf{\Delta }}_2=0.512$ in Equation (47). Then $\tau =0.0015625$, which represents the time interval for each iteration of the system evolution. According to Equation (52), the effect of the action of $W_2$ can be determined as

$$W_2:|00\rangle \to 0.8359|00\rangle +0.1941|01\rangle -0.4871|10\rangle +0.1624|11\rangle .$$

(58)

It is constructed as shown in Figure A1 with Equations (A1) and (A2). Construct the revolving gate $R_n\left(\theta \right)$ of the first auxiliary qubit according to Equation (A1), such that

$$R_n\left(\theta \right)|00\rangle \to \left(0.8581|0\rangle +0.5135|1\rangle \right)|0\rangle .$$

(59)

It is obtained that $R_n\left(\theta \right)$ is $R_y\left(1.078\right)$.

According to Equation (A3), two controlled operators $U_1$, $U_2$ are to be realized as follows

$$\begin{array}{cc}\hfill U_1& |0\rangle \to 0.9741|0\rangle +0.2262|1\rangle ,\hfill \\ \hfill U_2& |0\rangle \to -0.9486|0\rangle +0.3163|1\rangle .\hfill \end{array}$$

(60)

The controlled operator $U_1$ can be obtained as $R_y\left(0.4563\right)$ and $U_2$ as $R_y\left(5.639\right)$. The quantum circuit for each iteration of the solution Equation (57) so far is given in Figure 9.

And the specific quantum circuit of the R-operation in Figure 9 is shown in Figure 10.

The result of numerically simulating the first 10 iterations of this quantum circuit is shown in Figure 11. The orange curve represents the theoretical value. The blue points represent the values solved by the numerical simulation quantum algorithm.

5. Discussion

For a quantum algorithm that solves a d-dimensional partial differential equation (meaning that there are d spatial variables) with a spatial discretization number of N, the output is an approximation $C\left(f\right)$ of the function f with an error $ϵ$. In fact, for the problem of quantum algorithms solving partial differential equations, the number of discrete points N and the error $ϵ$ are interrelated [12,25]. The correlations are as follows

$$N=O\left(\mathrm{poly}\left(\frac{1}{{ϵ}^d}\right)\right).$$

(61)

For the preparation of the initial state, its complexity is $O\left(\mathrm{poly}log\left(N\right)\right)$. The complexity of each iteration of the algorithm in this paper will be given below. First, according to Lemma 5.5 and Lemma 7.1 in article [24], it can be obtained that the controlled gate ${\mathrm{C}}^{n-1}\left(U\right)$ for n qubits having $n-1$ control qubits can be split into CNOT gates and single-qubit gates for a total of $2^{n+1}-5$, where $n\ge 3$. For the quantum circuit in Figure 1, the total number of elementary quantum gates required is

$$\begin{array}{c}\hfill 3+O\left({\mathrm{C}}^2\left(U\right)\right)+\cdots +O\left({\mathrm{C}}^n\left(U\right)\right)\\ \hfill =2^{n+3}-5n-8=8N-5{log}_{2}N-8\simeq O\left(N\right)\end{array}$$

(62)

For the quantum circuit in Figure 5, the total number of elementary quantum gates required is

$$\begin{array}{c}\hfill 8+2\left[O\left({\mathrm{C}}^2\left(U\right)\right)+\cdots +O\left({\mathrm{C}}^{n+1}\left(U\right)\right)\right]\\ \hfill =2^{n+5}-10n-24=32N-10{log}_{2}N-24\simeq O\left(N\right)\end{array}$$

(63)

For the quantum circuit in Figure 8, the total number of elementary quantum gates required is

$$\begin{array}{c}\hfill 8+4\left[O\left({\mathrm{C}}^2\left(U\right)\right)+\cdots +O\left({\mathrm{C}}^{n+1}\left(U\right)\right)\right]\\ \hfill =2^{n+6}-20n-56=64N-20{log}_{2}N-56\simeq O\left(N\right)\end{array}$$

(64)

It can be found that the quantum-solving algorithm given in this paper has a quadratic acceleration for each iteration compared to the classical algorithm. However, the state of the auxiliary qubits needs to be selected after measurement at the end of each iteration. This result is probabilistic, and the overall success rate of the algorithm decreases exponentially as the number of iterations increases if the selection is made after each iteration. Therefore, ensuring an overall higher success rate requires the use of the quantum search algorithm [17,18,19,20] to amplify the amplitude of the target state before measurement. Under ideal conditions of the device, it is proved that the Grover–Long algorithm can achieve a 100% success rate in all cases [17,26,27]. Thus, using the Grover–Long algorithm under ideal conditions to amplify the amplitude, it is possible to obtain a 100% success rate every time. If the complexity of each iteration step is considered synthetically, then the complexity of each iteration step is $O\left(N\right)+O\left(\sqrt{M}\right)$, where M is the dimension of the auxiliary qubits space.

In the future, our quantum algorithms are expected to be combined with finite element methods to solve complex practical problems, such as those related to fluid dynamics [28].