Added Higher-order linear multistep methods using HHL (#4)

Added Higher-order linear multistep

Author: dgan181

Project.toml

@@ -4,6 +4,7 @@ authors = ["Roger-luo <hiroger@qq.com>"]
 version = "0.1.0"
 
 [deps]
+DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LuxurySparse = "d05aeea4-b7d4-55ac-b691-9e7fabb07ba2"
@@ -12,8 +13,8 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Yao = "5872b779-8223-5990-8dd0-5abbb0748c8c"
 
 [extras]
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
 test = ["Test", "OrdinaryDiffEq"]

src/lin_diffEq_HHL.jl

@@ -1,10 +1,14 @@
-export Array_QuEuler, prepare_init_state, solve_QuEuler
+export array_qudiff, prepare_init_state, LDEMSAlgHHL, bval, aval
+export QuEuler, QuLeapfrog, QuAB2, QuAB3, QuAB4
+export QuLDEMSProblem
 
+using DiffEqBase
 """
     Based on : arxiv.org/abs/1010.2745v2
 
-    QuArray_Euler(N_t,N,h,A)
-    prepare_init_state(b,x,h,N_t)
+    * array_qudiff(N_t,N,h,A) - generates matrix for k-step solver
+    * prepare_init_state(b,x,h,N_t) - generates inital states
+    * solve_qudiff - solver
 
     x' = Ax + b
 
@@ -15,72 +19,139 @@ export Array_QuEuler, prepare_init_state, solve_QuEuler
     * h - step size.
     * tspan - time span.
 """
-function prepare_init_state(tspan::Tuple,x::Vector,h::Float64,g::Function)
-  init_state = x;
-  N_t = round(2*(tspan[2] - tspan[1])/h + 3)
-  b = similar(g(1))
-  for i = 1:N_t
-   if i < (N_t+1)/2
-      b = g(h*i + tspan[1])
-      init_state = [init_state;h*b]
-    else
-      init_state = [init_state;zero(b)]
-   end
-
-  end
-  init_state = [init_state;zero(init_state)]
-  init_state
+
+"""
+    LDEMSAlgHHL
+    * step - step for multistep method
+    * α - coefficients for xₙ
+    * β - coefficent for xₙ'
+"""
+abstract type QuODEAlgorithm <: DiffEqBase.AbstractODEAlgorithm end
+#abstract type QuODEProblem{uType,tType,isinplace} <: DiffEqBase.AbstractODEProblem{uType,tType,isinplace} end
+abstract type LDEMSAlgHHL <: QuODEAlgorithm end
+
+struct QuLDEMSProblem{F,C,U,T} #<: QuODEProblem{uType,tType,isinplace}
+    A::F
+    b::C
+    u0::U
+    tspan::NTuple{2,T}
+
+    #function QuLDEMSProblem(A,b,u0,tspan)
+    #  new{typeof(u0),typeof(tspan),false,typeof(A),typeof(b)}(A,b,u0,tspan)
+    #end
 end
 
-function Array_QuEuler(tspan::Tuple,h::Float64,g::Function)
-  N_t = round(2*(tspan[2] - tspan[1])/h + 3)
-  I_mat = Matrix{Float64}(I, size(g(1)));
-  A_ = I_mat;
-  zero_mat = zero(I_mat);
-  tmp_A = Array{Float64,2}(undef,size(g(1)))
-
-  for j = 2: N_t +1
-    A_ = [A_ zero_mat]
-  end
-
-  for i = 2 : N_t+1
-    tmp_A =  g(i*h + tspan[1])
-    tmp_A = -1*(I_mat + h*tmp_A)
-    tA_ = copy(zero_mat)
-    if i<3
-      tA_ = [tmp_A I_mat]
-    else
-      for j = 1:i-3
-        tA_ = [tA_ zero_mat]
-      end
-      if i < (N_t + 1)/2 + 1
-        tA_ = [tA_ tmp_A I_mat]
-      else
-        tA_ = [tA_ -1*I_mat I_mat]
-      end
+"""
+    Explicit Linear Multistep Methods
+"""
+struct QuEuler{T}<:LDEMSAlgHHL
+    step::Int
+    α::Vector{T}
+    β::Vector{T}
+
+    QuEuler(::Type{T} = Float64) where {T} = new{T}(1,[1.0,],[1.0,])
+end
+
+struct QuLeapfrog{T}<:LDEMSAlgHHL
+    step::Int
+    α::Vector{T}
+    β::Vector{T}
+
+    QuLeapfrog(::Type{T} = Float64) where {T} = new{T}(2,[0, 1.0],[2.0, 0])
+end
+struct QuAB2{T}<:LDEMSAlgHHL
+    step::Int
+    α::Vector{T}
+    β::Vector{T}
+
+    QuAB2(::Type{T} = Float64) where {T} = new{T}(2,[1.0, 0], [1.5, -0.5])
+end
+struct QuAB3{T}<:LDEMSAlgHHL
+    step::Int
+    α::Vector{T}
+    β::Vector{T}
+    QuAB3(::Type{T} = Float64) where {T} = new{T}(3,[1.0, 0, 0], [23/12, -16/12, 5/12])
+end
+struct QuAB4{T}<:LDEMSAlgHHL
+    step::Int
+    α::Vector{T}
+    β::Vector{T}
+
+    QuAB4(::Type{T} = Float64) where {T} = new{T}(4,[1.0, 0, 0, 0], [55/24, -59/24, 37/24, -9/24])
+end
+
+function bval(alg::LDEMSAlgHHL,t,h,g::Function)
+    b = zero(g(1))
+    for i in 1:(alg.step)
+        b += alg.β[i]*g(t-(i-1)*h)
     end
-    if i< N_t+1
-      for j = i+1: N_t+1
-        tA_ = [tA_ zero_mat]
-      end
+    return b
+end
+
+function aval(alg::LDEMSAlgHHL,t,h,g::Function)
+    sz, = size(g(1))
+    A = Array{ComplexF64}(undef,sz,(alg.step + 1)*sz)
+    i_mat = Matrix{Float64}(I, size(g(1)))
+    A[1:sz,sz*(alg.step) + 1:sz*(alg.step + 1)] = i_mat
+    for i in 1:alg.step
+        A[1:sz,sz*(i - 1) + 1: sz*i] = -1*(alg.α[alg.step - i + 1]*i_mat + h*alg.β[alg.step - i + 1]*g(t - (alg.step - i)*h))
     end
-    A_ = [A_;tA_]
+    return A
+end
 
-  end
-    A_ = [zero(A_) A_;A_' zero(A_)]
-    A_
+function prepare_init_state(tspan::NTuple{2, Float64},x::Vector,h::Float64,g::Function,alg::LDEMSAlgHHL)
+    N_t = round(Int, (tspan[2] - tspan[1])/h + 1) #number of time steps
+    N = nextpow(2,2*N_t + 1) # To ensure we have a power of 2 dimension for matrix
+    sz, = size(g(1))
+    init_state = zeros(ComplexF64,2*(N)*sz)
+    #inital value
+    init_state[1:sz] = x
+    for i in 2:N_t
+        b = bval(alg,h*(i - 1) + tspan[1],h,g)
+        init_state[Int(sz*(i - 1) + 1):Int(sz*(i))] = h*b
+    end
+    return init_state
 end
 
-function solve_QuEuler(A::Function, b::Function, x::Vector,tspan::Tuple, h::Float64)
-
-  mat = Array_QuEuler(tspan,h,A)
-  state = prepare_init_state(tspan,x,h,b)
-  λ = maximum(eigvals(mat))
-  C_value = minimum(eigvals(mat) .|> abs)*0.1;
-  n_reg = 12;
-  mat = 1/(λ*2)*mat
-  state = state*1/(2*λ) |> normalize!
-  res = hhlsolve(mat,state, n_reg, C_value)
-  res = res/λ
-  res
+function array_qudiff(tspan::NTuple{2, Float64},h::Float64,g::Function,alg::LDEMSAlgHHL)
+    sz, = size(g(1))
+    i_mat = Matrix{Float64}(I, size(g(1)))
+    N_t = round(Int, (tspan[2] - tspan[1])/h + 1) #number of time steps
+    N = nextpow(2,2*N_t + 1) # To ensure we have a power of 2 dimension for matrix
+    A_ = zeros(ComplexF64, N*sz, N*sz)
+    # Generates First two rows
+    @inbounds A_[1:sz, 1:sz] = i_mat
+    @inbounds A_[sz + 1:2*sz, 1:sz] = -1*(i_mat + h*g(tspan[1]))
+    @inbounds A_[sz + 1:2*sz,sz+1:sz*2] = i_mat
+    #Generates additional rows based on k - step
+    for i in 3:alg.step
+        @inbounds A_[sz*(i - 1) + 1:sz*i, sz*(i - 3) + 1:sz*i] = aval(QuAB2(),(i-2)*h + tspan[1],h,g)
+    end
+    for i in alg.step + 1:N_t
+        @inbounds A_[sz*(i - 1) + 1:sz*(i), sz*(i - alg.step - 1) + 1:sz*i] = aval(alg,(i - 2)*h + tspan[1],h,g)
+    end
+    #Generates half mirroring matrix
+    for i in N_t + 1:N
+        @inbounds A_[sz*(i - 1) + 1:sz*(i), sz*(i - 2) + 1:sz*(i - 1)] = -1*i_mat
+        @inbounds A_[sz*(i - 1) + 1:sz*(i), sz*(i - 1) + 1:sz*i] = i_mat
+    end
+    A_ = [zero(A_) A_;A_' zero(A_)]
+    return A_
 end
+
+function DiffEqBase.solve(prob::QuLDEMSProblem{F,C,U,T}, alg::LDEMSAlgHHL, dt = (prob.tspan[2]-prob.tspan[1])/100, n_reg::Int = 12) where {F,C,U,T}
+    A = prob.A
+    b = prob.b
+    tspan = prob.tspan
+    x = prob.u0
+
+    mat = array_qudiff(tspan, dt, A, alg)
+    state = prepare_init_state(tspan, x, dt, b, alg)
+    λ = maximum(eigvals(mat))
+    C_value = minimum(eigvals(mat) .|> abs)*0.01;
+    mat = 1/(λ*2)*mat
+    state = state*1/(2*λ) |> normalize!
+    res = hhlsolve(mat,state, n_reg, C_value)
+    res = res/λ
+    return res
+end;

test/lin_diffEq_test.jl

@@ -4,29 +4,47 @@ using QuAlgorithmZoo
 using Test, LinearAlgebra
 using OrdinaryDiffEq
 
-function diffEq_problem(nbit::Int)
+function diffeq_problem(nbit::Int)
     siz = 1<<nbit
     A = (rand(ComplexF64, siz,siz))
-    A = (A+A')/2
+    A = (A + A')/2
     b = normalize!(rand(ComplexF64, siz))
     x = normalize!(rand(ComplexF64, siz))
     A, b, x
 end
 
-@testset "Linear_differential_equation_Euler_HHL" begin
-N = 1
-h = 0.02
-tspan = (0.0,0.6)
-M, v, x = diffEq_problem(N)
-A(t) = M
-b(t) = v
-res = solve_QuEuler(A, b, x, tspan, h)
+@testset "Linear_differential_equations_HHL" begin
+    N = 1
+    h = 0.1
+    tspan = (0.0,0.6)
+    N_t = round(Int, 2*(tspan[2] - tspan[1])/h + 3)
+    M, v, x = diffeq_problem(N)
+    A(t) = M
+    b(t) = v
+    n_reg = 12
+    f(u,p,t) = M*u + v;
+    prob = ODEProblem(f, x, tspan)
+    qprob = QuLDEMSProblem(A, b, x, tspan)
+    sol = solve(prob, Tsit5(), dt = h, adaptive = :false)
+    s = vcat(sol.u...)
 
-f(u,p,t) = M*u + v;
-prob = ODEProblem(f,x,tspan)
-sol = solve(prob,Euler(),dt = 0.02,adaptive = :false)
-s = vcat(sol.u...)
-N_t = round(2*(tspan[2] - tspan[1])/h + 3);
-r = res[Int((N_t+1)*2+2^N+1): Int((N_t+1)*2+2^N + N_t - 1)] # range of relevant values in the obtained state.
-@test isapprox.(s,r,atol = 0.05) |> all
-end
+    res = solve(qprob, QuEuler(), h, n_reg)
+    r = res[(N_t + 1)*2 + 2^N - 1: (N_t + 1)*2 + 2^N + N_t - 3] # range of relevant values in the obtained state.
+    @test isapprox.(s, r, atol = 0.5) |> all
+
+    res = solve(qprob, QuLeapfrog(), h, n_reg)
+    r = res[(N_t + 1)*2 + 2^N - 1: (N_t + 1)*2 + 2^N + N_t - 3] # range of relevant values in the obtained state.
+    @test isapprox.(s, r, atol = 0.3) |> all
+
+    res = solve(qprob, QuAB2(), h,n_reg)
+    r = res[(N_t + 1)*2 + 2^N - 1: (N_t + 1)*2 + 2^N + N_t - 3] # range of relevant values in the obtained state.
+    @test isapprox.(s, r, atol = 0.3) |> all
+
+    res = solve(qprob, QuAB3(), h,n_reg)
+    r = res[(N_t + 1)*2 + 2^N - 1: (N_t + 1)*2 + 2^N + N_t - 3] # range of relevant values in the obtained state.
+    @test isapprox.(s, r, atol = 0.3) |> all
+
+    res = solve(qprob, QuAB4(), h ,n_reg)
+    r = res[(N_t + 1)*2 + 2^N - 1: (N_t + 1)*2 + 2^N + N_t - 3] # range of relevant values in the obtained state.
+    @test isapprox.(s, r, atol = 0.3) |> all
+end;