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

"$VERSION"

## Non-Linear Differential Equation Solver

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

The algorithm makes use of two sub-routines :

- A function transformation routine - this is a mapping from
to a polynomial,*z*. The function transformation is effected using a Hamiltonian simulation. *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

To begin, we encode the vector ** z , **after normalizing it,

**in a state**

The tensor product,

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

For efficient simulation, the mapping has to be *sparse* in nature. In general, the functions * *

To actually carry out the simulation, we need to build a hermitian operator containing*von Neumann measurement* prescription).

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

We can construct

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) isapprox.(r_out, out[2:3]*sqrt(2), atol = 1e-3) |> all

For forward Euler, we perform multiple iterations of the function transformation sub-routine. If the functions

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) isapprox.(r_out,real(out), atol = 1e-3) |> all

## Simulating the Wave Equation (in 1D)

The wave equation is given by

in one dimension. As is practice, the spacial part is discretised to leave a differential equation in time alone. I'll set

where

We can

This implies,

If

where

The value of

vertx = 7 ege = 8 #Constructing incidence matrix B B = zeros(vertx,ege) 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] = ϕ[1:vertx] for i in 2:n r, N = taylorsolve(H,ϕ,k,t) ϕ = N*vec(state(r)) res[i] = ϕ[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 = (u[:,1]) u2 = (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] = (sol1.u[i][:,1]) end isapprox.(res1,s1, atol = 1e-2) |> all

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