Quantum homotopy perturbation method for nonlinear dissipative ordinary differential equations

Cheng Xue CAS Key Laboratory of Quantum Information, University of Science and Technology of China, Hefei, Anhui, 230026, P. R. China CAS Center For Excellence in Quantum Information and Quantum Physics, University of Science and Technology of China, Hefei, Anhui, 230026, P. R. China Yu-Chun Wu wuyuchun@ustc.edu.cn CAS Key Laboratory of Quantum Information, University of Science and Technology of China, Hefei, Anhui, 230026, P. R. China CAS Center For Excellence in Quantum Information and Quantum Physics, University of Science and Technology of China, Hefei, Anhui, 230026, P. R. China Institute of Artificial Intelligence, Hefei Comprehensive National Science Center, Hefei, Anhui, 230026, P. R. China Guo-Ping Guo gpguo@ustc.edu.cn CAS Key Laboratory of Quantum Information, University of Science and Technology of China, Hefei, Anhui, 230026, P. R. China CAS Center For Excellence in Quantum Information and Quantum Physics, University of Science and Technology of China, Hefei, Anhui, 230026, P. R. China Institute of Artificial Intelligence, Hefei Comprehensive National Science Center, Hefei, Anhui, 230026, P. R. China Origin Quantum Computing Company Limited, Hefei, Anhui, 230026, P. R. China

Abstract

While quantum computing provides an exponential advantage in solving linear differential equations, there are relatively few quantum algorithms for solving nonlinear differential equations. In our work, based on the homotopy perturbation method, we propose a quantum algorithm for solving $n$ -dimensional nonlinear dissipative ordinary differential equations (ODEs). Our algorithm first converts the original nonlinear ODEs into the other nonlinear ODEs which can be embedded into finite-dimensional linear ODEs. Then we solve the embedded linear ODEs with quantum linear ODEs algorithm and obtain a state $\epsilon$ -close to the normalized exact solution of the original nonlinear ODEs with success probability $\Omega(1)$ . The complexity of our algorithm is $O(g\eta T{\rm poly}(\log(nT/\epsilon)))$ , where $\eta$ , $g$ measure the decay of the solution. Our algorithm provides exponential improvement over the best classical algorithms or previous quantum algorithms in $n$ or $\epsilon$ .

Quantum algorithm, nonlinear dissipative ordinary differential equations, homotopy perturbation method

I Introduction

Nonlinear differential equations appear in many fields, such as fluid dynamics, biology, finance, etc. In general, the analytical solutions of nonlinear differential equations cannot be obtained effectively. Numerical methods are often used to solve nonlinear differential equations. However, when solving high-dimensional nonlinear differential equations, too many computing resources are required, which may exceed the computing power of classic computers. It is important to develop more efficient algorithms for solving nonlinear differential equations.

Quantum computing provides a promising way to speed up the solution of various equations. In recent years many quantum algorithms have been developed to solve various equations, such as system of linear equationsHarrow et al. (2009); Childs et al. (2017); Subaşı et al. (2019); Clader et al. (2013), Poisson equationCao et al. (2013), Dirac equationFillion-Gourdeau et al. (2017), heat equationLinden et al. (2020), linear ordinary differential equations (ODEs)Berry (2014); Berry et al. (2017); Xin et al. (2020); Childs and Liu (2020), linear partial ODEsArrazola et al. (2019); Childs et al. (2020) and so onMontanaro and Pallister (2016); Costa et al. (2019); Engel et al. (2019).

However, because of the linearity of quantum mechanics, solving nonlinear equations with quantum computing is challenging, some related algorithms are proposedLeyton and Osborne (2008); Qian et al. (2019); Lubasch et al. (2020); Liu et al. (2021); Lloyd et al. (2020); Budinski (2021); Chen et al. (2021); Xue et al. (2021); Kyriienko et al. (2021). An early quantum algorithm for solving nonlinear ordinary differential equations (ODEs) is proposed in Leyton and Osborne (2008), but the complexity of the algorithm increases exponentially with the evolution time. In Lubasch et al. (2020), the authors proposed a variational quantum algorithm for solving nonlinear differential equations and demonstrate the algorithm by solving 1-dimensional nonlinear Schrodinger equation. However, when the equations become complicated, whether the parameterized quantum circuit used in their work is capable of expressing the solution of the problem and the optimization problem of the parameterized quantum circuit has not yet been concluded. In Liu et al. (2021), a quantum algorithm for solving nonlinear dissipative ODEs is constructed based on Carleman linearizationCarleman et al. (1932); Kowalski and Steeb (1991). This approach embeds the nonlinear ODEs into linear ODEs and solves the linear ODEs with quantum algorithm. The complexity of their algorithm is $O(\frac{T^{2}q}{\epsilon}{\rm poly}(\log(nT/\epsilon)))$ and $q$ measures decay of the solution. In Lloyd et al. (2020), nonlinear ODEs are embedded in Hilbert space, and the evolution of nonlinear ODEs is approximated by the evolution of subsystems in large systems.

Homotopy perturbation methodHe (1999); Babolian et al. (2009); Chakraverty et al. (2019) is a semi-analytical technique for solving linear as well as nonlinear ordinary/partial differential equations. This method, which is a combination of homotopy in topology and classic perturbation techniques, provides us with a convenient way to obtain analytic or approximate solutions for a wide variety of problems arising in different fields, such as Duffing equationHe (2003), nonlinear wave equationsHe (2005) and so on.

In our work, we propose a quantum algorithm for solving time-independent quadratic nonlinear dissipative ODEs. The more general nonlinear ODEs can be reduced to the quadratic ODEs by introducing additional variablesKerner (1981); Forets and Pouly (2017). Our algorithm uses the homotopy perturbation method to transform the problem into a series of nonlinear ODEs which have a special structure. The transformed nonlinear ODEs can be embedded into linear ODEs with a technique similar to Carleman linearization. Then the linear ODEs are solved with quantum linear ODEs algorithm proposed in Berry et al. (2017). Finally, we measure some qubit registers and obtain a state $\epsilon$ -close to the normalized exact solution $u(T)/\|u(T)\|$ at evolution time $T$ .

Our work is similar with Liu’s workLiu et al. (2021), here we list some differences: (1) The truncation method is different, our work uses homotopy perturbation method and Liu et al. (2021) uses Carleman linearization. The convergence condition of homotopy perturbation method and Carleman linearization are different. In Liu et al. (2021), Liu et al. define $R=\frac{\|u_{in}\|\|F_{2}\|}{|Re(\lambda_{1})|}$ (Here we omit the inhomogeneity term $F_{0}$ in their definition) and the convergence condition is $R<1$ . In our work, we define $K=\frac{4\|u_{in}\|\|F_{2}\|}{|Re(\lambda_{1})|}=4R$ , the convergence condition is $K<1$ , the factor ’ $4$ ’ is caused by the homotopy perturbation method. By Corollary $\bm{1}$ in Liu et al. (2021) and Lemma 9, the truncation errors of Carleman linearization and homotopy perturbation method decrease exponentially with the truncation order $c$ as $\|u_{in}\|R^{c}(1-e^{Re(\lambda_{1})t})^{c}$ and $\frac{K^{c+2}}{1-K}$ , respectively. (2)The dependence of our algorithm on the error $\epsilon$ and evolution time $T$ in our algorithm is $O(T\rm{poly}(\log(1/\epsilon)))$ , which is $O(T^{2}/\epsilon)$ in Liu et al. (2021). Therefore the complexity of our algorithm is exponentially improved on $\epsilon$ compared to Liu et al. (2021), the cost of this improvement is a stronger constraint on the problem, i.e., $K<\sqrt{2}/2$ and $\frac{(c+1)|Re(\lambda_{1})|}{\|F_{2}\|}\leq 1$ .

This paper is organized as follows. Sect.II introduces the quadratic ODEs to be solved. Then we show the details of quantum homotopy perturbation method in Sect.III. Initial state preparation and oracle construction of matrix $A$ are discussed in Sect.IV. In the following three sections we analyze our method from different aspects: Sect.V gives an upper bound of the condition number of the linear system to be solved. Sect.VI analyzes the solution error of our method. Sect.VII gives a lower bound of success probability. Next, the main result of our work is proved in Sect.VIII. Finally, we conclude our work with a discussion of the result and some open problems in Sect.IX.

II Quadratic ODEs

We focus on an initial value problem described by the $n$ -dimensional quadratic ODEs. The problem to be solved is defined as

\frac{du}{dt}=F_{1}u+F_{2}u^{\otimes 2},\ u(0)=u_{in},

(1)

where $u$ , $u_{in}\in\mathbbm{R}^{n}$ , $F_{1}\in\mathbbm{R}^{n\times n}$ , $F_{2}\in\mathbbm{R}^{n\times n^{2}}$ are time independent sparse matrices. The sparsity of $F_{1}$ , $F_{2}$ is $s$ , which means the number of non-zero elements in each row or column of $F_{1},F_{2}$ does not exceed $s$ . We assume $F_{1}$ is a normal matrix and the eigenvalues $\lambda_{i}$ of $F_{1}$ satisfy $Re(\lambda_{n})\leq\dots\leq Re(\lambda_{1})<0$ . We also assume oracles $O_{F1}$ , $O_{F2}$ and $O_{u}$ are given, $O_{F1}$ , $O_{F2}$ are used to extract the non-zero position and value of $F_{1}$ , $F_{2}$ respectively, and $|u_{in}\rangle$ is prepared with $O_{u}$ . In specific, $O_{F1}$ , $O_{F2}$ and $O_{u}$ are defined as

$\displaystyle O_{F11}\|j\rangle\|k\rangle=\|j\rangle\|f_{1}(j,k)\rangle,$	$\displaystyle 0\leq j\leq n-1,\ 0\leq k\leq s-1,$	(2)
$\displaystyle O_{F12}\|j\rangle\|k\rangle\|z\rangle=\|j\rangle\|k\rangle\|z\oplus{F_{1}}_{j,k}\rangle,$	$\displaystyle 0\leq j,\ k\leq n-1,$
$\displaystyle O_{F13}\|j\rangle\|0\rangle=\|j\rangle\|g(j)\rangle,$	$\displaystyle 0\leq j\leq n-1,$
$\displaystyle O_{F21}\|j\rangle\|k\rangle=\|j\rangle\|f_{2}(j,k)\rangle,$	$\displaystyle 0\leq j\leq n-1,\ 0\leq k\leq s-1,$
$\displaystyle O_{F22}\|j\rangle\|k\rangle\|z\rangle=\|j\rangle\|k\rangle\|z\oplus{F_{2}}_{j,k}\rangle,$	$\displaystyle 0\leq j\leq n-1,\ 0\leq k\leq n^{2}-1,$
$\displaystyle O_{u}\|0\rangle=\|u_{in}/\\|u_{in}\\|\rangle,$

where $f_{1}(j,k)$ and $f_{2}(j,k)$ represent the column number of the $k$ -th non-zero element in $j$ -th row of $F_{1}$ , $F_{2}$ respectively, $g(j)$ satisfies $f_{1}(j,g(j))=j$ , here we treat the diagonal element of $F_{1}$ as a non-zero element. $O_{F13}$ is used to construct an oracle of a matrix related to $F_{1}$ , the details are shown in the proof of Lemma 5.

We define a parameter $K$ which characterizes the nonlinearity of Eq.(1),

K:=\frac{4\|u_{in}\|\|F_{2}\|}{|Re(\lambda_{1})|}.

(3)

We assume $K\geq\|u_{in}\|$ , if this is not satisfied, we rescale $u$ to $\zeta u$ with a suitable constant $\zeta$ which keeps $\frac{4\|u_{in}\|\|F_{2}\|}{|Re(\lambda_{1})|}$ unchanged and makes $K\geq\|u_{in}\|$ . In this paper we use spectral norm, it means $\|\cdot\|=\|\cdot\|_{2}$ .

III Quantum Homotopy Perturbation Method

In this section, we give the whole process of quantum homotopy perturbation method for solving Eq.(1). It contains four steps:

(1)

Using homotopy perturbation method to transform Eq.(1) into Eq.(III.1), a series of nonlinear ODEs with variable $\nu_{i}$ .
(2)

Embedding Eq.(III.1) into Eq.(8), linear ODEs with variable $\vec{y}$ .
(3)

Solving Eq.(8) with quantum algorithm proposed in Berry et al. (2017).
(4)

Measurement.

The following four subsections introduce the details of these four steps.

III.1 Homotopy perturbation method

Firstly, we introduce the process of homotopy perturbation methodHe (1999); Babolian et al. (2009); Chakraverty et al. (2019) for solving Eq.(1). Using homotopy method, we construct homotopy $\nu(t,p):\mathbbm{R^{+}}\times[0,1]\to\mathbbm{R}^{n}$ , which satisfies

H(\nu,p)=\frac{d\nu}{dt}-F_{1}\nu-pF_{2}\nu^{\otimes 2}=0,\ \nu(0,p)=u_{in}.

(4)

Assuming $\nu$ is represented as

\nu=\nu_{0}+p\nu_{1}+p^{2}\nu_{2}+\dots+p^{c}\nu_{c}.

(5)

Substituting Eq.(5) into Eq.(4), then equating the terms with identical powers of $p$ , we have the following equations:

		$\displaystyle\frac{d\nu_{0}}{dt}-F_{1}\nu_{0}=0,\ \nu_{0}(0)=u_{in},$
		$\displaystyle\frac{d\nu_{1}}{dt}-F_{1}\nu_{1}-F_{2}\nu_{0}\otimes\nu_{0}=0,\ \nu_{1}(0)=0,$
		$\displaystyle\frac{d\nu_{2}}{dt}-F_{1}\nu_{2}-F_{2}(\nu_{0}\otimes\nu_{1}+\nu_{1}\otimes\nu_{0})=0,\ \nu_{2}(0)=0,$
		$\displaystyle\dots$
		$\displaystyle\frac{d\nu_{c}}{dt}-F_{1}\nu_{c}-F_{2}\sum_{j=0}^{c-1}{\nu_{j}\otimes\nu_{c-1-j}}=0,\ \nu_{c}(0)=0.$		(6)

When $p=1$ in Eq.(4), we have

\tilde{u}=\nu_{0}+\nu_{1}+\nu_{2}+\dots+\nu_{c}.

(7)

The difference between $\tilde{u}$ and $u$ is analyzed in Sect.VI.1.

III.2 Linear embedding

Secondly, Eq.(III.1) is embedded into the linear ODEs defined in Eq.(8):

\frac{d\vec{y}}{dt}=A\vec{y},\ \vec{y}(0)=y_{in},

(8)

where $\vec{y}=[\vec{y}_{0},\vec{y}_{1},\dots,\vec{y}_{c}]$ , $\vec{y}_{i}$ satisfies

\vec{y}_{i}=\left\{\begin{array}[]{ll}[\nu_{0}+\nu_{1}+\dots+\nu_{c}],&i=0,\\ \left[\vec{y}_{i,0},\vec{y}_{i,1},\dots,\vec{y}_{i,\beta_{i}-1}\right],&1\leq i\leq c,\end{array}\right.

(9)

where $\beta_{i}$ represents the number of items in $\vec{y}_{i}$ , $\vec{y}_{i,j}$ represents $j$ -th item of $\vec{y}_{i}$ , it is represented as $\vec{y}_{i,j}=\otimes_{k=0}^{i}{\nu_{a_{i,j,k}}}$ , $a_{i,j,k}$ satisfies

a_{i,j,k}\geq 0,\ i+1\leq\sum_{k=0}^{i}{(a_{i,j,k}+1)}\leq c+1.

(10)

By Eq.(10), $\beta_{i}$ satisfies

\beta_{i}=\left\{\begin{array}[]{ll}1,&i=0,\\ \sum_{k=i}^{c}{\binom{k}{i}},&1\leq i\leq c.\end{array}\right.

(11)

We set $\vec{y}_{i,0}=\nu_{0}^{\otimes{i+1}},\ i=1,2,\dots,c$ , then by Eq.(III.1), $y_{in}$ is written as

y_{in}=[[u_{in}],[u_{in}^{\otimes 2},0,\dots,0],[u_{in}^{\otimes 3},0,\dots,0]\dots,[u_{in}^{\otimes c+1}]].

(12)

We define $\vec{a}_{i,j}=[a_{i,j,0},a_{i,j,1},\dots,a_{i,j,i}]$ , the mapping $\vec{a}_{i,j}\mapsto j$ is one-to-one mapping, so we can construct the following two operations

O_{a1}|\vec{a}_{i,j}\rangle|0\rangle=|\vec{a}_{i,j}\rangle|j\rangle,

(13)

O_{a2}|i\rangle|j\rangle|0\rangle=|i\rangle|j\rangle|\vec{a}_{i,j}\rangle

(14)

with $O(c)$ -qubit quantum arithmetic circuit, and the gate complexity is $O({\rm poly}(c))$ Nielsen and Chuang (2002). $O({\rm poly}(c))$ will not influence the complexity of our algorithm, so the complexity of $O_{a1}$ and $O_{a2}$ can be ignored in the following analysis. The dimension of $\vec{y}_{i}$ is $n^{i+1}\beta_{i}$ , so the dimension of $\vec{y}$ is

N=\sum_{i=0}^{c}{n^{i+1}\beta_{i}}=(n+1)^{c+1}-1-cn\approx(n+1)^{c+1}.

(15)

Next we analyze the structure of matrix $A$ . We have

		$\displaystyle\frac{d\vec{y}_{i,j}}{dt}=(\sum_{j=0}^{i}{I_{n}^{\otimes j}\otimes F_{1}\otimes I_{n}^{\otimes i-j}})\vec{y}_{i,j}$
		$\displaystyle+\sum_{k=0}^{i}{(I_{n}^{k}\otimes F_{2}\otimes I_{n}^{\otimes{i-k}})\nu_{a_{i,j,0}}\otimes\dots\otimes\nu_{a_{i,j,k-1}}\otimes(\sum_{l=0}^{a_{i,j,k}-1}{\nu_{l}\otimes\nu_{a_{i,j,k}-1-l}})\otimes\nu_{a_{i,j,k+1}}\otimes\dots\otimes\nu_{a_{i,j,i}}},$		(16)

where $\nu_{a_{i,j,0}}\otimes\dots\otimes\nu_{a_{i,j,k-1}}\otimes\nu_{l}\otimes\nu_{a_{i,j,k}-1-l}\otimes\nu_{a_{i,j,k+1}}\otimes\dots\otimes\nu_{a_{i,j,i}}\in\vec{y}_{i+1}$ , so Eq.(8) can be written as

\frac{\mathrm{d}}{\mathrm{d}t}\left(\begin{array}[]{c}\vec{y}_{0}\\ \vec{y}_{1}\\ \vdots\\ \vec{y}_{c-1}\\ \vec{y}_{c}\end{array}\right)=\left(\begin{array}[]{cccccc}A_{0,0}&A_{0,1}&&&\\ &A_{1,1}&A_{1,2}&&\\ &&\ddots&\ddots&\\ &&&A_{c-1,c-1}&A_{c-1,c}\\ &&&&A_{c,c}\end{array}\right)\left(\begin{array}[]{c}\vec{y}_{0}\\ \vec{y}_{1}\\ \vdots\\ \vec{y}_{c-1}\\ \vec{y}_{c}\end{array}\right),

(17)

$A_{i,i}$ is $n^{i+1}\beta_{i}$ dimensional square matrix, which is represented as

A_{i,i}=I_{\beta_{i}}\otimes(\sum_{j=0}^{i}{I_{n}^{j}\otimes F_{1}\otimes I_{n}^{\otimes i-j}}),

(18)

$A_{i,i+1}$ is $n^{i+1}\beta_{i}\times n^{i+2}\beta_{i+1}$ dimensional matrix, its elements are determined by Eq.(III.2). $|y(t)\rangle$ is defined to represent $\vec{y}(t)$ :

|y(t)\rangle=\sum_{i=0}^{c}{\sum_{j=0}^{\beta_{i}-1}{|i,j\rangle|y_{i,j}(t)\rangle}}.

(19)

III.3 Quantum linear ODEs algorithm

Thirdly, Eq.(8) is solved with the quantum algorithm proposed in Berry et al. (2017). $\vec{y}(t)$ is written as

\vec{y}(t)=e^{At}\vec{y}(0).

(20)

We define $T_{k}(z):=\sum_{j=0}^{k}{\frac{z^{j}}{j!}}$ . When $k$ is large enough and the evolution time $h$ is relatively short (for example, $h\leq 1/\|A\|$ ), we have $\vec{y}(h)\approx T_{k}(Ah)\vec{y}(0)$ . This approximate solution can be used as the initial condition for the next approximation, repeating this procedure $m$ steps we have the approximation of $\vec{y}(mh)$ .

Next we introduce the details of the algorithm proposed in Berry et al. (2017). Let $m,k,p\in\mathbbm{Z}^{+}$ and define

	$\displaystyle C_{m,k,p}(A):=$	$\displaystyle\sum_{j=0}^{d}\|j\rangle\langle j\|\otimes I-\sum_{i=0}^{m-1}\sum_{j=1}^{k}\|i(k+1)+j\rangle\langle i(k+1)+j-1\|\otimes A/j$		(21)
		$\displaystyle-\sum_{i=0}^{m-1}\sum_{j=0}^{k}\|(i+1)(k+1)\rangle\langle i(k+1)+j\|\otimes I-\sum_{j=d-p+1}^{d}\|j\rangle\langle j-1\|\otimes I,$		(21)

where $d:=m(k+1)+p$ , $I$ is an $N$ -dimensional unit matrix. We consider the linear system

C_{m,k,p}(Ah)|x\rangle=|0\rangle|y_{in}\rangle,

(22)

where $|y_{in}\rangle\in\mathbbm{C}^{N},h\in\mathbbm{R}^{+}$ . After evolving $m$ steps, the approximate solution of $k$ -order Taylor series is obtained, and the solution remains unchanged at $p$ steps. The solution of Eq.(22) is represented as $|x\rangle=C_{m,k,p}(Ah)^{-1}|0\rangle\left|y_{in}\right\rangle$ , it can also be written as

|x\rangle=\sum_{i=0}^{m-1}\sum_{j=0}^{k}|i(k+1)+j\rangle\left|x_{i,j}\right\rangle+\sum_{j=0}^{p}|m(k+1)+j\rangle\left|x_{m,j}\right\rangle.

(23)

By Eq.(21), $|x_{i,j}\rangle$ satisfies

\begin{array}[]{l}|x_{0,0}\rangle=|y_{in}\rangle,\\ |x_{i,0}\rangle=\sum_{j=0}^{k}|x_{i-1,j}\rangle,\quad 1\leq i\leq m,\\ |x_{i,1}\rangle=Ah|x_{i,0}\rangle,\quad 0\leq i<m,\\ |x_{i,j}\rangle=(Ah/j)|x_{i,j-1}\rangle,\quad 0\leq i<m,2\leq j\leq k,\\ |x_{m,j}\rangle=|x_{m,j-1}\rangle,\quad 1\leq j\leq p.\end{array}

(24)

Then we have

$\displaystyle\left\|x_{0,0}\right\rangle$	$\displaystyle=\left\|y_{in}\right\rangle,$	(25)
$\displaystyle\left\|x_{0,j}\right\rangle$	$\displaystyle=\left((Ah)^{j}/j!\right)\left\|x_{0,0}\right\rangle,\quad 1\leq j\leq k,$
$\displaystyle\left\|x_{1,0}\right\rangle$	$\displaystyle=T_{k}(Ah)\left\|x_{0,0}\right\rangle\approx\exp(Ah)\left\|y_{in}\right\rangle,$
$\displaystyle\left\|x_{1,j}\right\rangle$	$\displaystyle=\left((Ah)^{j}/j!\right)\left\|x_{1,0}\right\rangle,\quad 1\leq j\leq k,$
$\displaystyle\left\|x_{2,0}\right\rangle$	$\displaystyle=T_{k}(Ah)\left\|x_{1,0}\right\rangle\approx\exp(2Ah)\left\|y_{in}\right\rangle,$
	$\displaystyle\vdots$
$\displaystyle\left\|x_{m-1,0}\right\rangle$	$\displaystyle=T_{k}(Ah)\left\|x_{m-2,0}\right\rangle\approx\exp(Ah(m-1))\left\|y_{in}\right\rangle,$
$\displaystyle\left\|x_{m-1,j}\right\rangle$	$\displaystyle=\left((Ah)^{j}/j!\right)\left\|x_{m-1,0}\right\rangle,\quad 1\leq j\leq k,$
$\displaystyle\left\|x_{m,0}\right\rangle$	$\displaystyle=T_{k}(Ah)\left\|x_{m-1,0}\right\rangle\approx\exp(Ahm)\left\|y_{in}\right\rangle,$
$\displaystyle\left\|x_{m,j}\right\rangle$	$\displaystyle=\left\|x_{m,0}\right\rangle\approx\exp(Ahm)\left\|y_{in}\right\rangle,\quad 1\leq j\leq p,$

$|x_{i,0}\rangle$ is the approximate solution of the system at time $ih$ , $i=\in\{0,1,2,\dots,m\}$ . $|x_{m,0}\rangle=|x_{m,1}\rangle=\dots=|x_{m,p}\rangle$ is the approximate solution of $\vec{y}(t)=e^{At}y_{in}$ at $t=mh$ .

III.4 Measurement

Finally, we measure some qubit registers of $|x\rangle$ and get a state $\epsilon$ -close to the normalized solution of Eq.(1). The measurement is divided into two steps: (1) Measure the first qubit register of $|x\rangle$ which is defined in Eq.(23), if the result is $|m(k+1)+j\rangle,j=0,1,\dots,p$ , we have $|y(t)\rangle$ in the second qubit register of $|x\rangle$ . (2) Measure the first qubit register of $|y(t)\rangle$ which is defined in Eq.(19), if the result is $|0,0\rangle$ , we get a state $\epsilon$ -close to $|u(t)/\|u(t)\|\rangle$ in the qubit second register of $|y(t)\rangle$ .

This measurement step is probabilistic, the success probability is analyzed in Sect.VII.

IV State Preparation and Oracle Construction

In this section, we give the preparation of $|y_{in}\rangle$ and oracle construction of $A$ .

IV.1 State preparation

We first discuss the preparation of $|y_{in}\rangle$ , the result is shown in Lemma 1.

Lemma 1.

By Eq.(12), $|y_{in}\rangle$ is defined as

|y_{in}\rangle=\frac{1}{\sqrt{M}}\sum_{i=0}^{c}{|i,0\rangle|u_{in}^{\otimes i+1}\rangle},

(26)

where $M=\sum_{j=i}^{c}{\|u_{in}\|^{2(i+1)}}$ . Given $O_{u}$ defined in Eq.(2), $|y_{in}\rangle$ can be prepared by querying $O_{u}$ $O(c)$ times.

Proof.

First we prepare

|\psi\rangle=\frac{1}{\sqrt{M}}\sum_{i=0}^{c}{\|u_{in}\|^{2(i+1)}|i,0\rangle}.

(27)

Then we execute controlled $O_{u}$ operation $C-O_{u}=\sum_{i=0}^{c}{|i,0\rangle\langle i,0|\otimes O_{u}^{\otimes i+1}}$ on $|\psi\rangle$ ,

|y_{in}\rangle=\frac{1}{\sqrt{M}}\sum_{i=0}^{c}{|i,0\rangle|u_{in}^{\otimes i+1}\rangle}.

(28)

The query complexity of $O_{u}$ is $O(c)$ .

IV.2 Oracle construction of $A$

Before introducing oracle construction of $A$ , we analyze some features of $A$ , including sparsity, upper bound of $\|A\|$ and eigenvalue of $A$ . The results are shown in Lemma 2, Lemma 3 and Lemma 4.

Lemma 2.

The sparsity of the matrix $A$ is $O(sc^{2})$ .

Proof.

The sparsity $A_{i,i}$ is $(i+1)sc$ . The sparsity of $A_{0,1}$ is $sc(c+1)/2$ , when $i\geq 1$ , the sparsity of $A_{i,i+1}$ is $\max\{s(c-i),s(i+1)\}$ . Therefore, the sparsity of matrix $A$ is $O(sc^{2})$ .

Lemma 3.

$\|A\|$ satisfies $\|A\|\leq(c+1)(\|F_{1}\|+\|F_{2}\|)$ .

Proof.

By the definition of $A_{i,i}$ , we have

\|A_{i,i}\|\leq(c+1)\|F_{1}\|,\quad i\in[c+1]_{0},

(29)

$[c+1]_{0}=[0,1,\dots,c]$ and in this paper, for any $i\in\mathbbm{N}$ , $[i+1]_{0}=[0,1,\dots,i]$ . Then we analyze upper bound of $\|A_{i,i+1}\|$ . We have $A_{0,1}A_{0,1}^{T}=\frac{c(c+1)}{2}F_{2}F_{2}^{T}$ , then $\|A_{0,1}\|$ satisfies

\displaystyle\|A_{0,1}\|=\sqrt{\frac{c(c+1)}{2}}\|F_{2}\|<(c+1)\|F_{2}\|.

(30)

When $i\geq 1$ , $A_{i,i+1}$ can be regarded as $\beta_{i}\times\beta_{i+1}$ dimensional block matrix, each block unit is an $n^{i+1}\times n^{i+2}$ dimensional matrix and has the form $I^{\otimes j}\otimes F_{2}\otimes I^{\otimes i-j}$ , $j\in[i+1]_{0}$ . From the structure of $A_{i,i+1}$ , the number of non-zero block unit in each row or column of $A_{i,i+1}$ is no more than $c$ , so $A_{i,i+1}$ can be divided into at most $c$ matrices which have only one non-zero block unit in each row or column. Therefore,

\|A_{i,i+1}\|\leq c\|F_{2}\|,\quad i\in[1,2,\dots,c-1].

(31)

Combining Eq.(29), Eq.(30) and Eq.(31), $\|A\|$ satisfies

$\displaystyle\\|A\\|\leq$	$\displaystyle\\|{\rm diag}(A_{0,0},A_{1,1},\dots,A_{c,c})\\|+\\|{\rm diag}(A_{0,1},A_{1,2},\dots,A_{c-1,c})\\|$
$\displaystyle=$	$\displaystyle\max_{i=0}^{c}\{\\|A_{i,i}\\|\}+\max_{i=0}^{c-1}\{\\|A_{i,i+1}\\|\}$
$\displaystyle\leq$	$\displaystyle(c+1)(\\|F_{1}\\|+\\|F_{2}\\|).$	(32)

Lemma 4.

The eigenvalue $\gamma_{i}$ of matrix $A$ satisfies $Re(\gamma_{i})<0,i\in[N]_{0}$ .

Proof.

From the structure of $A$ , the eigenvalues of $A$ are the sets of eigenvalues of $A_{i,i}$ for $i\in[c+1]_{0}$ . The eigenvalue of $A_{i,i}$ is the sum of $i+1$ eigenvalues of $F_{1}$ . For any $j\in[1,2,\dots,n]$ , eigenvalue of $F_{1}$ $\lambda_{j}$ satisfies $Re(\lambda_{j})<0$ , so the real part of all eigenvalues of $A_{i,i}$ is less than 0. Therefore, the eigenvalue $\gamma_{i}$ of $A$ satisfies $Re(\gamma_{i})<0,i\in[N]_{0}$ .

Next, we introduce the way to construct oracle $O_{A}$ of $A$ . $O_{A}$ gives the way to extract non-zero element position and value of $A$ . We first give Lemma 5.

Lemma 5.

Let matrix $B(m)=\sum_{j=0}^{m}{I_{n}^{j}\otimes F_{1}\otimes I_{n}^{\otimes m-j}}$ . Oracles $O_{m,1}$ and $O_{m,2}$ are defined as

		$\displaystyle O_{m,1}\|\vec{j}\rangle\|l\rangle=\|\vec{j}\rangle\|g_{m}(\vec{j},l)\rangle,$		(33)
		$\displaystyle O_{m,2}\|\vec{j}\rangle\|\vec{k}\rangle\|z\rangle=\|\vec{j}\rangle\|\vec{k}\rangle\|z\oplus B(m)_{\vec{j},\vec{k}}\rangle,$		(33)

where $\vec{j}=j_{m}j_{m-1}\dots j_{1}j_{0}$ , $j_{i}\in[n]_{0}$ . $g_{m}(\vec{j},l)$ represents the column number of $l$ -th non-zero element in $\vec{j}$ -th row of $B(m)$ , $g_{m}(\vec{j},l)$ is also written as $g_{m}(\vec{j},l)=k_{m}k_{m-1}\dots k_{0}$ , $k_{i}\in[n]_{0}$ . Then $O_{m,1}$ , $O_{m,2}$ can be constructed by querying $O_{F1}$ $m+1$ times.

Proof.

When $m=0$ , $B(0)=F_{1}$ , $O_{0,1}$ , $O_{0,2}$ can be constructed by querying $O_{F1}$ once. Assuming when $m\geq 1$ , Oracles $O_{m-1,1}$ , $O_{m-1,2}$ of $B(m-1)$ are constructed by querying $O_{F1}$ $m$ times. $B(m)$ is also written as

B(m)=I\otimes B(m-1)+F_{1}\otimes I^{\otimes m}.

(34)

Generally, we treat the diagonal element of $F_{1}$ as a non-zero element, so the sparsity of $B(m)$ is $(m+1)c-m$ . $g_{m}(\vec{j},k)$ can be represented as

g_{m}(\vec{j},k)=\left\{\begin{array}[]{ll}f_{1}(j_{m},k)j_{m-1}\dots j_{0},&0\leq k<g(j_{m}),\\ j_{m}g_{m-1}(j_{m-1}\dots j_{0},k-g(j_{m})),&0\leq k-g(j_{m})\leq m(c-1)+1,\\ f_{1}(j_{m},k-m(c-1))j_{m-1}\dots j_{0},&g(j_{m})+1\leq k-m(c-1)\leq c-1,\end{array}\right.

(35)

where $f_{1}(j,k)$ and $g(j)$ are defined in Eq.(2). Then $O_{m,1}$ can be constructed with Eq.(35), the construction process needs to query $O_{m-1,1}$ , $O_{F11}$ , $O_{F13}$ once each.

On the other hand, the element of $B(m)$ is written as

B(m)_{\vec{j},\vec{k}}=\delta_{j_{m}k_{m}}B(m-1)_{j_{m-1}\dots j_{0},k_{m-1}\dots k_{0}}+\delta_{j_{m-1}\dots j_{0},k_{m-1}\dots k_{0}}B_{j_{m},k_{m}}.

(36)

By Eq.(36), $O_{m,2}$ can be constructed by querying $O_{m-1,2}$ and $O_{F12}$ once each. Therefore $O_{m,1}$ and $O_{m,2}$ can be constructed by querying $O_{m-1,1}$ , $O_{m-1,2}$ , $O_{F11}$ , $O_{F12}$ and $O_{F13}$ once each.

From the above analysis, $O_{m,1}$ and $O_{m,2}$ can be constructed by querying $O_{F1}$ $m+1$ times.

Lemma 6.

The oracle $O_{A}$ of $A$ can be constructed by querying $O_{F1}$ $O(c)$ times and querying $O_{F2}$ $O(1)$ times.

Proof.

To construct $O_{A}$ , we need to construct oracles of $A_{i,i}$ and $A_{i,i+1}$ . We first consider $A_{i,i}$ , we define $B(i)=\sum_{j=0}^{i}{I_{n}^{j}\otimes F_{1}\otimes I_{n}^{\otimes c-i-j}}$ , by Lemma 5, the oracle of $B(i)$ can be constructed by querying $O_{F1}$ $i+1$ times. By Eq.(18), the oracle of $A_{i,i}$ can be constructed by querying oracle of $B(i)$ once.

Next we consider $A_{i,i+1}$ . When $i=0$ , $A_{0,1}=[F_{2},F_{2},\dots,F_{2}]$ , the oracle of $A_{0,1}$ is constructed by querying $O_{F2}$ once. When $i\geq 1$ , $A_{i,i+1}$ is also regarded as $\beta_{i}\times\beta_{i+1}$ dimensional block matrix $D(i)$ , the dimension of each block unit is $n^{i+1}\times n^{i+2}$ . The number of non-zero block unit in $j$ -th row of $D(i)$ is $\sum_{k=0}^{i}{a_{i,j,k}}$ . Consider the $l$ -th non-zero bolck unit in $j$ -th row of $D(i)$ , $l$ is represented as $l=\sum_{k=0}^{i_{1}}{a_{i,j,k}}+i_{2}$ , where $i_{1}\in[i]_{0},i_{2}\in\mathbbm{N}$ , the $l$ -th non-zero block unit is $I^{\otimes i_{1}}\otimes F_{2}\otimes I^{\otimes i-i_{1}}$ , and the corresponding $\vec{a}_{i+1,j^{{}^{\prime}}}$ is represented as

\vec{a}_{i+1,j^{{}^{\prime}}}=[a_{i,j,0},\dots,a_{i,j,i_{1}},i_{2},a_{i,j,i_{1}+1}-1-i_{2},a_{i,j,i+2},\dots a_{i,j,i}],

(37)

$j^{{}^{\prime}}$ can be obtained from $\vec{a}_{i+1,j^{{}^{\prime}}}$ with $O_{a1}$ . The oracle of the non-zero block unit $I^{\otimes i_{1}}\otimes F_{2}\otimes I^{\otimes i-i_{1}}$ is constructed by querying $O_{F2}$ once. By realizing the above process with quantum circuit, we construct an oracle that extracts the non-zero position of $A_{i,i+1}$ . The specific implementation process is

	$\displaystyle\|j,b_{i},b_{i-1},\dots,b_{i_{1}},\dots,b_{0}\rangle\|l,l_{b}\rangle\|0,0\rangle$
$\displaystyle\xrightarrow{O_{a2}}$	$\displaystyle\|j,b_{i},b_{i-1},\dots,b_{i_{1}},\dots,b_{0}\rangle\|l,l_{b}\rangle\|i_{1},i_{2}\rangle$
$\displaystyle\xrightarrow{O_{a1}}$	$\displaystyle\|j,b_{i},b_{i-1},\dots,b_{i_{1}},\dots,b_{0}\rangle\|j^{{}^{\prime}},l_{b}\rangle\|i_{1},i_{2}\rangle$
$\displaystyle\xrightarrow{O_{F21}}$	$\displaystyle\|j,b_{i},b_{i-1},\dots,b_{i_{1}},\dots,b_{0}\rangle\|j^{{}^{\prime}},b_{i},b_{i-1},\dots,f_{2}(b_{i_{1}},l_{b}),\dots b_{0}\rangle\|i_{1},i_{2}\rangle$
$\displaystyle\xrightarrow{uncompute}$	$\displaystyle\|j,b_{i},b_{i-1},\dots,b_{i_{1}},\dots,b_{0}\rangle\|j^{{}^{\prime}},b_{i},b_{i-1},\dots,f_{2}(b_{i_{1}},l_{b}),\dots b_{0}\rangle\|0,0\rangle.$	(38)

There are some ancilla qubits to represent $|\vec{a}_{i,j}\rangle$ and $|\vec{a}_{i+1,j^{{}^{\prime}}}\rangle$ in the process shown in Eq.(Proof). For simplicity, we ignore the compute and uncompute process of $|\vec{a}_{i,j}\rangle$ and $|\vec{a}_{i+1,j^{{}^{\prime}}}\rangle$ .

The oracle that extracts the non-zero value of $A_{i,i+1}$ can also be constructed in a similar process. For any $j\in[\beta_{i}]_{0}$ , $k\in[\beta_{i+1}]_{0}$ , we can judge whether $D(i)_{j,k}$ is a non-zero block unit and use $l_{j,k}$ to represent the judgement. If $D(i)_{j,k}$ is a non-zero block unit, we can find $i_{1}$ and represent it as $I^{\otimes i_{1}}\otimes F_{2}\otimes I^{\otimes i-i_{1}}$ , then we use $O_{F22}$ to extract the elements of the non-zero block unit. The whole process is shown as

	$\displaystyle\|\vec{j}\rangle\|\vec{k}\rangle\|0,0\rangle\|z\rangle=\|j,b_{i},b_{i-1},\dots,b_{0}\rangle\|k,b_{i}^{{}^{\prime}},b_{i-1}^{{}^{\prime}},\dots,b_{0}^{{}^{\prime}}\rangle\|0,0\rangle\|z\rangle$
$\displaystyle\xrightarrow{O_{a1},O_{a2}}$	$\displaystyle\|\vec{j}\rangle\|\vec{k}\rangle\|l_{j,k},i_{1}\rangle\|z\rangle$
$\displaystyle\xrightarrow{O_{F22}}$	$\displaystyle\|\vec{j}\rangle\|\vec{k}\rangle\|l_{j,k},i_{1}\rangle\|z\oplus l_{j,k}\times{F_{2}}_{b_{i_{1}},b_{i_{1}}^{{}^{\prime}}}\prod_{l\in[i],l\neq i_{1}}{\delta_{b_{l}b_{l}^{{}^{\prime}}}}\rangle$
$\displaystyle\xrightarrow{uncompute}$	$\displaystyle\|\vec{j}\rangle\|\vec{k}\rangle\|0,0\rangle\|z\oplus D(i)_{\vec{j},\vec{k}}\rangle.$	(39)

The query complexity of $O_{F22}$ in the above process is $O(1)$ . Therefore, oracle of $A_{i,i+1}$ can be constructed by querying $O_{F2}$ $O(1)$ times.

After constructing the oracles of $A_{i,i}$ and $A_{i,i+1}$ , the oracle of $A$ can be directly constructed by querying oracles of $A_{i,i}$ and $A_{i,i+1}$ once. So the oracle $O_{A}$ can be constructed by querying $O_{F1}$ $O(c)$ times and querying $O_{F2}$ $O(1)$ times.

V Condition Number

In this section, we give an upper bound of the condition number of $C_{m,k,p}(Ah)$ defined in Sect.III.3. We first analyze the upper bound of $\|e^{Aht}\|$ , we have the following lemma.

Lemma 7.

Consider the matrix $A$ defined in Sect.III.2, when

\frac{(c+1)\|F_{2}\|}{|Re(\lambda_{1})|}\leq 1,

(40)

for $t\geq 0$ , $\|e^{At}\|$ satisfies

\|e^{At}\|\leq c+1.

(41)

Proof.

We consider $A$ as a $c+1$ -dimensional block matrix. $A$ is divided into

A=B+C,

(42)

$B$ contains $A_{i,i},i=0,1,\dots,c$ and $C$ contains $A_{i,i+1},i=0,1,\dots,c-1$ . We analyze the upper bound of $\|e^{At}\|$ according to the method introduced in Van Loan (1977), $e^{At}$ is written as

e^{(B+C)t}=e^{Bt}+\int_{0}^{t}{e^{B(t-t_{0})}Ce^{(B+C)t_{0}}dt_{0}}.

(43)

Using this formula to expand $e^{(B+C)t_{0}}$ we obtain

e^{(B+C)t}=e^{Bt}+\int_{0}^{t}{e^{B(t-t_{0})}Ce^{Bt_{0}}dt_{0}}+\int_{0}^{t}{\int_{0}^{t_{0}}{e^{B(t-t_{0})}Ce^{B(t_{0}-t_{1})}Ce^{(B+C)t_{1}}dt_{1}}dt_{0}}.

(44)

Clearly, a repetition of this process gives

e^{(B+C)t}=e^{Bt}+\sum_{k=0}^{c-1}{A_{k}(t)}+R_{c}(t),

(45)

where

A_{k}(t)=\int_{0}^{t}\int_{0}^{t_{0}}\cdots\int_{0}^{t_{k-1}}e^{B\left(t-t_{0}\right)}Ce^{B\left(t_{0}-t_{1}\right)}C\cdots Ce^{Bt_{k}}dt_{k}\cdots dt_{0}

(46)

and

R_{c}(t)=\int_{0}^{t}\int_{0}^{t_{0}}\cdots\int_{0}^{t_{c-1}}e^{B\left(t-t_{0}\right)}C\cdots Ce^{B\left(t_{c-1}-t_{c}\right)}Ce^{(B+C)t_{c}}dt_{c}\cdots dt_{0}.

(47)

Noting that

\|A_{k}(t)\|\leq\|e^{Bt}\|\|C\|^{k+1}\times\frac{t^{k+1}}{(k+1)!}=\|e^{Bt}\|\frac{(\|Ct\|)^{k+1}}{(k+1)!},k=0,1,\dots,c-1,

(48)

\|e^{Bt}\|=e^{Re(\lambda_{1})t},\ \|C\|\leq(c+1)\|F_{2}\|,

(49)

and the the matrix $[e^{B(t-t_{0})}C]\dots[e^{B(t_{c-1}-t_{c})}C]$ is zero beacuse it is the product of $c+1$ , $(c+1)\times(c+1)$ strictly upper triangular block matrices and thus, $R_{c}(t)=0$ . Hence,

\|e^{At}\|=\|e^{(B+C)t}\|\leq e^{Re(\lambda_{1})t}\sum_{k=0}^{c}{\frac{((c+1)\|F_{2}\|t)^{k}}{k!}}.

(50)

By Lemma 14, Eq.(40) and Eq.(50),

\|e^{At}\|\leq c+1.

(51)

Then the upper bound of the condition number $\kappa_{C}$ of $C_{m,k,p}(Ah)$ is analyzed in Lemma 8.

Lemma 8.

Consider the matrix $C_{m,k,p}(Ah)$ defined in Sect.III.3. Let $h\in\mathbbm{R^{+}}$ and satisfies $\|Ah\|\leq 1$ , $c,m,k,p\in\mathbbm{Z}^{+}$ . When $\frac{(c+1)\|F_{2}\|}{|Re(\lambda_{1})|}\leq 1$ , $k\geq 5$ and $\frac{2}{(k+1)!}m(c+1)(c+2)\leq 1$ , the condition number $\kappa_{C}$ of $C_{m,k,p}(Ah)$ satisfies

\kappa_{C}\leq 2e\sqrt{k}(m(k+1)+p)(c+2),

(52)

where $e$ is mathematical constant.

Proof.

First we analyze the upper bound of $\|C_{m,k,p}(Ah)^{-1}\|$ , we have

C_{m,k,p}(Ah)|x\rangle=|B\rangle,

(53)

where $|B\rangle=\sum_{l=0}^{d}{b_{l}|l\rangle}$ , $|l\rangle$ represents an $N$ -dimensional state. We define $|b^{l}\rangle:=b_{l}|l\rangle$ ,

\sum_{l=0}^{d}{\||b^{l}\rangle\|^{2}}=\||B\rangle\|^{2}.

(54)

For any $l\in[d]_{0}$ , we define

|x^{l}\rangle:=C_{m,k,p}(Ah)^{-1}|b^{l}\rangle=\sum_{n=0}^{d}{|x_{n}^{l}\rangle}.

(55)

We consider two cases: $0\leq l<m(k+1)$ and $m(k+1)\leq l\leq d$ .

When $0\leq l<m(k+1)$ , assuming $l=a(k+1)+b$ , $0\leq a<m$ , $0\leq b\leq k$ . Then based on definition of $x$ , we have

$\displaystyle\|x_{i,j}^{l}\rangle$	$\displaystyle=0,$	$\displaystyle 0\leq i<a,0\leq j\leq k,$	(56)
$\displaystyle\|x_{a,j}^{l}\rangle$	$\displaystyle=0,$	$\displaystyle 0\leq j<b,$
$\displaystyle\|x_{a,j}^{l}\rangle$	$\displaystyle=b!(Ah)^{j-b}/j!\|b^{l}\rangle,$	$\displaystyle b\leq j\leq k$
$\displaystyle\|x_{a+1,0}^{l}\rangle$	$\displaystyle=T_{b,k}(Ah)\|b^{l}\rangle,$
$\displaystyle\|x_{a+1,j}^{l}\rangle$	$\displaystyle=(Ah)^{j}/j!\|x_{a+1,0}^{l}\rangle,$
$\displaystyle\|x_{a+2,0}^{l}\rangle$	$\displaystyle=T_{k}(Ah)\|x_{a+1,0}^{l}\rangle=T_{k}(Ah)T_{b,k}(Ah)\|b^{l}\rangle,$
	$\displaystyle\vdots$
$\displaystyle\|x_{m,0}^{l}\rangle$	$\displaystyle=T_{k}(Ah)\|x_{m-1,0}^{l}\rangle=\left(T_{k}(Ah)\right)^{m-a-1}T_{b,k}(Ah)\|b^{l}\rangle,$
$\displaystyle\|x_{m,j}^{l}\rangle$	$\displaystyle=\|x_{m,0}\rangle=\left(T_{k}(Ah)\right)^{m-a-1}T_{b,k}(Ah)\|b^{l}\rangle,$	$\displaystyle 1\leq j\leq p,$

where $|x_{i,j}^{l}\rangle=|x_{i(k+1)+j}^{l}\rangle$ and $T_{b,k}(Ah):=\sum_{j=b}^{k}{\frac{b!(Ah)^{j-b}}{j!}}$ . By Lemma 15, we have

\|e^{Ah(m-a-1)}-T_{k}^{m-a-1}(Ah)\|\leq\frac{2}{(k+1)!}(m-a-1)(c+1)(c+2),

(57)

by Lemma 7, $\|e^{Ah(m-a-1)}\|\leq c+1$ , then we have

\|\left(T_{k}(Ah)\right)^{m-a-1}\|\leq c+2.

(58)

Therefore,

$\displaystyle\\|(Ah)^{j-b}\\|$	$\displaystyle\leq\\|(Ah)\\|^{j-b}\leq 1,\quad\quad\quad\quad\quad\quad\quad b\leq j\leq k,$	(59)
$\displaystyle\\|T_{b,k}(Ah)\\|$	$\displaystyle\leq T_{b,k}(\\|Ah\\|)\leq e^{\\|Ah\\|}\leq e,$
$\displaystyle\\|T_{k}(Ah)\\|$	$\displaystyle\leq T_{k}(\\|Ah\\|)\leq e^{\\|Ah\\|}\leq e,$
$\displaystyle\\|\left(T_{k}(Ah)\right)^{m-a-1}T_{b,k}(Ah)\\|$	$\displaystyle\leq\\|\left(T_{k}(Ah)\right)^{m-a-1}\\|\\|T_{b,k}(Ah)\\|\leq e(c+2).$

By Eq.(56) and Eq.(59), $\|x_{i,j}^{l}\|$ satisfies

$\displaystyle\\|\|x_{a,j}^{l}\rangle\\|$	$\displaystyle\leq\frac{b!}{j!}\\|(Ah)^{j-b}\\|\\|\|b^{l}\rangle\\|\leq\\|\|b^{l}\rangle\\|,$	$\displaystyle b\leq j\leq k$	(60)
$\displaystyle\\|\|x_{a+1,0}^{l}\rangle\\|$	$\displaystyle\leq\\|T_{b,k}(Ah)\\|\\|\|b^{l}\rangle\\|\leq e\\|\|b^{l}\rangle\\|,$
$\displaystyle\\|\|x_{i,j}^{l}\rangle\\|$	$\displaystyle\leq e(c+2)\\|\|b^{l}\rangle\\|,$	$\displaystyle a+1\leq i\leq m-1,0\leq j\leq k$
$\displaystyle\\|\|x_{m,j}^{l}\rangle\\|$	$\displaystyle\leq e(c+2)\\|\|b^{l}\rangle\\|,$	$\displaystyle 0\leq j\leq p$

therefore,

\|C_{m,k,p}(Ah)^{-1}|b^{l}\rangle\|^{2}\leq(m(k+1)+p)(e(c+2)))^{2}\||b^{l}\rangle\|^{2}.

(61)

When $m(k+1)\leq l\leq d$ , assuming $l=m(k+1)+b$ , $0\leq b\leq p$ , $x_{i,j}$ satisfies

$\displaystyle\|x_{i,j}^{l}\rangle$	$\displaystyle=0,$	$\displaystyle 0\leq i<m,0\leq j\leq k,$	(62)
$\displaystyle\|x_{m,j}^{l}\rangle$	$\displaystyle=0,$	$\displaystyle 0\leq j<b,$
$\displaystyle\|x_{m,j}^{l}\rangle$	$\displaystyle=\|b^{l}\rangle,$	$\displaystyle b\leq j\leq p,$

then

\|C_{m,k,p}(Ah)^{-1}|b^{l}\rangle\|^{2}\leq p\||b^{l}\rangle\|^{2}.

(63)

From Eq.(61) and Eq.(63), for any $|B\rangle$ , $\|C_{m,k,p}(Ah)^{-1}|B\rangle\|$ satisfies

	$\displaystyle\\|C_{m,k,p}(Ah)^{-1}\|B\rangle\\|^{2}\leq$	$\displaystyle(m(k+1)+p)\sum_{l=0}^{d}{\\|C_{m,k,p}(Ah)^{-1}\|b^{l}\rangle\\|^{2}}$
	$\displaystyle\leq$	$\displaystyle(m(k+1)+p)^{2}(e(c+2))^{2}\\|\|B\rangle\\|^{2},$		(64)

then

\|C_{m,k,p}(Ah)^{-1}\|\leq(m(k+1)+p)e(c+2).

(65)

From Lemma 4 in Berry et al. (2017), we have

\|C_{m,k,p}(Ah)\|\leq 2\sqrt{k}.

(66)

Therefore, by Eq.(66) and Eq.(65), we have $\kappa_{C}\leq 2e\sqrt{k}(m(k+1)+p)(c+2)$ .

VI Solution Error

In this section, we analyze the solution error of our algorithm. The error mainly comes from two aspects: (1) Homotopy perturbation method truncation error. The solution $\tilde{u}(t)$ defined in Eq.(7) is an approximate solution of Eq.(1), the error bound is determined by the truncation order $c$ , in Sect.VI.1, we analyze the convergence condition of Eq.(7) and give the error bound. (2) Linear ODEs solution error. We solve the linear ODEs defined in Eq.(8) with the quantum algorithm proposed in Berry et al. (2017). This algorithm also generates intermediate error, we analyze the error bound in Sect.VI.2.

VI.1 Homotopy perturbation method truncation error

We first analyze homotopy perturbation method truncation error, the result is shown in Lemma 9.

Lemma 9.

Consider the nonlinear ODEs defined in Eq.(1), the solution of Eq.(1) obtained by homotopy perturbation method is written as $\tilde{u}(t)=\sum_{i=0}^{c}{\nu_{i}(t)}$ , $u(t)$ represents the exact solution of Eq.(1). When $K<1$ and $c>\log_{1/K}{\frac{1}{\epsilon(1-K)}}$ , $\tilde{u}(t)$ satisfies

\|u(t)-\tilde{u}(t)\|\leq\epsilon.

(67)

Proof.

$u(t)$ can be represented as $u(t)=\sum_{i=0}^{\infty}{\nu_{i}(t)}$ , Eq.(67) is transformed into

\|\sum_{j=c+1}^{\infty}{\nu_{j}(t)}\|<\epsilon.

(68)

To prove Eq.(68), we analyze the upper bound of $\|\nu_{i}\|$ defined in Eq.(III.1). $\nu_{0}(t)=e^{F_{1}t}u_{in}$ , we have

\|\nu_{0}(t)\|\leq\|e^{F_{1}t}\|\times\|u_{in}\|\leq\|u_{in}\|.

(69)

We define $K_{1}=\frac{\|u_{in}\|\|F_{2}\|}{|Re(\lambda_{1})|}$ and assume when $j\leq i$ , $\|\nu_{j}\|$ satisfies

\|\nu_{j}\|\leq\alpha_{j}K_{1}^{j}\|u_{in}\|,\ j=0,1,\dots,i.

(70)

Then $\|\nu_{i+1}(t)\|$ satisfies

$\displaystyle\\|\nu_{i+1}(t)\\|$	$\displaystyle=\\|\int_{0}^{t}{e^{F_{1}(t-\tau)}F_{2}\sum_{j=0}^{i}{\nu_{j}(\tau)\otimes\nu_{i-j}(\tau)}d\tau}\\|$
	$\displaystyle\leq\int_{0}^{t}{\\|e^{F_{1}(t-\tau)}F_{2}\\|\times(\sum_{j=0}^{i}{\alpha_{j}\alpha_{i-j}})K_{1}^{i}\\|u_{in}\\|^{2}d\tau}$
	$\displaystyle\leq(\sum_{j=0}^{i}{\alpha_{j}\alpha_{i-j}})K_{1}^{i}\\|u_{in}\\|^{2}\\|F_{2}\\|\int_{0}^{t}{\\|e^{Re(\lambda_{1})(t-\tau)}\\|d\tau}$
	$\displaystyle\leq(\sum_{j=0}^{i}{\alpha_{j}\alpha_{i-j}})K_{1}^{i}\\|u_{in}\\|^{2}\\|F_{2}\\|\frac{1-e^{Re(\lambda_{1})t}}{\|Re(\lambda_{1})\|}$
	$\displaystyle\leq(\sum_{j=0}^{i}{\alpha_{j}\alpha_{i-j}})K_{1}^{i+1}\\|u_{in}\\|.$	(71)

By Eq.(70), Eq.(Proof), $\alpha_{i}$ can be defined as

\alpha_{i+1}=\sum_{j=0}^{i}{\alpha_{j}\alpha_{i-j}},\ \alpha_{0}=1.

(72)

Eq.(72) is the catalan sequenceKoshy (2008) and satisfies

\alpha_{i}=\frac{1}{i+1}\binom{2i}{i}\approx\frac{4^{i}}{i^{3/2}\sqrt{\pi}}<4^{i}.

(73)

Combining Eq.(3), Eq.(70), Eq.(Proof), Eq.(72) and Eq.(73), for any $i\in\mathbbm{N}$ , $\|\nu_{i}(t)\|$ has the upper bound

\|\nu_{i}(t)\|<(4K_{1})^{i}\|u_{in}\|\leq K^{i+1}.

(74)

Substituting Eq.(74) into Eq.(68), we have

\sum_{i=c+1}^{\infty}{K^{i+1}}<\epsilon.

(75)

Therefore, when $K<1$ and $c>\log_{1/K}{\frac{1}{\epsilon(1-K)}}$ , we have $\|u(t)-\tilde{u}(t)\|\leq\epsilon$ .

VI.2 Linear ODEs solution error

Then we analyze the error of solving the linear ODEs defined in Eq.(8), we have a similar conclusion with Theorem $\bm{6}$ in Berry et al. (2017), the difference comes from the upper bound of $\|e^{Aj}\|$ or $\|T_{k}^{j}(A)\|$ in our work is different from their work. Our result is shown in Lemma 10.

Lemma 10.

Consider the linear ODEs defined in Eq.(8) and the system of linear equations defined in Eq.(22). Let $h\in\mathbbm{R}^{+}$ satisfies $\|Ah\|\leq 1$ . When $\frac{2}{(k+1)!}m(c+1)(c+2)\leq 1$ , we have

\||y(jh)\rangle-|x_{j,0}\rangle\|\leq\frac{2j(c+1)(c+2)\|y_{in}\|}{(k+1)!},j=0,1,\dots,m.

(76)

Proof.

The exact solution of Eq.(8) at $t=jh$ is $|y(jh)\rangle$ , it is written as

|y(jh)\rangle=e^{Ajh}|y(0)\rangle.

(77)

The solution of Eq.(22) $|x_{j,0}\rangle$ satisfies

|x_{j,0}\rangle=T_{k}^{j}(Ah)|x_{0,0}\rangle,

(78)

where $T_{k}(Ah)=\sum_{l=0}^{k}{\frac{(Ah)^{l}}{l!}}$ , the initial value satisfies $|y(0)\rangle=|x_{0,0}\rangle=|y_{in}\rangle$ . By Lemma 7 and Lemma 15, $\||y(jh)\rangle-|x_{j,0}\rangle\|$ satisfies

\displaystyle\||y(jh)\rangle-|x_{j,0}\rangle\|\leq\|e^{Ajh}-T_{k}^{j}(Ah)\|\|y_{in}\|\leq\frac{2j(c+1)(c+2)\|y_{in}\|}{(k+1)!},\ j=0,1,\dots,m.

(79)

VII Success Probability

As introduced in Sect.III.4, there are two probabilistic steps in our method. This section gives a lower bound of the success probability of these two steps. The results are shown in Lemma 11 and Lemma 12.

Firstly, we analyze the success probability of getting $|x_{m,i}\rangle,i=0,1,\dots,p$ when measuring $|x\rangle$ . By setting appropriate conditions, Lemma 11 gives the same conclusion as Theorem $\bm{7}$ in Berry et al. (2017).

Lemma 11.

Consider the system of linear equations defined in Eq.(22). Let $h\in\mathbbm{R}^{+}$ satisfies $\|Ah\|\leq 1$ . $g=\max_{t\in[0,mh]}\||y(t)\rangle\|/\||y(mh)\rangle\|$ , $m,k,p\in\mathbbm{Z}^{+}$ . When $(k+1)!\geq 50m(c+1)(c+2)g$ , we have

\frac{\||x_{m,0}\rangle\|^{2}}{\||x\rangle\|^{2}}\geq\frac{1}{p+77mg^{2}}.

(80)

Proof.

As introduced before, $|x_{m,j}\rangle=|x_{m,0}\rangle$ for $j\in[p]_{0}$ , we define

|x_{good}\rangle:=\sum_{j=0}^{p}|m(k+1)+j\rangle|x_{m,j}\rangle

(81)

and

|x_{bad}\rangle:=\sum_{i=0}^{m-1}{\sum_{j=0}^{k}{|i(k+1)+j\rangle|x_{i,j}\rangle}}.

(82)

We see $|x\rangle=|x_{good}\rangle+|x_{bad}\rangle$ and $\langle x_{good}|x_{bad}\rangle=0$ , then

	$\displaystyle\\|\|x\rangle\\|^{2}=$	$\displaystyle\\|\|x_{good}\rangle\\|^{2}+\\|\|x_{bad}\rangle\\|^{2}$
	$\displaystyle=$	$\displaystyle(p+1)\\|\|x_{m,0}\rangle\\|^{2}+\\|\|x_{bad}\rangle\\|^{2}.$		(83)

Next we give a lower bound of $\||x_{m,0}\rangle\|$ and an upper bound of $\||x_{bad}\rangle\|$ . We define $q=\||y(mh)\rangle\|$ , by Lemma 10 and $(k+1)!\geq 50m(c+1)(c+2)g$ , we have

\||x_{i,0}\rangle-|y(ih)\rangle\|\leq 0.04q,\ 0\leq l\leq m.

(84)

By the definition of $g$ , $\||y(ih)\rangle\|\leq gq$ for any $i\in[m-1]_{0}$ , then we have

\||x_{i,0}\rangle\|\leq(g+0.04)q\leq 1.04gq,\ 0\leq i\leq m-1,

(85)

and

0.96q\leq\||x_{m,0}\rangle\|\leq 1.04q.

(86)

For any $i\in[m]_{0}$ , $|x_{i,j}\rangle=(Ah/j)|x_{i,j-1}\rangle$ , we have

|x_{i,j}\rangle=\frac{(Ah)^{j-1}}{j!}|x_{i,1}\rangle,\ 2\leq j\leq k.

(87)

We have $\|Ah\|\leq 1$ , therefore,

\||x_{i,j}\rangle\|\leq\frac{\||x_{i,1}\rangle\|}{j!},\ 2\leq j\leq k.

(88)

Next, based on $|x_{i+1,0}\rangle=|x_{i,0}\rangle+\sum_{j=1}^{k}{|x_{i,j}\rangle}$ , we have

$\displaystyle 2.08gq\geq$	$\displaystyle\\|\|x_{i+1,0}\rangle\\|+\\|\|x_{i,0}\rangle\\|$
$\displaystyle\geq$	$\displaystyle\\|\|x_{i+1,0}\rangle-\|x_{i,0}\rangle\\|$
$\displaystyle\geq$	$\displaystyle\\|\|x_{i,1}\rangle\\|-\sum_{j=2}^{k}{\\|\|x_{i,j}\rangle\\|}$
$\displaystyle\geq$	$\displaystyle\left(1-\sum_{j=2}^{k}{\frac{1}{j!}}\right)\\|\|x_{i,1}\rangle\\|$
$\displaystyle\geq$	$\displaystyle(3-e)\\|\|x_{i,1}\rangle\\|,$	(89)

then

\||x_{i,1}\rangle\|\leq\frac{2.08gq}{3-e},\ 0\leq i\leq m-1.

(90)

and

\||x_{i,j}\rangle\|\leq\frac{2.08gq}{j!(3-e)},\ 0\leq i\leq m-1,\ 1\leq j\leq k.

(91)

$\||x_{bad}\rangle\|$ satisfies

$\displaystyle\\|\|x_{bad}\rangle\\|^{2}$	$\displaystyle=\sum_{i=0}^{m-1}{\\|\|x_{i,0}\rangle\\|^{2}}+\sum_{i=0}^{m-1}{\sum_{j=1}^{k}{\\|\|x_{i,j}\rangle\\|^{2}}}$
	$\displaystyle\leq 1.04^{2}mg^{2}q^{2}+m\sum_{j=1}^{k}{\frac{(2.08gq)^{2}}{(j!)^{2}(3-e)^{2}}}$
	$\displaystyle\leq 70.9mg^{2}q^{2}.$	(92)

The last step of Eq.(Proof) is derived from the inequality $\sum_{j=1}^{k}{\frac{1}{(j!)^{2}}}\leq I_{0}(2)-1<1.28$ , where $I_{0}(2)<2.28$ a modified Bessel function of the first kindBerry et al. (2017). Combining Eq.(86) and Eq.(Proof), we have

	$\displaystyle\frac{\\|\|x_{m,0}\rangle\\|^{2}}{\\|\|x\rangle\\|^{2}}$	$\displaystyle\geq\frac{(0.96q)^{2}}{p(0.96q)^{2}+70.9mg^{2}q^{2}}$
		$\displaystyle\geq\frac{1}{p+77mg^{2}}.$		(93)

Secondly, we analyze the success probability of the second probabilistic step. After the first measurement, the desired state is the state $|y(T)\rangle$ defined in Eq.(19). Then we measure the first qubit register of $|y(T)\rangle$ , if the result is $|0,0\rangle$ , we have a state $\epsilon$ -close to $|u(T)/\|u(T)\|\rangle$ in the second qubit register of $|y(T)\rangle$ . The lower bound of the success probability in this measurement is analyzed in Lemma 12.

Lemma 12.

Let $T\in\mathbbm{R}^{+}$ , $\eta^{{}^{\prime}}=K/\|\tilde{u}(T)\|$ . When $K<\sqrt{2}/2$ , we have

\frac{\||y_{0}(T)\rangle\|^{2}}{\||y(T)\rangle\|^{2}}\geq\frac{1-2K^{2}}{1-2K^{2}+2(\eta^{{}^{\prime}})^{2}}.

(94)

Proof.

We rearrange the components of $\vec{y}$ ,

\vec{y}=[\vec{y}_{0}^{{}^{\prime}},\vec{y}_{1}^{{}^{\prime}},\dots,\vec{y}_{c}^{{}^{\prime}}],

(95)

where $\vec{y}_{0}^{{}^{\prime}}=\vec{y}_{0}=\sum_{j=0}^{c}{\nu_{j}}$ , when $i\geq 1$ , the component of $\vec{y}_{i}^{{}^{\prime}}$ is represented as $\nu_{a_{i,0}^{{}^{\prime}}}\otimes\nu_{a_{i,1}^{{}^{\prime}}}\otimes\dots\otimes\nu_{a_{i,k}^{{}^{\prime}}}$ , $a_{i,j}^{{}^{\prime}}$ satisfies

k\geq 1,\ a_{i,j}^{{}^{\prime}}\geq 0,\ \sum_{j=0}^{k}{(a_{i,j}^{{}^{\prime}}+1)}=i+1.

(96)

By Eq.(96), the number of elements in $\vec{y}_{i}^{{}^{\prime}}$ is $2^{i}-1$ . By Eq.(74) and Eq.(96), for any $j\in[2^{i}-1]_{0}$ , we have $\|\vec{y}_{i,j}^{{}^{\prime}}\|\leq K^{i+1}$ . Therefore,

\|\vec{y}_{i}^{{}^{\prime}}\|^{2}\leq(2^{i}-1)K^{2(i+1)}<(2K^{2})^{i}.

(97)

By Eq.(97) and $\|\vec{y}_{0}\|=\|\vec{y}_{0}^{{}^{\prime}}\|=\|\tilde{u}(T)\|=K/\eta^{{}^{\prime}}$ , we have

$\displaystyle\frac{\\|\|y_{0}\rangle\\|^{2}}{\\|\|y\rangle\\|^{2}}=$	$\displaystyle\frac{(K/\eta^{{}^{\prime}})^{2}}{(K/\eta^{{}^{\prime}})^{2}+\sum_{i=1}^{c}{(2K^{2})^{i}}}$	(98)
$\displaystyle\geq$	$\displaystyle\frac{(K/\eta^{{}^{\prime}})^{2}}{(K/\eta^{{}^{\prime}})^{2}+\frac{2K^{2}}{1-2K^{2}}}$	(99)
$\displaystyle=$	$\displaystyle\frac{1-2K^{2}}{1-2K^{2}+2(\eta^{{}^{\prime}})^{2}}.$	(100)

VIII Main Result

In this section, we give the main result of our work.

Theorem 1.

Given $n$ -dimensional nonlinear dissipative ODEs $\frac{du}{dt}=F_{1}u+F_{2}u^{\otimes 2}$ defined in Eq.(1) and construct the linear ODEs $\frac{d\vec{y}}{dt}=A\vec{y}$ defined in Eq.(8). Let $T>0$ , $g=\max_{t\in[0,T]}\||y(t)\rangle\|/\||y(T)\rangle\|$ , $\eta=\|u_{in}\|/\|u(T)\|$ , $c=\left\lceil\log_{1/K}{\frac{4\|u_{in}\|}{(1-K)\epsilon\eta}}\right\rceil$ . When $K<\sqrt{2}/2$ and $\frac{(c+1)\|F_{2}\|}{|Re(\lambda_{1})|}\leq 1$ , there exists a quantum algorithm to produce $|u_{out}(T)\rangle$ which satisfies $\|u_{out}(T)-u(T)/\|u(T)\|\|\leq\epsilon$ with $\Omega(1)$ success probability. The query complexity of the algorithm for $O_{F1}$ , $O_{F2}$ , $O_{u}$ is

O\left(\frac{g\eta sT(\|F_{1}\|+\|F_{2}\|)}{\sqrt{1-2K^{2}}\|u_{in}\|}{\rm poly}\left(\log\left(\frac{g\eta sT\|F_{1}\|\|F_{2}\|}{\epsilon(1-2K^{2})\|u_{in}\|}\right)\right)\right).

(101)

The gate complexity of this algorithm is larger than its query complexity by a factor of

O\left({\rm poly}\left(\log\left(\frac{ng\eta sT\|F_{1}\|\|F_{2}\|}{\epsilon(1-2K^{2})\|u_{in}\|}\right)\right)\right).

(102)

Proof.

Whole process. First we show the whole process of our algorithm. We define $\eta^{{}^{\prime}}=\eta K/\|u_{in}\|$ and set

\epsilon\leq 0.1\sqrt{1-2K^{2}}/\eta^{{}^{\prime}},\ \epsilon_{1}=\frac{\epsilon K}{4\eta^{{}^{\prime}}},\ \delta=\epsilon\frac{\sqrt{1-2K^{2}}}{30\sqrt{78m}g\eta^{{}^{\prime}}},

(103)

and construct the $N$ -dimensional linear ODEs defined in Eq.(8)

d\vec{y}/dt=A\vec{y},\ \vec{y}(0)=y_{in}.

(104)

We also set $h=T/\lceil T\|A\|\rceil$ , $m=p=T/h=\lceil T\|A\|\rceil$ , $k=\lfloor\frac{2\log(\Omega)}{\log(\log(\Omega))}\rfloor$ , where $\Omega=50m(c+1)(c+2)g/\delta$ , so $k$ satisfies $(k+1)!\geq\Omega$ . Then we construct the linear system defined in Eq.(22) and solve the linear system with the algorithm proposed in Childs et al. (2017). The normalized solution of Eq.(22) is represented as

|x\rangle=\sum_{l=0}^{d}{|l\rangle|x_{l}\rangle},

(105)

where $d=m(k+1)+p$ . $|x\rangle$ can also be represented as

|\bar{x}\rangle=\sum_{l=0}^{d}{\alpha_{l}|l\rangle|\bar{x}_{l}\rangle},

(106)

where $|\bar{x}_{l}\rangle$ is a normalized state. Assuming the output state of quantum linear system algorithmChilds et al. (2017) is

|\bar{x}^{{}^{\prime}}\rangle=\sum_{l=0}^{d}{\alpha_{l}^{{}^{\prime}}|l\rangle|\bar{x}_{l}^{{}^{\prime}}\rangle}.

(107)

By Theorem $\bm{5}$ in Childs et al. (2017), we make $|\bar{x}^{{}^{\prime}}\rangle$ satisfies

\||\bar{x}\rangle-|\bar{x}^{{}^{\prime}}\rangle\|\leq\delta.

(108)

Then we execute measurement step discussed in Sect.III.4 and get a state $\epsilon$ -close to $|u(T)/\|u(T)\|\rangle$ with success probability $O((1-2K^{2})/(g\eta^{{}^{\prime}})^{2})$ . We can amplify the success probability to $\Omega(1)$ with quantum amplitude amplification algorithmBrassard et al. (2002) by running quantum linear system algorithm $O(g\eta^{{}^{\prime}}/\sqrt{1-2K^{2}})$ times.

Proof of correctness. Then we analyze the error bound and the impact of error on success probability.

Assuming $u(T)$ represents the exact solution. We define $|\tilde{u}(T)\rangle$ as

|\tilde{u}(T)\rangle=|\sum_{i=0}^{c}{\nu_{i}(T)}\rangle=|y_{0,0}(T)\rangle.

(109)

By Lemma 9 and our choice of parameter $c$ , we have

\||u(T)\rangle-|\tilde{u}(T)\rangle\|\leq\epsilon_{1}.

(110)

Then by Lemma 16 and $\||u(T)\rangle\|=K/\eta^{{}^{\prime}}$ , we have

\||\bar{u}(T)\rangle-|\bar{\vec{y}}_{0,0}(T)\rangle\|\leq\frac{2\eta^{{}^{\prime}}\epsilon_{1}}{K}=\epsilon/2.

(111)

Let $S:=\{m(k+1),m(k+1)+1,\dots,m(k+1)+p\}$ . By Lemma 10 and our choice of parameter $k$ , for any $l\in S$ , we have

\||x_{l}\rangle-|y(T)\rangle\|\leq\delta\||y(T)\rangle\|.

(112)

By Eq.(112) and Lemma 16, we have

\||\bar{x}_{l}\rangle-|\bar{y}(T)\rangle\|\leq 2\delta.

(113)

By Lemma 11, for any $l\in S$ , we have

\alpha_{l}=\frac{\||x_{l}\rangle\|}{\||x\rangle\|}\geq\frac{1}{\sqrt{78m}g}\geq\delta.

(114)

By Eq.(108), Eq.(114) and Lemma 17, we have

\||\bar{x}_{l}\rangle-|\bar{x}^{{}^{\prime}}_{l}\rangle\|\leq\frac{2\delta}{\alpha_{l}-\delta}.

(115)

Then, combine Eq.(113) and Eq.(115), we have

	$\displaystyle\\|\|\bar{x}^{{}^{\prime}}_{l}\rangle-\|\bar{y}(T)\rangle\\|$	$\displaystyle\leq\\|\|\bar{x}_{l}^{{}^{\prime}}\rangle-\|\bar{x}_{l}\rangle\\|+\\|\|\bar{x}_{l}\rangle-\|\bar{y}(T)\rangle\\|$
		$\displaystyle\leq 2\delta\left(1+\frac{1}{\alpha_{l}-\delta}\right).$		(116)

We default $\alpha_{l}$ is small enough, such as $\alpha_{l}<0.5$ , then by Eq.(114) and $\delta=\epsilon\frac{\sqrt{1-2K^{2}}}{30\sqrt{78m}g\eta^{{}^{\prime}}}$ , we have

2\delta(1+\frac{1}{\alpha_{l}-\delta})\leq 2\delta\frac{\sqrt{3}}{\alpha_{l}}\leq\frac{\epsilon\sqrt{1-2K^{2}}}{5\sqrt{3}\eta^{{}^{\prime}}}.

(117)

$|y(T)\rangle$ and $|\bar{x}_{l}^{{}^{\prime}}\rangle$ are also written as

|\bar{y}(T)\rangle=\sum_{w=0}^{c}{\chi_{w}|\bar{y}_{w}(T)\rangle}

(118)

and

|\bar{x}_{l}^{{}^{\prime}}\rangle=\sum_{w=0}^{c}{\chi_{w}^{{}^{\prime}}|\bar{x}_{l,w}^{{}^{\prime}}\rangle}.

(119)

By Lemma 12, we have

\chi_{0}=\frac{\||y_{0}(T)\rangle\|}{\||y(T)\rangle\|}\geq\sqrt{\frac{1-2K^{2}}{(1-2K^{2})+2(\eta^{{}^{\prime}})^{2}}}\geq\frac{\sqrt{1-2K^{2}}}{\sqrt{3}\eta^{{}^{\prime}}}.

(120)

By Eq.(Proof), Eq.(117), Eq.(120) and Lemma 17, we have

\||\bar{x}_{l,0}^{{}^{\prime}}\rangle-|\bar{y}_{0,0}(T)\rangle\|\leq\frac{2(\||\bar{x}^{{}^{\prime}}_{l}\rangle-|\bar{y}(T)\rangle\|)}{\chi_{0}-(\||\bar{x}^{{}^{\prime}}_{l}\rangle-|\bar{y}(T)\rangle\|)}\leq\frac{2\epsilon}{5-\epsilon}\leq\epsilon/2.

(121)

We notice $|\bar{x}_{l,0}^{{}^{\prime}}\rangle$ is the output state of our algorithm, we have $|u_{out}(T)\rangle=|\bar{x}_{l,0}^{{}^{\prime}}\rangle$ . Combining Eq.(111) and Eq.(121), we have

\||\bar{u}(T)\rangle-|u_{out}(T)\rangle\|\leq\||\bar{u}(T)\rangle-|\bar{y}_{0,0}(T)\rangle\|+\||\bar{x}_{l,0}^{{}^{\prime}}\rangle-|\bar{y}_{0,0}(T)\rangle\|\leq\epsilon.

(122)

On the other hand, caused by solution error, the success probability also changes. By Lemma 18, Eq.(120) and $\epsilon\leq 0.1\sqrt{1-2K^{2}}/\eta^{{}^{\prime}}$ , we have

\chi_{0}^{{}^{\prime}}\geq\chi_{0}-\frac{1}{\sqrt{3}\eta^{{}^{\prime}}}\sqrt{1-2K^{2}}\times\frac{1}{5}\epsilon\geq\chi_{0}-\epsilon/2\geq\frac{\sqrt{1-2K^{2}}}{2\eta^{{}^{\prime}}}

(123)

and

\alpha_{l}^{{}^{\prime}}\geq\alpha_{l}-\delta\geq\left(\frac{1}{\sqrt{78}}-\frac{1}{30\sqrt{78}}\right)\frac{1}{\sqrt{m}g}\geq\frac{1}{10\sqrt{m}g},

(124)

then

\sum_{l\in S}{|\alpha_{l}^{{}^{\prime}}|^{2}}\geq\frac{p}{100mg^{2}}=\frac{1}{100g^{2}},

(125)

Therefore, the success probability of our algorithm is

p=(\chi_{0}^{{}^{\prime}})^{2}\times\left(\sum_{l\in S}{|\alpha_{l}^{{}^{\prime}}|^{2}}\right)\geq\frac{1-2K^{2}}{400(g\eta^{{}^{\prime}})^{2}}.

(126)

Complexity analysis. Finally, we analyze the complexity of our algorithm.

By Lemma 2 and Lemma 3, the sparsity of $A$ is $c^{2}s$ and $\|A\|\leq(c+1)(\|F_{1}\|+\|F_{2}\|)$ . The sparsity $s_{C}$ of $C_{m,k,p}(Ah)$ satisfies

s_{C}<k+c^{2}s.

(127)

By Lemma 8, the condition number $\kappa_{C}$ of $C_{m,k,p}(Ah)$ satisfies

\kappa_{C}\leq 2e\sqrt{k}((m(k+1)+p))(c+2).

(128)

By Theorem 5 in Childs et al. (2017) and Eq.(108), the query complexity of quantum linear system algorithm to oracle of $C_{m,k,p}(Ah)$ and $|0\rangle|y_{in}\rangle$ is

O(s_{C}\kappa_{C}{\rm poly}({\rm log}(s_{C}\kappa_{C}/\delta))).

(129)

By the definition of $C_{m,k,p}(Ah)$ , the oracle of $C_{m,k,p}(Ah)$ can be constructed by querying oracle $O_{A}$ once, by Lemma 6, $O_{A}$ is constructed by querying $O_{F1}$ $O(c)$ times and $O_{F2}$ $O(1)$ times. By Lemma 1, $|0\rangle|y_{in}\rangle$ can be prepared by querying $O_{u}$ $O(c)$ times.

Substituting Eq.(127), Eq.(128) into Eq.(129) and considering the choice of all parameters, the query complexity of solving Eq.(22) for $O_{F1}$ , $O_{F2}$ and $O_{u}$ is

O\left(sT(\|F_{1}\|+\|F_{2}\|){\rm poly}\left(\log\left(\frac{g\eta sT\|F_{1}\|\|F_{2}\|}{\epsilon(1-2K^{2})\|u_{in}\|}\right)\right)\right).

(130)

Using amplitude amplification algorithmBrassard et al. (2002), we repeat the above process $O(1/\sqrt{p})$ times and get $|u_{out}(T)\rangle=|\bar{x}_{l,0}^{{}^{\prime}}(T)\rangle$ which satisfies $\||u_{out}(T)\rangle-u(T)/\|u(T)\|\|\leq\epsilon$ with $\Omega(1)$ success probability.

The query complexity of the whole process for $O_{F1}$ , $O_{F2}$ and $O_{u}$ is

O\left(\frac{g\eta sT(\|F_{1}\|+\|F_{2}\|)}{\sqrt{1-2K^{2}}\|u_{in}\|}{\rm poly}\left(\log\left(\frac{g\eta sT\|F_{1}\|\|F_{2}\|}{\epsilon(1-2K^{2})\|u_{in}\|}\right)\right)\right).

(131)

The gate complexity of this algorithm is larger than its query complexity by a factor of

O\left({\rm poly}\left(\log\left(\frac{ng\eta sT\|F_{1}\|\|F_{2}\|}{\epsilon(1-2K^{2})\|u_{in}\|}\right)\right)\right).

(132)

IX Conclusion and Discussion

In this paper, we presented a quantum homotopy perturbation method for solving nonlinear dissipative ODEs. The gate complexity of our algorithm is $O(g\eta T{\rm poly}(\log(nT/\epsilon)))$ . The complexity of the optimal classical algorithm for solving Eq.(1) is at least linear with $n$ , the complexity of the algorithm proposed in Liu et al. (2021) is linear with $1/\epsilon$ , so our algorithm provides exponential improvement over the best classical algorithms or previous quantum algorithms in $n$ or $\epsilon$ . $\eta$ and $g$ also affect the complexity of our algorithm, $\eta$ measures the decay of $u(T)$ and increases exponentially as $T$ increases, $g$ measures the decay of $\vec{y}(T)$ defined in Eq.(8) and also increases exponentially as $T$ increases. Our algorithm is effective when $T$ is relatively small which makes $\eta$ and $g$ small enough. $\eta$ and $g$ are also affected by $F_{2}$ , when $T$ is relatively small, the trend of $\eta$ and $g$ increasing exponentially with $T$ may not be obvious due to the influence of $F_{2}$ , this case makes our algorithm perform better. Our algorithm has the potential to accelerate the solution of various nonlinear equations, and can be applied to nonlinear problems in various fields, such as fluid dynamics, biology, finance, etc, thereby accelerating the research progress of nonlinear science.

Our algorithm only discusses time-independent homogeneous quadratic nonlinear ODEs. When solving time-dependent nonlinear ODEs, the algorithm proposed in Berry et al. (2017) is not suitable, an alternative way is to use the algorithm proposed in Berry (2014) to solve the linear ODEs, then the dependence of complexity on error $\epsilon$ becomes $O(\rm{poly}(1/\epsilon))$ . Is it possible to optimize the complexity of time-dependent quadratic nonlinear ODEs to $O(\rm poly(\log(1/\epsilon)))$ is an open question.

On the other hand, homotopy analysis methodShijun (1998) and its derivativesIlhan et al. (2021); Veeresha et al. (2021a, b) are similar to homotopy perturbation method. Whether we can use quantum computing to accelerate the execution process of homotopy analysis method and thus construct a quantum homotopy analysis method is also a question to be investigated further.

Furthermore, how to induce nonlinearity in quantum computing is a basic problem when solving nonlinear equations with quantum algorithm. A common method is producing multiple copies of the original system, some nonlinear quantum algorithms contain copy processLeyton and Osborne (2008); Lubasch et al. (2020); Lloyd et al. (2020). In Joseph (2020), a linearization technique of nonlinear classical dynamics based on Koopman-von Neumann method is proposed. Engel et al. (2021) summarizes three classical linear embedding techniques, including Carleman embedding(Carleman linearization is also called Carleman embedding)Carleman et al. (1932); Kowalski and Steeb (1991), coherent states embeddingKowalski and Steeb (1991); Kowalski (1994) and position-space embeddingKoopman (1931), and then puts forward the prospects of these linear embedding techniques to construct effective quantum algorithms. An open question is whether there are other ways to induce nonlinearity in quantum computing.

Acknowledgements

This work was supported by the National Key Research and Development Program of China (Grant No. 2016YFA0301700), the National Natural Science Foundation of China (Grants Nos. 11625419), the Strategic Priority Research Program of the Chinese Academy of Sciences (Grant No. XDB24030600), and the Anhui Initiative in Quantum Information Technologies (Grants No. AHY080000).

Appendix

In this appendix, we give some lemmas used in proving some conclusions of our work. Lemma 16, Lemma 17 and Lemma 18 are given in Berry et al. (2017), we just list them again.

Lemma 13.

Let $\gamma\geq 1$ , $m\in N^{+}$ , when $t\geq 0$ , we have

\sum_{j=0}^{m-1}{\frac{t^{j}}{j!}{e^{-\gamma t}}}\leq m.

(133)

Proof.

We consider two cases: (1) $0<t\leq m-1$ ; (2) $t\geq m$ .

When $0<t\leq m-1$ , we can find $i\in[m-1]_{0}$ which satisfies $t\in(i,i+1]$ , then we have

\sum_{j=0}^{m-1}{\frac{t^{j}}{j!}{e^{-\gamma t}}}<m\frac{t^{i}}{i!}e^{-\gamma t}.

(134)

We define

g(t)=m\frac{t^{i}}{i!}e^{-\gamma t},\ 0<t\leq m-1.

(135)

It is obvious that $g(t)\leq g(\frac{i}{\gamma})$ for $0<t\leq m-1$ . Using Stirling’s formula $i!\approx\sqrt{2\pi i}(\frac{i}{e})^{i}$ , we have

g(t)\leq m\frac{\left(\frac{i}{\gamma}\right)^{i}}{\sqrt{2\pi i}(\frac{i}{e})^{i}}e^{-i}\leq m\left(\frac{1}{\gamma}\right)^{i}\leq m.

(136)

Therefore,

\sum_{j=0}^{m-1}{\frac{t^{j}}{j!}{e^{-\gamma t}}}\leq m.

(137)

When $t\geq m$ , we have

\sum_{j=0}^{m-1}{\frac{t^{j}}{j!}{e^{-\gamma t}}}<m\frac{t^{m-1}}{(m-1)!}e^{-\gamma t}.

(138)

Similar with the case $0<t\leq m-1$ , we also have

m\frac{t^{m-1}}{(m-1)!}e^{-\gamma t}\leq m\frac{\left(\frac{m-1}{\gamma}\right)^{m-1}}{\sqrt{2\pi(m-1)}\left(\frac{m-1}{e}\right)^{m-1}}e^{-(m-1)}\leq m\left(\frac{1}{\gamma}\right)^{m-1}\leq m.

(139)

Therefore, for any $t\geq 0$ , we have $\sum_{j=0}^{m-1}{\frac{t^{j}}{j!}{e^{-\gamma t}}}\leq m$ .

Lemma 14.

Let $\gamma,\beta\in\mathbbm{R}^{+}$ , $m\in N^{+}$ , when $t\geq 0$ and $\gamma/\beta\geq 1$ , we have

\sum_{j=0}^{m-1}{\frac{(\beta t)^{j}}{j!}{e^{-\gamma t}}}\leq m.

(140)

Proof.

We define $t^{{}^{\prime}}=\beta t$ , $\gamma^{{}^{\prime}}=\gamma/\beta$ , then Eq.(140) is written as

\sum_{j=0}^{m-1}{\frac{(t^{{}^{\prime}})^{j}}{j!}{e^{-\gamma^{{}^{\prime}}t^{{}^{\prime}}}}}\leq m.

(141)

By Lemma 13, we directly obtain Eq.(141).

Lemma 15.

Given a matrix $M$ satisfies $\|M\|\leq 1$ and $\|e^{Mt}\|\leq\Delta$ for any $t\geq 0$ . Let $k,m\in N^{+}$ and satisfy $\frac{2}{(k+1)!}m\Delta(\Delta+1)\leq 1$ , $T_{k}(M)=\sum_{l=0}^{k}{\frac{(Ah)^{l}}{l!}}$ . Then for any $l\in[m]_{0}$ , we have

\|e^{Ml}-T_{k}^{l}(M)\|\leq\frac{2l\Delta(\Delta+1)}{(k+1)!}.

(142)

Proof.

When $l=0$ , $\|e^{Ml}-T_{k}^{l}(M)\|=0$ . When $l=1$ ,

\|e^{M}-T_{k}(M)\|=\|\sum_{j=k+1}^{\infty}{\frac{M^{j}}{j!}}\|\leq\sum_{j=k+1}^{\infty}{\frac{\|M\|^{j}}{j!}}\leq\sum_{j=k+1}^{\infty}{\frac{1}{j!}}\leq\frac{2}{(k+1)!}\leq\frac{2}{(k+1)!}\Delta(\Delta+1).

(143)

Assuming when $l\leq l^{{}^{\prime}}$ , we have

\|e^{Ml}-T_{k}^{l}(M)\|\leq\frac{2}{(k+1)!}l\Delta(\Delta+1),

(144)

then by $\|e^{Ml}\|\leq\Delta$ , we have

\|T_{k}^{l}(M)\|\leq\Delta+\frac{2}{(k+1)!}l\Delta(\Delta+1)\leq\Delta+1.

(145)

When $l=l^{{}^{\prime}}+1$ ,

$\displaystyle\\|e^{M(l^{{}^{\prime}}+1)}-T_{k}^{l^{{}^{\prime}}+1}(M)\\|=$	$\displaystyle\\|(e^{M}-T_{k}(M))(\sum_{j=0}^{l^{{}^{\prime}}}{e^{Mj}T_{k}^{l-j}(M)})\\|$
$\displaystyle\leq$	$\displaystyle\frac{2\Delta}{(k+1)!}\times(\sum_{j=0}^{l^{{}^{\prime}}}{\\|T_{k}^{l-j}(M)\\|})$
$\displaystyle\leq$	$\displaystyle\frac{2}{(k+1)!}(l^{{}^{\prime}}+1)\Delta(\Delta+1).$	(146)

Therefore, for any $l\in[m]_{0}$ , we have $\|e^{Ml}-T_{k}^{l}(M)\|\leq\frac{2}{(k+1)!}l\Delta(\Delta+1)$ .

Lemma 16.

(Berry et al. (2017)). Let $|\psi\rangle$ and $|\varphi\rangle$ be two vectors such that $\||\psi\rangle\|\geq\alpha>0$ and $\||\psi\rangle-|\varphi\rangle\|\leq\beta$ . Then

\left\|\frac{|\psi\rangle}{\||\psi\rangle\|}-\frac{|\varphi\rangle}{\||\varphi\rangle\|}\right\|\leq\frac{2\beta}{\alpha}.

(147)

Lemma 17.

(Berry et al. (2017)). Let $|\psi\rangle=\alpha|0\rangle|\psi_{0}\rangle+\sqrt{1-\alpha^{2}}|1\rangle|\psi_{1}\rangle$ and $|\varphi\rangle=\beta|0\rangle|\varphi_{0}\rangle+\sqrt{1-\beta^{2}}|1\rangle|\varphi_{1}\rangle$ , where $|\psi_{0}\rangle$ , $\psi_{1}\rangle$ , $|\varphi_{0}\rangle$ , $|\varphi_{1}\rangle$ are unit vectors, and $\alpha,\beta\in[0,1]$ . Suppose $\||\psi\rangle-|\varphi\rangle\|\leq\delta<\alpha$ . Then $\||\psi_{0}\rangle-|\varphi_{0}\rangle\|\leq\frac{2\delta}{\alpha-\delta}$ .

Lemma 18.

References

Harrow et al. (2009) A. W. Harrow, A. Hassidim, and S. Lloyd, Physical review letters 103, 150502 (2009).
Childs et al. (2017) A. M. Childs, R. Kothari, and R. D. Somma, SIAM Journal on Computing 46, 1920 (2017).
Subaşı et al. (2019) Y. Subaşı, R. D. Somma, and D. Orsucci, Physical review letters 122, 060504 (2019).
Clader et al. (2013) B. D. Clader, B. C. Jacobs, and C. R. Sprouse, Physical review letters 110, 250504 (2013).
Cao et al. (2013) Y. Cao, A. Papageorgiou, I. Petras, J. Traub, and S. Kais, New Journal of Physics 15, 013021 (2013).
Fillion-Gourdeau et al. (2017) F. Fillion-Gourdeau, S. MacLean, and R. Laflamme, Physical Review A 95, 042343 (2017).
Linden et al. (2020) N. Linden, A. Montanaro, and C. Shao, arXiv preprint arXiv:2004.06516 (2020).
Berry (2014) D. W. Berry, Journal of Physics A: Mathematical and Theoretical 47, 105301 (2014).
Berry et al. (2017) D. W. Berry, A. M. Childs, A. Ostrander, and G. Wang, Communications in Mathematical Physics 356, 1057 (2017).
Xin et al. (2020) T. Xin, S. Wei, J. Cui, J. Xiao, I. Arrazola, L. Lamata, X. Kong, D. Lu, E. Solano, and G. Long, Physical Review A 101, 032307 (2020).
Childs and Liu (2020) A. M. Childs and J.-P. Liu, Communications in Mathematical Physics 375, 1427 (2020).
Arrazola et al. (2019) J. M. Arrazola, T. Kalajdzievski, C. Weedbrook, and S. Lloyd, Physical Review A 100, 032306 (2019).
Childs et al. (2020) A. M. Childs, J.-P. Liu, and A. Ostrander, arXiv preprint arXiv:2002.07868 (2020).
Montanaro and Pallister (2016) A. Montanaro and S. Pallister, Physical Review A 93, 032324 (2016).
Costa et al. (2019) P. C. Costa, S. Jordan, and A. Ostrander, Physical Review A 99, 012323 (2019).
Engel et al. (2019) A. Engel, G. Smith, and S. E. Parker, Physical Review A 100, 062315 (2019).
Leyton and Osborne (2008) S. K. Leyton and T. J. Osborne, arXiv preprint arXiv:0812.4423 (2008).
Qian et al. (2019) P. Qian, W.-C. Huang, and G.-L. Long, arXiv preprint arXiv:1903.05608 (2019).
Lubasch et al. (2020) M. Lubasch, J. Joo, P. Moinier, M. Kiffner, and D. Jaksch, Physical Review A 101, 010301 (2020).
Liu et al. (2021) J.-P. Liu, H. Ø. Kolden, H. K. Krovi, N. F. Loureiro, K. Trivisa, and A. M. Childs, Proceedings of the National Academy of Sciences 118 (2021).
Lloyd et al. (2020) S. Lloyd, G. De Palma, C. Gokler, B. Kiani, Z.-W. Liu, M. Marvian, F. Tennie, and T. Palmer, arXiv preprint arXiv:2011.06571 (2020).
Budinski (2021) L. Budinski, arXiv preprint arXiv:2103.03804 (2021).
Chen et al. (2021) Z.-Y. Chen, C. Xue, S.-M. Chen, B.-H. Lu, Y.-C. Wu, J.-C. Ding, S.-H. Huang, and G.-P. Guo, arXiv preprint arXiv:2102.03557 (2021).
Xue et al. (2021) C. Xue, Y. Wu, and G. Guo, in SPIN (World Scientific, 2021), p. 2140004.
Kyriienko et al. (2021) O. Kyriienko, A. E. Paine, and V. E. Elfving, Physical Review A 103, 052416 (2021).
Carleman et al. (1932) T. Carleman et al., Acta Mathematica 59, 63 (1932).
Kowalski and Steeb (1991) K. Kowalski and W.-H. Steeb, Nonlinear dynamical systems and Carleman linearization (World Scientific, 1991).
He (1999) J.-H. He, Computer methods in applied mechanics and engineering 178, 257 (1999).
Babolian et al. (2009) E. Babolian, A. Azizi, and J. Saeidian, Mathematical and computer modelling 50, 213 (2009).
Chakraverty et al. (2019) S. Chakraverty, N. Mahato, P. Karunakar, and T. D. Rao, Advanced numerical and semi-analytical methods for differential equations (Wiley Online Library, 2019).
He (2003) J.-H. He, Applied Mathematics and computation 135, 73 (2003).
He (2005) J.-H. He, Chaos, Solitons & Fractals 26, 695 (2005).
Kerner (1981) E. H. Kerner, Journal of Mathematical Physics 22, 1366 (1981).
Forets and Pouly (2017) M. Forets and A. Pouly, Explicit error bounds for carleman linearization (2017), eprint 1711.02552.
Nielsen and Chuang (2002) M. A. Nielsen and I. Chuang, Quantum computation and quantum information (2002).
Van Loan (1977) C. Van Loan, SIAM Journal on Numerical Analysis 14, 971 (1977).
Koshy (2008) T. Koshy, Catalan numbers with applications (Oxford University Press, 2008).
Brassard et al. (2002) G. Brassard, P. Hoyer, M. Mosca, and A. Tapp, Contemporary Mathematics 305, 53 (2002).
Shijun (1998) L. Shijun, Applied Mathematics and Mechanics 19, 957 (1998).
Ilhan et al. (2021) E. Ilhan, P. Veeresha, and H. M. Baskonus, Chaos, Solitons & Fractals 152, 111347 (2021).
Veeresha et al. (2021a) P. Veeresha, H. M. Baskonus, and W. Gao, Axioms 10, 123 (2021a).
Veeresha et al. (2021b) P. Veeresha, E. Ilhan, D. Prakasha, H. M. Baskonus, and W. Gao, Computer Modeling in Engineering & Sciences 127, 1013 (2021b).
Joseph (2020) I. Joseph, Physical Review Research 2, 043102 (2020).
Engel et al. (2021) A. Engel, G. Smith, and S. E. Parker, Physics of Plasmas 28, 062305 (2021).
Kowalski (1994) K. Kowalski, Methods of Hilbert spaces in the theory of nonlinear dynamical systems (World Scientific, 1994).
Koopman (1931) B. O. Koopman, Proceedings of the national academy of sciences of the united states of america 17, 315 (1931).

$\displaystyle\left\|x_{0,0}\right\rangle$	$\displaystyle=\left\|y_{in}\right\rangle,$	(25)
$\displaystyle\left\|x_{0,j}\right\rangle$	$\displaystyle=\left((Ah)^{j}/j!\right)\left\|x_{0,0}\right\rangle,\quad 1\leq j\leq k,$
$\displaystyle\left\|x_{1,0}\right\rangle$	$\displaystyle=T_{k}(Ah)\left\|x_{0,0}\right\rangle\approx\exp(Ah)\left\|y_{in}\right\rangle,$
$\displaystyle\left\|x_{1,j}\right\rangle$	$\displaystyle=\left((Ah)^{j}/j!\right)\left\|x_{1,0}\right\rangle,\quad 1\leq j\leq k,$
$\displaystyle\left\|x_{2,0}\right\rangle$	$\displaystyle=T_{k}(Ah)\left\|x_{1,0}\right\rangle\approx\exp(2Ah)\left\|y_{in}\right\rangle,$
	$\displaystyle\vdots$
$\displaystyle\left\|x_{m-1,0}\right\rangle$	$\displaystyle=T_{k}(Ah)\left\|x_{m-2,0}\right\rangle\approx\exp(Ah(m-1))\left\|y_{in}\right\rangle,$
$\displaystyle\left\|x_{m-1,j}\right\rangle$	$\displaystyle=\left((Ah)^{j}/j!\right)\left\|x_{m-1,0}\right\rangle,\quad 1\leq j\leq k,$
$\displaystyle\left\|x_{m,0}\right\rangle$	$\displaystyle=T_{k}(Ah)\left\|x_{m-1,0}\right\rangle\approx\exp(Ahm)\left\|y_{in}\right\rangle,$
$\displaystyle\left\|x_{m,j}\right\rangle$	$\displaystyle=\left\|x_{m,0}\right\rangle\approx\exp(Ahm)\left\|y_{in}\right\rangle,\quad 1\leq j\leq p,$

$\displaystyle\\|A\\|\leq$	$\displaystyle\\|{\rm diag}(A_{0,0},A_{1,1},\dots,A_{c,c})\\|+\\|{\rm diag}(A_{0,1},A_{1,2},\dots,A_{c-1,c})\\|$
$\displaystyle=$	$\displaystyle\max_{i=0}^{c}\{\\|A_{i,i}\\|\}+\max_{i=0}^{c-1}\{\\|A_{i,i+1}\\|\}$
$\displaystyle\leq$	$\displaystyle(c+1)(\\|F_{1}\\|+\\|F_{2}\\|).$	(32)

$\displaystyle\|x_{i,j}^{l}\rangle$	$\displaystyle=0,$	$\displaystyle 0\leq i<a,0\leq j\leq k,$	(56)
$\displaystyle\|x_{a,j}^{l}\rangle$	$\displaystyle=0,$	$\displaystyle 0\leq j<b,$
$\displaystyle\|x_{a,j}^{l}\rangle$	$\displaystyle=b!(Ah)^{j-b}/j!\|b^{l}\rangle,$	$\displaystyle b\leq j\leq k$
$\displaystyle\|x_{a+1,0}^{l}\rangle$	$\displaystyle=T_{b,k}(Ah)\|b^{l}\rangle,$
$\displaystyle\|x_{a+1,j}^{l}\rangle$	$\displaystyle=(Ah)^{j}/j!\|x_{a+1,0}^{l}\rangle,$
$\displaystyle\|x_{a+2,0}^{l}\rangle$	$\displaystyle=T_{k}(Ah)\|x_{a+1,0}^{l}\rangle=T_{k}(Ah)T_{b,k}(Ah)\|b^{l}\rangle,$
	$\displaystyle\vdots$
$\displaystyle\|x_{m,0}^{l}\rangle$	$\displaystyle=T_{k}(Ah)\|x_{m-1,0}^{l}\rangle=\left(T_{k}(Ah)\right)^{m-a-1}T_{b,k}(Ah)\|b^{l}\rangle,$
$\displaystyle\|x_{m,j}^{l}\rangle$	$\displaystyle=\|x_{m,0}\rangle=\left(T_{k}(Ah)\right)^{m-a-1}T_{b,k}(Ah)\|b^{l}\rangle,$	$\displaystyle 1\leq j\leq p,$

$\displaystyle\\|(Ah)^{j-b}\\|$	$\displaystyle\leq\\|(Ah)\\|^{j-b}\leq 1,\quad\quad\quad\quad\quad\quad\quad b\leq j\leq k,$	(59)
$\displaystyle\\|T_{b,k}(Ah)\\|$	$\displaystyle\leq T_{b,k}(\\|Ah\\|)\leq e^{\\|Ah\\|}\leq e,$
$\displaystyle\\|T_{k}(Ah)\\|$	$\displaystyle\leq T_{k}(\\|Ah\\|)\leq e^{\\|Ah\\|}\leq e,$
$\displaystyle\\|\left(T_{k}(Ah)\right)^{m-a-1}T_{b,k}(Ah)\\|$	$\displaystyle\leq\\|\left(T_{k}(Ah)\right)^{m-a-1}\\|\\|T_{b,k}(Ah)\\|\leq e(c+2).$

$\displaystyle\\|\|x_{a,j}^{l}\rangle\\|$	$\displaystyle\leq\frac{b!}{j!}\\|(Ah)^{j-b}\\|\\|\|b^{l}\rangle\\|\leq\\|\|b^{l}\rangle\\|,$	$\displaystyle b\leq j\leq k$	(60)
$\displaystyle\\|\|x_{a+1,0}^{l}\rangle\\|$	$\displaystyle\leq\\|T_{b,k}(Ah)\\|\\|\|b^{l}\rangle\\|\leq e\\|\|b^{l}\rangle\\|,$
$\displaystyle\\|\|x_{i,j}^{l}\rangle\\|$	$\displaystyle\leq e(c+2)\\|\|b^{l}\rangle\\|,$	$\displaystyle a+1\leq i\leq m-1,0\leq j\leq k$
$\displaystyle\\|\|x_{m,j}^{l}\rangle\\|$	$\displaystyle\leq e(c+2)\\|\|b^{l}\rangle\\|,$	$\displaystyle 0\leq j\leq p$

Quantum homotopy perturbation method for nonlinear dissipative ordinary differential equations

Abstract

I Introduction

II Quadratic ODEs

III Quantum Homotopy Perturbation Method

III.1 Homotopy perturbation method

III.2 Linear embedding

III.3 Quantum linear ODEs algorithm

III.4 Measurement

IV State Preparation and Oracle Construction

IV.1 State preparation

Lemma 1.

Proof.

IV.2 Oracle construction of AA

Lemma 2.

Proof.

Lemma 3.

Proof.

Lemma 4.

Proof.

Lemma 5.

Proof.

Lemma 6.

Proof.

V Condition Number

Lemma 7.

Proof.

Lemma 8.

Proof.

VI Solution Error

VI.1 Homotopy perturbation method truncation error

Lemma 9.

Proof.

VI.2 Linear ODEs solution error

Lemma 10.

Proof.

VII Success Probability

Lemma 11.

Proof.

Lemma 12.

Proof.

VIII Main Result

Theorem 1.

Proof.

IX Conclusion and Discussion

Acknowledgements

Appendix

Lemma 13.

Proof.

Lemma 14.

Proof.

Lemma 15.

Proof.

Lemma 16.

Lemma 17.

Lemma 18.

References

IV.2 Oracle construction of $A$