提高 Julia 中对数尺度的绘图质量

Question

我有一个求解两个 ODE 的函数并加入该解决方案，该解决方案在线性尺度上给出了一个非常好的图，但根据我使用的参数，它在对数尺度上的质量下降很大。在下面的代码中，我绘制了两组参数的解，其中您可以看到第一组不平滑，而第二组还可以。如果我尝试使用第二个 ODE 中的 saveat 选项获得更平滑的可视化，solclass = solve(probclass,Tsit5(),saveat=0.001)，则在绘制第二组时会出现错误：ArgumentError: range step cannot be zero。 除了手动更改 saveat 选项之外，还有其他方法可以获得平滑的线性日志？另外，我尝试过使用其他一些后端，但它们在绘制解决方案时出错。

using DifferentialEquations, Plots, RecursiveArrayTools

function alpha_of_phi!(s2,d,a0,ϕ₀)
    # α in the quantum phase 
    function quantum!(dv,v,p,ϕ)
        s2,d=p
        α = v[1]
        dv[1] = (ϕ*s2*sin(2*d*α)+2*d*sinh(s2*α*ϕ))/(-α*s2*sin(2*d*α)+2*d*cos(2*d*α)+2*d*cosh(s2*α*ϕ))
    end
    # When dα/dϕ = 1, we reach the classical regime, and stop the integration
    condition(v,ϕ,integrator) = (ϕ*s2*sin(2*d*v[1])+2*d*sinh(s2*v[1]*ϕ))/(-v[1]*s2*sin(2*d*v[1])+2*d*cos(2*d*v[1])+2*d*cosh(s2*v[1]*ϕ))==1.0
    affect!(integrator) = terminate!(integrator)
    cb = DiscreteCallback(condition,affect!)
    # Initial Condition at the bounce
    α₀ = [a0]
    classspan = (0,ϕ₀)
    probquant = ODEProblem(quantum!,α₀,classspan,(s2,d))
    solquant = solve(probquant,Tsit5(),callback=cb)

    # α in the classical phase
    function classic!(du,u,p,ϕ)
        αc = u[1]
        dαc = u[2]
        du[1] = dαc
        du[2] = 3*(-dαc^3/sqrt(2)+dαc^2+dαc/sqrt(2)-1)
    end
    # IC retrieved from end of quantum phase
    init = [last(solquant);1.0]
    classspan = (last(solquant.t),ϕ₀)
    probclass = ODEProblem(classic!,init,classspan)
    solclass = solve(probclass,Tsit5())

    # α(ϕ) for ϕ>0
    solu = append!(solquant[1,:],solclass[1,:]) # α
    solt = append!(solquant.t,solclass.t) # ϕ

    # α(ϕ) for ϕ<0
    soloppu = reverse(solu)
    soloppt = -reverse(solt)
    pop!(soloppu)
    pop!(soloppt)

    # Join the two solutions
    soltotu = append!(soloppu,solu)
    soltott = append!(soloppt,solt)
    soltot = DiffEqArray(soltotu,soltott)
end

plot(alpha_of_phi!(10000.0,-0.0009,0.0074847,2.0),yaxis=:log)
plot!(alpha_of_phi!(16.0,-0.1,0.00001,2.0))
```

Answer 1

如果您正在绘制由solve 直接返回的解，则微分方程的绘图配方为plot 启用一个名为plotdensity 的可选关键字参数，它可以让您选择绘制的点数，从而选择平滑度，如docs 中所述，例如：

plot(sol,plotdensity=10000)

但是，此关键字似乎需要解决方案对象，而不是 DiffEqArray。因此，您最好的选择确实是手动设置saveat。对于这种方法，saveat = 0.01 似乎足以获得完全平滑的线条。但是，这仍然会给出您描述的“范围步长不能为零”的错误。

虽然我对您正在求解的系统没有深入了解，但对结果的检查显示，在没有 saveat 的情况下运行 alpha_of_phi!(16.0,-0.1,0.00001,2.0) 的结果中有重复的时间步长，这表明经典模拟运行的范围超过没时间。换句话说，这暗示last(solquant.t) 很可能等于或大于ϕ₀ 与这些参数，导致时间跨度为零。如果是这样，当您向saveat 请求某个有限时间在该时间跨度(last(solquant.t), ϕ₀) 时，这将完全可以理解地失败。

因此，如果我们只是重写你的函数来检查这个条件，那么在这个假设上工作

using DifferentialEquations, Plots

function alpha_of_phi!(s2,d,a0,ϕ₀)
    # α in the quantum phase
    function quantum!(dv,v,p,ϕ)
       s2,d=p
       α = v[1]
       dv[1] = (ϕ*s2*sin(2*d*α)+2*d*sinh(s2*α*ϕ))/(-α*s2*sin(2*d*α)+2*d*cos(2*d*α)+2*d*cosh(s2*α*ϕ))
    end
    # When dα/dϕ = 1, we reach the classical regime, and stop the integration
    condition(v,ϕ,integrator) = (ϕ*s2*sin(2*d*v[1])+2*d*sinh(s2*v[1]*ϕ))/(-v[1]*s2*sin(2*d*v[1])+2*d*cos(2*d*v[1])+2*d*cosh(s2*v[1]*ϕ))==1.0
    affect!(integrator) = terminate!(integrator)
    cb = DiscreteCallback(condition,affect!)
    # Initial Condition at the bounce
    α₀ = [a0]
    classspan = (0,ϕ₀)
    probquant = ODEProblem(quantum!,α₀,classspan,(s2,d))
    solquant = solve(probquant,Tsit5(),callback=cb,saveat=0.01)

    # α in the classical phase
    function classic!(du,u,p,ϕ)
       αc = u[1]
       dαc = u[2]
       du[1] = dαc
       du[2] = 3*(-dαc^3/sqrt(2)+dαc^2+dαc/sqrt(2)-1)
    end
    if last(solquant.t) < ϕ₀
       # IC retrieved from end of quantum phase
       init = [last(solquant);1.0]
       classspan = (last(solquant.t),ϕ₀)
       probclass = ODEProblem(classic!,init,classspan)
       solclass = solve(probclass,Tsit5(),saveat=0.01)

       # α(ϕ) for ϕ>0
       solu = append!(solquant[1,:],solclass[1,:]) # α
       solt = append!(solquant.t,solclass.t) # ϕ
    else
       solu = solquant[1,:] # α
       solt = solquant.t # ϕ
    end

    # α(ϕ) for ϕ<0
    soloppu = reverse(solu)
    soloppt = -reverse(solt)
    pop!(soloppu)
    pop!(soloppt)

    # Join the two solutions
    soltotu = append!(soloppu,solu)
    soltott = append!(soloppt,solt)
    soltot = DiffEqArray(soltotu,soltott)
end

plot(alpha_of_phi!(10000.0,-0.0009,0.0074847,2.0),yaxis=:log)
plot!(alpha_of_phi!(16.0,-0.1,0.00001,2.0))

那么我们似乎在做生意！