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

Question

Recently, I got started with Julia's (v1.0.3) DifferentialEquations.jl package. I tried solving a simple ODE system, with the same structure as my real model, but much smaller. Depending on the solver which I use, the example either solves or throws an error. Consider this MWE, a Chemical Engineering model of a consecutive / parallel reaction in a 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)

This works perfectly. However, if a choose a different solver, say Rosenbrock23() or Rodas4(), the ODE is not solved and I get the following error:

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

I won't paste the whole stacktrace here, since it is very long, but you can easily reproduce this by changing sol = solve(prob) into sol = solve(prob, Rosenbrock23()). It seems to me that the error occurs when the solver tries to derive Jacobians, but I have no clue why. And why does the default solver work, but others don't?

Please, could anyone tell me why this error occurs and how it can be fixed?

Should you not have `rate[2]=k_2*x[2]*x[2]`? – Lutz Lehmann Jan 29 '19 at 13:59 — Lutz Lehmann, Jan 29 '19 at 13:59

score 2 · Accepted Answer · edited Jul 19 '20 at 18:16

2

Automatic differentiation works by passing Dual types through your function, instead of the floats you would normally use it with. So the problem arises because you fix the internal value rate to be of type Vector{Float64} (see the third point here, and this advice). Fortunately, that's easy to fix (and even better looking, IMHO):

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

That works with both Rosenbrock23 and Rodas4.

Alternatively, you can turn off AD with Rosenbrock23(autodiff=false) (which, I think, will use finite differences instead), or supply a Jacobian.

edited Jul 19 '20 at 18:16

Chris Rackauckas

18,645
3
50
81

answered Jan 29 '19 at 14:15

phipsgabler

20,535
4
40
60

2

It's worth noting too that turning off autodiff, like `Rosenbrock23(autodiff=false)`, will work as well. Or supplying a Jacobian function. – Chris Rackauckas Jan 30 '19 at 08:01
Small error in your answer, since the fourth rate is actually `k_4*x[2]`, not `k_4*x[4]`, and the second rate should be, as pointed out by @LutzL, `k_2*x[2]*x[2]`. So the overall expression reads `rate = [k_1*x[1], k_2*x[2]*x[2], k_3*x[3], k_4*x[2]]`. Now, autodiff works! – winkmal Jan 30 '19 at 09:54
I have two more issues: First, I would like to include `rate` as part of the solution object, so it can be viewed same as `u` and `t`. Is this possible somehow? Second, when I use `ParameterizedFunctions.jl`, autodiff does not work either, but the error message is different. I should probably post a new question for the second issue and reference this one. – winkmal Jan 30 '19 at 09:59
Yes, please do ask a new question. – phipsgabler Jan 30 '19 at 10:03
Done, [see here](https://stackoverflow.com/q/54438451/4027729). How about Issue #1? Is it possible to store intermediate variables (other than `u` and `t`) in the `sol` object? If so, how? – winkmal Jan 30 '19 at 14:42
Maybe by using the [integrator interface](http://docs.juliadiffeq.org/stable/basics/integrator.html#Integrator-Interface-1), and maintaining a separate array. I'm actually the wrong person to ask that, though, never having used the package. But it might be a question worth asking here as well, or on Discourse. – phipsgabler Jan 30 '19 at 15:37

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

1 Answers1

Linked