How to fix TypeError: in setindex! in DifferentialEquations.jl

Question

最近，我开始使用 Julia 的 (v1.0.3) DifferentialEquations.jl 软件包。我尝试求解一个简单的 ODE 系统，其结构与我的真实模型相同，但要小得多。根据我使用的求解器，该示例要么解决要么抛出错误。考虑这个 MWE，CSTR 中连续/平行反应的化学工程模型：

using DifferentialEquations
using Plots

# Modeling a consecutive / parallel reaction in a CSTR
# A --> 2B --> C, C --> 2B, B --> D
# PETERSEN-Matrix
#   No.     A       B       C       D       Rate
#   1      -1       2                       k1*A
#   2              -2       1               k2*B*B
#   3               2      -1               k3*C
#   4              -1               1       k4*B

function fpr(dx, x, params, t)
    k_1, k_2, k_3, k_4, q_in, V_liq, A_in, B_in, C_in, D_in = params
    # Rate equations
    rate = Array{Float64}(undef, 4)
    rate[1] = k_1*x[1]
    rate[2] = k_2*x[2]*x[2]
    rate[3] = k_3*x[3]
    rate[4] = k_4*x[2]

    dx[1] = -rate[1] + q_in/V_liq*(A_in - x[1])
    dx[2] = 2*rate[1] - 2*rate[2] + 2*rate[3] - rate[4] + q_in/V_liq*(B_in - x[2])
    dx[3] = rate[2] - rate[3] + q_in/V_liq*(C_in - x[3])
    dx[4] = rate[4] + q_in/V_liq*(D_in - x[4])
end 

u0 = [1.5, 0.1, 0, 0]
params = [1.0, 1.5, 0.75, 0.15, 3, 15, 0.5, 0, 0, 0]
tspan = (0.0, 15.0)
prob = ODEProblem(fpr, u0, tspan, params)
sol = solve(prob)
plot(sol)

这非常有效。但是，如果选择不同的求解器，例如 Rosenbrock23() 或 Rodas4()，则 ODE 未求解并且出现以下错误：

ERROR: LoadError: TypeError: in setindex!, in typeassert, expected Float64,
got ForwardDiff.Dual{Nothing,Float64,4}

我不会在这里粘贴整个堆栈跟踪，因为它很长，但您可以通过将 sol = solve(prob) 更改为 sol = solve(prob, Rosenbrock23()) 轻松重现它。在我看来，错误发生在求解器试图推导雅可比矩阵时，但我不知道为什么。为什么默认求解器可以工作，而其他求解器不能？

拜托，谁能告诉我为什么会出现此错误以及如何修复它？

Answer 1

自动微分的工作原理是将 Dual 类型传递给您的函数，而不是您通常使用的浮点数。所以问题就出现了，因为你将内部值rate固定为Vector{Float64}类型（见第三点here, and this advice）。幸运的是，这很容易修复（恕我直言，甚至更好看）：

julia> function fpr(dx, x, params, t)
           k_1, k_2, k_3, k_4, q_in, V_liq, A_in, B_in, C_in, D_in = params
           # Rate equations
           # should actually be rate = [k_1*x[1], k_2*x[2]*x[2], k_3*x[3], k_4*x[2]], as per @LutzL's comment
           rate = [k_1*x[1], k_2*x[2], k_3*x[3], k_4*x[2]]

           dx[1] = -rate[1] + q_in/V_liq*(A_in - x[1])
           dx[2] = 2*rate[1] - 2*rate[2] + 2*rate[3] - rate[4] + q_in/V_liq*(B_in - x[2])
           dx[3] = rate[2] - rate[3] + q_in/V_liq*(C_in - x[3])
           dx[4] = rate[4] + q_in/V_liq*(D_in - x[4])
       end

这适用于 Rosenbrock23 和 Rodas4。

或者，您可以使用 Rosenbrock23(autodiff=false) 关闭 AD（我认为这将使用有限差分），或者提供雅可比行列式。

How to fix TypeError: in setindex! in DifferentialEquations.jl

How to fix TypeError: in setindex! in DifferentialEquations.jl

typeerror

ode

differential-equations

julia