如何解决Rosenbrock23与ParameterizedFunctions.jl DSL的ODE兼容报错？

I've run into similar compatibility issues with ParameterizedFunctions and rigid solvers like Rosenbrock23 before. The root problem is that the old @ode_def DSL generates code that doesn't play nicely with the type expectations of Rosenbrock23 (even with autodiff=false). Here are two reliable fixes:

Option 1: Migrate to ModelingToolkit (Recommended) #

ModelingToolkit is the modern replacement for ParameterizedFunctions in the DifferentialEquations ecosystem—it's more flexible, type-stable, and works seamlessly with rigid solvers. Here's your model rewritten with MTK:

using DifferentialEquations
using Plots
using ModelingToolkit

# Define the consecutive/parallel CSTR reaction system symbolically
@variables t A(t) B(t) C(t) D(t)
@parameters k_1 k_2 k_3 k_4 q_in V_liq A_in B_in C_in D_in

D = Differential(t)

# Write out the ODE equations
eqs = [
    D(A) ~ -k_1*A + q_in/V_liq*(A_in - A),
    D(B) ~ 2*k_1*A - 2*k_2*B^2 + 2*k_3*C - k_4*B + q_in/V_liq*(B_in - B),
    D(C) ~ k_2*B^2 - k_3*C + q_in/V_liq*(C_in - C),
    D(D) ~ k_4*B + q_in/V_liq*(D_in - D)
]

# Build and simplify the system
@named sys = ODESystem(eqs)
sys = structural_simplify(sys)

# Set up initial conditions, parameters, and time span
u0 = [
    A => 1.5,
    B => 0.1,
    C => 0.0,
    D => 0.0
]

params = [
    k_1 => 1.0,
    k_2 => 1.5,
    k_3 => 0.75,
    k_4 => 0.15,
    q_in => 3.0,
    V_liq => 15.0,
    A_in => 0.5,
    B_in => 0.0,
    C_in => 0.0,
    D_in => 0.0
]

tspan = (0.0, 15.0)

# Solve and plot
prob = ODEProblem(sys, u0, tspan, params)
sol = solve(prob, Rosenbrock23(), autodiff=false)
plot(sol)

Option 2: Manually Define the ODE Function #

If you prefer to stick with a non-symbolic approach, writing the ODE function manually avoids the type issues introduced by the ParameterizedFunctions DSL:

using DifferentialEquations
using Plots

# Manually implement the CSTR ODE function
function fpr!(du, u, p, t)
    # Unpack parameters and state variables
    k_1, k_2, k_3, k_4, q_in, V_liq, A_in, B_in, C_in, D_in = p
    A, B, C, D = u

    # Compute derivatives
    du[1] = -k_1*A + q_in/V_liq*(A_in - A)
    du[2] = 2*k_1*A - 2*k_2*B*B + 2*k_3*C - k_4*B + q_in/V_liq*(B_in - B)
    du[3] = k_2*B*B - k_3*C + q_in/V_liq*(C_in - C)
    du[4] = k_4*B + q_in/V_liq*(D_in - D)
end

# Setup and solve
u0 = [1.5, 0.1, 0.0, 0.0]
params = [1.0, 1.5, 0.75, 0.15, 3.0, 15.0, 0.5, 0.0, 0.0, 0.0]
tspan = (0.0, 15.0)

prob = ODEProblem(fpr!, u0, tspan, params)
sol = solve(prob, Rosenbrock23(), autodiff=false)
plot(sol)

Why the Original Code Failed #

The @ode_def macro generates code under the hood that doesn't handle parameter passing in a way that Rosenbrock23 expects—even with autodiff=false, it tries to convert an array to a scalar at some point during solver setup or derivative calculations. ParameterizedFunctions is essentially deprecated at this point, so switching to ModelingToolkit or using a manual function is the best long-term fix.

内容的提问来源于stack exchange，提问作者winkmal