朱莉娅：Runge Kutta方法的计算时间

问题描述

我正在尝试计算矩阵形式的ODE， 即（∂/∂t）ψ=-iHψ，其中ψ是一个向量，H是一个（时间独立的）矩阵。
我有两个问题。

1。 DifferentialEquations.jl和Runge Kutta代码之间的计算时间差异
我用两种方法计算了上式。
DifferentialEquations.jl

    using Linearalgebra    
    using OrdinaryDiffEq   
    using DifferentialEquations

    #matrix
    function Ham(Lx,Ly,Lz) 
     4*SINk₁(Lx,Lz) + 3*im*SINk₂(Lx,Lz) + COSk₃(Lx,Lz) -im*EYE(Lx,Lz)
    end

    # Initial conditions
    ψ0 = Complex{Float64}[] 
    σ = sqrt(5/2)

    for iz = 1:Lz　　　　　　　　
     for ix = 1:Lx
      for iy=1:Ly             
        gauss = (1/(((sqrt(2*π))*σ)^3))exp(-((ix-5)^2 + (iy-5)^2 + (iz-5)^2)/(2*(σ)^2))
        push!(ψ0,gauss) 
      end
     end
    end

    #normalization
    ψ0 = ψ0./norm(ψ0)

    #time span
    tspan = (0.0,20.0) 

    #ODE
    function time_evolution(ψdot::Array{Complex{Float64},1},ψ::Array{Complex{Float64},para::Float64,t::Float64)
     ψdot.=-im.*Ham(Lx,Lz)*ψ
    end

    #solve the equation
    prob = ODEProblem(time_evolution,ψ0,tspan)
    @time sol = solve(prob)

其中SINk₁（Lx，Ly，Lz），SINk²（Lx，Ly，Lz），COSk₃（Lx，Ly，Lz），EYE（Lx，Ly，Lz）是

 function EYE(Lx,Lz)
     N = Lx*Ly*Lz
     mat = Matrix{Complex{Float64}}(I,N,N) 
     return mat
 end


 function SINk₁(Lx,Lz)
  N = Lx*Ly*Lz

mat = zeros(Complex{Float64},N) 

for ix = 1:Lx
    for iy = 1:Ly 
        for iz = 1:Lz 
            for dx in -1:1             
             jx = ix + dx 
                for dy in -1:1
                jy = iy + dy
                    for dz in -1:1
                    jz = iz + dz
                    
                    
                        ii = (iz-1)*Lx*Ly + (ix-1)*Ly + iy 
                        jj = (jz-1)*Lx*Ly + (jx-1)*Ly + jy

                            if 1 <= jx <= Lx && 1 <= jy <= Ly && 1 <= jz <= Lz
                                                     
                                if dx == +1 && dy == 0 && dz == 0
                                    mat[ii,jj] += -(im/2)　
                                end
                
                                if dx == -1 && dy == 0 && dz == 0
                                    mat[ii,jj] += im/2
                                end
                                   
                            end    
                
                    end
                end
            end
        end
    end
end

return mat
end


function SINk₂(Lx,N) 

for ix = 1:Lx
    for iy = 1:Ly 
        for iz = 1:Lz 
            for dx in -1:1             
             jx = ix + dx 
                for dy in -1:1
                jy = iy + dy
                    for dz in -1:1
                    jz = iz + dz
                    
                    
                        ii = (iz-1)*Lx*Ly + (ix-1)*Ly + iy 
                        jj = (jz-1)*Lx*Ly + (jx-1)*Ly + jy

                            if 1 <= jx <= Lx && 1 <= jy <= Ly && 1 <= jz <= Lz
                                                     
                                if dx == 0 && dy == +1 && dz == 0
                                    mat[ii,jj] += -im/2
                                end
                
                                if dx == 0 && dy == -1 && dz == 0
                                    mat[ii,jj] += im/2
                                end
                                   
                            end    
                
                    end
                end
            end
        end
    end
end

return mat
end


function COSk₃(Lx,Lz)
 N = Lx*Ly*Lz

mat = zeros(Complex{Float64},N) 

for ix = 1:Lx
    for iy = 1:Ly 
        for iz = 1:Lz 
            for dx in -1:1             
             jx = ix + dx 
                for dy in -1:1
                jy = iy + dy
                    for dz in -1:1
                    jz = iz + dz
                    
                    
                        ii = (iz-1)*Lx*Ly + (ix-1)*Ly + iy 
                        jj = (jz-1)*Lx*Ly + (jx-1)*Ly + jy

                            if 1 <= jx <= Lx && 1 <= jy <= Ly && 1 <= jz <= Lz
                                                     
                                if dx == 0 && dy == 0 && dz == +1
                                    mat[ii,jj] += 1/2　
                                end
                
                                if dx == 0 && dy == 0 && dz == -1
                                    mat[ii,jj] += 1/2
                                end
                                   
                            end    
                
                    end
                end
            end
        end
    end
end

return mat
end

对于矩阵大小Lx = Ly = Lz = 10，这大约需要50秒：

  49.591043 seconds (28.10 k allocations: 141.020 GiB,21.55% gc time)

Runge Kutta代码

function rk4(ψ::Array{Complex{Float64},H::Array{Complex{Float64},2},δt::Float64)

f(a) = -im*H*a

k1 = f(ψ)
k2 = f(ψ + (1/2)*δt*k1)
k3 = f(ψ + (1/2)*δt*k2)
k4 = f(ψ + δt*k3)
ψ = ψ + δt*(1/6)*(k1 + 2*k2 + 2*k3 + k4)

return ψ

end


function RK4(ti::Int64,tf::Int64,interval::Int64,ψ0::Array{Complex{Float64},1})

#time span
δt = (tf-ti)/interval

size = length(ψ0)

#container of solutions
list_t = zeros(Float64,interval+1)
list_ψ = zeros(Complex{Float64},interval+1,size)

#initial condition
list_t[1] = ti
list_ψ[1,1:size] = ψ0

 #store solutions
for i in 2:interval+1
    list_ψ[i,1:size] = rk4(list_ψ[i-1,1:size],H,δt)
    list_t[i] = list_t[i-1]+δt
end

return list_t,list_ψ

end


#matrix
H = Ham(Lx,Lz) 

#time span
ti = 0
tf = 20
interval = 60

@time list_t1,list_ψ1 = RK4(ti,tf,interval,ψ0)

对于矩阵大小Lx = Ly = Lz = 10和相同的δt，这大约需要2.0秒：

 1.882741 seconds (1.45 k allocations: 3.592 GiB,24.03% gc time)

是什么原因导致差异？

2。 Julia和Python之间的计算时间差异

即使在Runge Kutta代码中，大型矩阵也要花费太多时间。
例如，Lx = Ly = Lz = 20，即8000×8000矩阵，大约需要130秒。

但是，在python中，我听说Runge Kutta 4阶方法即使对于100000 * 100000矩阵也只需要几秒钟。
为什么这么快？
茱莉亚有相同的计算时间吗？

（2020/10/24）
我在茱莉亚使用adaptive=false检查了计算时间。 由于安装BenchmarkTools由于某种原因而失败，因此我使用了@time宏。
对于上述Ham（Lx，Ly，Lz），ψ0和Lx = Ly = Lz = 10，

RK4

#time span
ti = 0
tf = 20
interval = 60  #dt=1/3

@time list_t1,ψ0)

0.983133 seconds (1.40 M allocations: 3.657 GiB,18.14% gc time)

DifferentialEquations.jl

控制准确性

#time span
tspan = (0.0,20.0) 

prob = ODEProblem(time_evolution,tspan)
@time sol = solve(prob)

32.249755 seconds (28.10 k allocations: 141.020 GiB,29.34% gc time)

无法控制准确性

#time span
tspan = (0.0,20.0)

prob = ODEProblem(time_evolution,tspan)
@time sol = solve(prob,RK4(),dt=1/3,adaptive=false)

7.910114 seconds (6.38 k allocations: 35.919 GiB,29.14% gc time)

由于固定dt，在DifferentialEquations.jl中它肯定变得更快。
但是，即使使用相同的dt，它也不适合RK4。
我不知道为什么...

无论我怎么努力，它都不能像python一样快？

（2020/10/29）
根据我朋友的建议，改进了Runge Kutta代码。

function rk4(ψ::Array{Complex{Float64},δt::Float64)
    
    f(a) = -im*mul!(similar(a),a)
    
    k1 = f(ψ)
    k2 = f(ψ + (1/2)*δt*k1)
    k3 = f(ψ + (1/2)*δt*k2)
    k4 = f(ψ + δt*k3)
    ψ = ψ + δt*(1/6)*(k1 + 2*k2 + 2*k3 + k4)
    
    return ψ
    
end

然后，使用H = sparse(Ham(Lx,Lz))对相同的设置产生0.007810 seconds (1.33 k allocations: 19.391 MiB)。

解决方法

您没有比较相同的东西。 DifferentialEquations.jl是一种自适应时间步长，可以达到公差。另一个RK4是固定时间步长，会出现任何错误。当然，计算精度更高几个数量级的解决方案会更昂贵！他们使用的时间步数大不相同，因此计算成本当然是不同的。试试：

using BenchmarkTools
@btime sol = solve(prob,RK4(),dt=...,adaptive=false)

并且如果两者都使用相同的dt，则它们应该相同。

为什么这么快？

因为它不控制精度，所以它确实很快给出了一个不受控制的错误量的答案，您可以强制DifferentialEquations.jl进行此操作，但是由于它试图解决ODE，因此默认情况下不会这样做。