Как исправить ошибку TypeError: в setindex! в DifferentialEquations.jl

Question

Недавно я начал с Джулии (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} (см. третий пункт здесь, и этот совет). К счастью, это легко исправить (и даже лучше, ИМХО):

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,

Кроме того, вы можете отключить AD с Rosenbrock23(autodiff=false) (который, я думаю, вместо этого будет использовать конечные различия), или поставит якобиана.