High order differentiation on rational equations

[Jerrychow0505]

I'm working on an ODE system, and using the function Symbolics.expand_derivatives. This works fine when the system does not include rational equations. However, when the system is like:

@parameters t k1 k2 k3 k4
@variables r(t) w(t)
D = Differential(t)

eqs = [D(r) ~ k1*r*w/(k2+w) - k3*r,
       D(w) ~ (-k1/k4)*(k1*r*w/(k2+w))]

Symbolics.expand_derivatives only works when differentiation is under 3rd derivatives. When the order is higher, like 4. It returns the error:

ERROR: BoundsError: attempt to access 2-element Vector{Any} at index [3]
Stacktrace:
 [1] getindex(A::Vector{Any}, i1::Int64)
   @ Base ./array.jl:861
 [2] expand_derivatives(O::Term{Real, Nothing}, simplify::Bool; occurances::Term{Real, Nothing}) (repeats 5 times)
   @ Symbolics ~/.julia/packages/Symbolics/jplLc/src/diff.jl:215
 [3] expand_derivatives (repeats 2 times)
   @ ~/.julia/packages/Symbolics/jplLc/src/diff.jl:142 [inlined]

The problematic equation is:

Differential(t)((-3(k1^2)*Differential(t)(w(t))*Differential(t)(Differential(t)(r(t))) - (k1^2)*w(t)*Differential(t)(Differential(t)(Differential(t)(r(t)))) - (k1^2)*r(t)*Differential(t)(Differential(t)(Differential(t)(w(t)))) - 3(k1^2)*Differential(t)(r(t))*Differential(t)(Differential(t)(w(t)))) / (k4*(k2 + w(t))) + (-(k1^2)*(2k2 + 2w(t))*(Differential(t)(w(t))^3)*r(t) - 2(k1^2)*(Differential(t)(w(t))^3)*r(t)*w(t) - (k1^2)*(2k2 + 2w(t))*(Differential(t)(w(t))^2)*w(t)*Differential(t)(r(t)) - 2(k1^2)*(2k2 + 2w(t))*r(t)*w(t)*Differential(t)(w(t))*Differential(t)(Differential(t)(w(t)))) / (k4*((k2 + w(t))^4)) + (2(k1^2)*(Differential(t)(w(t))^2)*Differential(t)(r(t)) + (k1^2)*w(t)*Differential(t)(w(t))*Differential(t)(Differential(t)(r(t))) + (k1^2)*r(t)*w(t)*Differential(t)(Differential(t)(Differential(t)(w(t)))) + 2(k1^2)*w(t)*Differential(t)(r(t))*Differential(t)(Differential(t)(w(t))) + 3(k1^2)*r(t)*Differential(t)(w(t))*Differential(t)(Differential(t)(w(t)))) / (k4*((k2 + w(t))^2)) + (((k1^2)*r(t)*Differential(t)(w(t)) + (k1^2)*w(t)*Differential(t)(r(t)))*Differential(t)(Differential(t)(w(t)))) / (k4*((k2 + w(t))^2)) - k4*((-(k1^2)*w(t)*Differential(t)(Differential(t)(r(t))) - (k1^2)*r(t)*Differential(t)(Differential(t)(w(t))) - 2(k1^2)*Differential(t)(r(t))*Differential(t)(w(t))) / ((k4^2)*((k2 + w(t))^2)))*Differential(t)(w(t)) - 4k4*((k2 + w(t))^3)*((-(k1^2)*(2k2 + 2w(t))*(Differential(t)(w(t))^2)*r(t)*w(t)) / ((k4^2)*((k2 + w(t))^8)))*Differential(t)(w(t)) - k4*(((k1^2)*(Differential(t)(w(t))^2)*r(t) + (k1^2)*w(t)*Differential(t)(r(t))*Differential(t)(w(t)) + (k1^2)*r(t)*w(t)*Differential(t)(Differential(t)(w(t)))) / ((k4^2)*((k2 + w(t))^4)))*(2k2 + 2w(t))*Differential(t)(w(t)) - k4*(2k2 + 2w(t))*((((k1^2)*r(t)*Differential(t)(w(t)) + (k1^2)*w(t)*Differential(t)(r(t)))*Differential(t)(w(t))) / ((k4^2)*((k2 + w(t))^4)))*Differential(t)(w(t)))

The code I use to generate the derivatives is:

using DifferentialEquations
using Symbolics


# create the differentiation equations 
function derivat(e, H)
    tmp = e
    for h in 1:H
        for i in 1:length(tmp)
            println(D(tmp[i].rhs))
            println()
            tmp[i] = D(tmp[i].lhs) ~ Symbolics.expand_derivatives(D(tmp[i].rhs))     
        end
    end
    return eq
end 


@parameters t k1 k2 k3 k4
@variables r(t) w(t)
D = Differential(t)

eqs = [D(r) ~ k1*r*w/(k2+w) - k3*r,
       D(w) ~ (-k1/k4)*(k1*r*w/(k2+w))]

h = 4
eq1 = derivat(eqs, h)

---

Really don't know what happens here, Hope this can be solved.

[orebas]

I am getting this error as well. It may be related to issue 1126.