Divyanshu Gupta / Sep 24 2019
Remix of Julia by Nextjournal

JSoC ' 19 : Non-Linear Differential Equation Solver and Simulating of the Wave Equation


Non-Linear Differential Equation Solver

The problem at hand is a set of differential equation of the form:

dzidt=fi(z)i=1,2...n.\frac{d{z_{i}}}{dt} = f_{i}(\textit{\textbf{z}}) \Big|_{i = 1,2...n} .

The algorithm makes use of two sub-routines :

  1. A function transformation routine - this is a mapping from z to a polynomial, . The function transformation is effected using a Hamiltonian simulation. The Taylor Truncation method is used here.
  2. A differential equations solver - Forward Euler is used for this purpose. We make use of the mapping,.

The non-linear differential equation solver included is for solving quadratic-type differential equation alone. Though, in principle it can be extended to higher degrees, in that case simulating the hamiltonian is quite costly.

Function transformation

The functions being quadratic can be expressed as a sum of monomials with coefficients .

fi=k,l=0nαiklzkzl,f_i = \sum_{k,l = 0 \rightarrow n} \alpha_i^{kl} z_k z_l,

is equal to .

To begin, we encode the vector z , after normalizing it, in a state

ϕ=120+12z,|\phi\rangle = \frac{1}{\sqrt{2}} |0\rangle + \frac{1}{\sqrt{2}} |z\rangle,

The tensor product, gives us the a set of all possible monomials in a quadratic equation.

ϕϕ=120+12k,l=0nzkzlkl.|\phi\rangle|\phi\rangle = \frac{1}{2} |0\rangle + \frac{1}{2} \sum_{k,l = 0}^n z_k z_l|k\rangle |l\rangle.

What's required now is an operator that assigns corresponding coefficients to each monomial. We define an operator ,

A=i,k,l=0naikli0kl,A = \sum_{i,k,l = 0}^{n}a_i^{kl}|i0⟩⟨kl|,

, for , and is zero otherwise.

acting on gives us the desired result

Aϕϕ=i,k,l=0naiklzkzli0.A |\phi\rangle|\phi\rangle = \sum_{i,k,l = 0}^{n}a_i^{kl}z_k z_l|i⟩|0⟩.

For efficient simulation, the mapping has to be sparse in nature. In general, the functions will not be measure preserving i.e. they do not preserve the norm of their arguments. In that case, the operator needs to be adjusted by appropriately multiplying its elements byor .

To actually carry out the simulation, we need to build a hermitian operator containing. A well-known trick is to write the hamiltonian (this is von Neumann measurement prescription). is simulated (using Taylor Truncation method). The resulting state is post-selected on to precisely get the what we are looking for.

Usage and Example

using QuDiffEq
using Random, Test, Yao, OrdinaryDiffEq, LinearAlgebra

We want to effect the following function transformation

function f(du,u,p,t)
    du[1] = -3*u[1]^2 + u[2]
    du[2] = -u[2]^2 - u[1]*u[2]
f (generic function with 1 method)

We can construct with the aforementioned prescription

N = 2
k = 3
ϵ = 1e-3
siz = nextpow(2, N + 1)
A = zeros(ComplexF64,2^(siz),2^(siz))
A[1,1] = ComplexF64(1)
A[5,3] = ComplexF64(1)
A[5,6] = ComplexF64(-3)
A[9,11] = ComplexF64(-1)
A[9,7] = ComplexF64(-1);

Comparing the result of the simulation with the actual value,

x = [0.6, 0.8];
tspan = (0.0,0.4)
qprob = QuLDEProblem(A, x, tspan)  
r, N = func_transform(qprob.A, qprob.b, k, ϵ)
out = N*vec(state(r))
r_out = zero(x)
f(r_out, x,1,1)
@test isapprox.(r_out, out[2:3]*sqrt(2), atol = 1e-3) |> all
Test Passed

For forward Euler, we perform multiple iterations of the function transformation sub-routine. If the functions are not measure preserving, we need to alterat every step.

I define the algorithm as QuNLDE(). It takes the arguments k, to select the order in the Taylor series, and ϵ, to set the precision in the Hamiltonian simulation. The algorithm is used as:


Note : In a real quantum computation, for m steps, the every step of the iteration needs to be performed times to obtain a success probability of . But for the sake of simplicity, I will forgo that.

Comparing the result with Euler() ,

prob = ODEProblem(f, x, tspan)
sol = solve(prob, Euler(), dt = 0.1, adaptive = false)
r_out = transpose(hcat(sol.u...))
out = solve(qprob, QuNLDE(3), dt = 0.1)
@test isapprox.(r_out,real(out), atol = 1e-3) |> all
Test Passed

Simulating the Wave Equation (in 1D)

The wave equation is given by

d2ϕdt2=c2d2ϕdx2,\frac{d^2 \phi}{dt^2} = c^2\frac{d^2 \phi}{dx^2},

in one dimension. As is practice, the spacial part is discretised to leave a differential equation in time alone. I'll set . The spacial part can be approximated by the following,

[d2dx2]ϕjϕj+12ϕj+ϕj1a2:=1a2[Lϕ]j,\Big[\frac{d^2}{dx^2}\Big]\phi_j \rightarrow \frac{\phi_{j+1} - 2\phi_j + \phi_{j-1}}{a^2} := \frac{-1}{a^2}[L \phi]_j,

where is a matrix acting on the vector and is small. For a given hermitian , the Schrodinger’s equation can be written as,

ddt(ϕVϕE)=iH(ϕVϕE).()\frac{d}{dt}\begin{pmatrix}\phi_V \\ \phi_E \end{pmatrix} = -iH\begin{pmatrix}\phi_V \\ \phi_E \end{pmatrix}.\:\:\:\:\:\:\:\:\:\:\:\: (*)

We can to be,

H=1a(0BB0).H = \frac{1}{a}\begin{pmatrix} 0 & B \\ B^{\dagger} & 0\end{pmatrix}.

This implies,

d2dt2(ϕVϕE)=H2(ϕVϕE)=1a2(BB00BB)(ϕVϕE).\frac{d^2}{dt^2}\begin{pmatrix}\phi_V \\ \phi_E \end{pmatrix} = -H^2\begin{pmatrix}\phi_V \\ \phi_E \end{pmatrix} = \frac{-1}{a^2}\begin{pmatrix} BB^{\dagger} & 0\\ 0 & B^{\dagger}B \end{pmatrix} \begin{pmatrix}\phi_V \\ \phi_E \end{pmatrix}.

If , the equation in gives us the wave equation. The here is generally referred to as the incidence matrix. A hamiltonian simulation of will, thus, result in the simulation of the wave equation itself. At the very beginning,store the initial positions and velocities of points on the graph (discretised space). We have,

ϕ0=ϕV,dϕ0dt=iaBϕE,\begin{aligned} \phi_0 & = \phi_V ,\\ \frac{d\phi_0}{dt} & = -\frac{i}{a}B \phi_E, \end{aligned}

where are the initial values. In the general case, can be constructed by solving the equation , with , . One may use HHL for this process.

The value of is determined by the graph and the boundary conditions. The following example is for a line segment, that is divided into seven segments, that follows stationary () Dirichlet boundary condition.

vertx = 7
ege = 8

#Constructing incidence matrix B
B = zeros(vertx,ege)
@inbounds for i in 1:vertx
    B[i,i] = -1
    B[i,i+1] = 1
B[1,1] = 1

de = 0.5
sn = sin.((0.0:de:(vertx-1)*de)*2*pi/((vertx-1)*de))
u0 = ComplexF64.(sn)

# Intial conditions (stationary to begin with)
u1 = [u0; zero(u0); 0.0; 0.0]
u_ = Float64[u0 zero(u0)]

k = 2   # order in Taylor expansion
t = 1e-2 # time step
B_t = transpose(B)
n = 11  # number of steps
a = 1e-1 # spactial discretization

function make_hamiltonian(B,B_t,a)
    vertx,ege = size(B)
    n = nextpow(2,vertx+ege)
    H = zeros(n,n)
    H[1:vertx,vertx+1:vertx+ege] = B
    H[vertx+1:vertx+ege,1:vertx] = B_t
    H = -im/a*H
    return H

function do_pde(ϕ,B,B_t,k,t,a,n)
    vertx, = size(B)
    H = make_hamiltonian(B,B_t,a)
    res = Array{Array{ComplexF64,1},1}(undef,n)
    res[1] = @view ϕ[1:vertx]
    for i in 2:n
    r, N = taylorsolve(H,ϕ,k,t)
    ϕ = N*vec(state(r))
    res[i] = @view ϕ[1:vertx]
    res_real = real(res)
    return res_real

res1 = do_pde(u1,B,B_t,k,t,a,n)

const D1 = -1/(a*a)*B*B_t # Laplacian Operator
function f(du,u,p,t)
    buffer, D = p
    u1 = @view(u[:,1])
    u2 = @view(u[:,2])
    mul!(buffer, D, u2)
    Du = buffer

    du[:,1] = Du
    du[:,2] = u1

tspan = (0.0,0.1)
prob1 = ODEProblem(f,u_,tspan,(zero(u_[:,1]), D1))
sol1 = solve(prob1,Tsit5(),dt=0.01,adaptive = false)
s1 = Array{Array{Float64,1},1}(undef,n)
for i in 1:n
    s1[i] = @view (sol1.u[i][:,1])

@test isapprox.(res1,s1, atol = 1e-2) |> all
Test Passed


  1. A quantum algorithm to solve nonlinear differential equations : https://arxiv.org/abs/0812.4423
  2. Quantum Algorithm for Simulating the Wave Equation: https://arxiv.org/abs/1711.05394