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

"$VERSION"

0.7s

Julia

"1.1.0"

Non-Linear Differential Equation Solver

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

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

The algorithm makes use of two sub-routines :

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.
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 .

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

|\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.

|\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 = \sum_{i,k,l = 0}^{n}a_i^{kl}|i0⟩⟨kl|,

, for , and is zero otherwise.

acting on gives us the desired result

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

0.6s

Julia

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]
end

1.2s

Julia

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);

0.4s

Julia

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

20.2s

Julia

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:

solve(::QuODEProblem,::QuNLDE)

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

6.1s

Julia

Test Passed

Simulating the Wave Equation (in 1D)

The wave equation is given by

\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,

\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,

\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 = \frac{1}{a}\begin{pmatrix} 0 & B \\ B^{\dagger} & 0\end{pmatrix}.

This implies,

\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,

\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.

3.1s

Julia

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
end
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
end

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]
    end
    res_real = real(res)
    return res_real
end

#Dirichlet
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
end

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])
end

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

3.1s

Julia

Test Passed

References

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